Lecturer: Brad Cenko
Jupyter Notebook Author: Brad Cenko, Dipankar Bhattacharya & Cameron Hummels
This is a Jupyter notebook lesson taken from the GROWTH Summer School 2020. For other lessons and their accompanying lectures, please see: http://growth.caltech.edu/growth-school-2020.html
Learn how to analyze x-ray data
See GROWTH school webpage for detailed instructions on how to install these modules and packages. Nominally, you should be able to install the python modules with pip install <module>
. The external astromatic packages are easiest installed using package managers (e.g., rpm
, apt-get
).
None.
from astropy.io import fits
import numpy as np
import matplotlib.pyplot as plt
import math
import scipy as sc
import os
Data from X-ray instruments are typically stored as event files - binary FITS tables with a list of photons detected (including time, location on the detector, and photon energy). Here we will start off by reading in one of these event files (in this case referred to as a "Level 2" event file because some cleaning has been done to filter for astrophysical photons) from the Cadmium Zinc Telluride Imager (CZTI) on the AstroSat satellite.
filename = os.path.join('data', 'AS1A02_005T01_9000000948_06884cztM0_level2_common_clean.evt')
dataHDU = fits.open(filename)
Here's what one of these event files looks like.
dataHDU[1].header
You can see that this was an observation of the bright AGN Mrk421 ("OBJECT" keyword) obtained on January 5 2017 ("DATE-OBS" keyword), with a total elapsed time of 6975 s ("TELAPSE" keyword) and an exposure time of 3943 s ("EXPOSURE" keyword). Now let's look at the actual data in the binary table. The attribute "dtype" describes the columns in a FITS record array.
dataHDU[1].data.dtype
CZTI contains four quadrants, here we obtain the time value for each quadrant and combine the data for them. Since the FITS tables are a list of detected photons, combining the list of times effectively produces a list of the times of all photons observed by the detector.
times = np.concatenate((dataHDU[1].data['Time'], dataHDU[2].data['Time'],
dataHDU[3].data['Time'], dataHDU[4].data['Time']))
To generate a light curve, we can just create a histogram of the times of photon arrival from the four quadrants. Defining the bin size (in this case 5 s) will significantly impact the appearance of the light curve.
# Define timebins
binsize = 5
tbins = np.arange(times.min(), times.max(), binsize)
# Make histogoram
counts, bins = np.histogram(times, bins=tbins)
bins = (bins[1:] + bins[:-1])/2
# Plot
plt.plot(bins, counts/binsize, drawstyle='steps-mid')
plt.xlabel('Time (s)')
plt.ylabel('$counts\ s^{-1}$')
plt.show()
Student Exercises: From other high-energy detectors, we know that a GRB occured at UTC 2017 Jan 5 06:14:06 (GRB170105A), corresponding to a mission time of ~ 221292850.
Part 1: Did the CZTI see a GRB at this time? Plot the four-quadrant light curve around this time window with a variety of different bin sizes to see if there is any evidence for a GRB at this time.
Part 2: Estimate the duration of GRB170105A in the CZTI bandpass (i.e., how long was there signal above the background level)?
We will look at observations obtained by the Large Area X-ray Proportional Counter (LAXPC) instrument on AstroSat to see if we can measure pulsations from the Crab pulsar. We need following files for this tutorial:
The files can be open in the same format as described in the example below.
GTI = Good Time Interval
fevents=fits.open("data/ObsID406_02741_event.fits")
fevents_bary=fits.open("data/ObsID406_02741_laxpc_bary.fits")
fgti=fits.open('data/ObsID406_02741_laxpc1_bary.gti')
From the files we can obtain informations like the length, header etc.
fevents_bary[1].header
For accurate timing analysis, we need to put all the observations times onto a common reference system. This is typically referenced to the frame of the Sun / Solar System, and is called a barycenter correction. Here the correction has been applied for us by the AstroSat pipeline. Observe the time difference between the data from the event file and barycenter correction file.
time=fevents[1].data['TIME']
time_bary=fevents_bary[1].data['TIME']
time_diff=time[0]- time_bary[0]
time_diff
In addition to the barycenter correction, we need to account for the fact that only a fraction of the data is obtained during so-called "Good Time Intervals (GTIs)". This could be due to time periods of elevated background (e.g., passage through the South Atlantic Anomaly), or simply because the target location is occulted by Earth.
First lets plot the light curve without applying gti cuts. This plot will consist of both gti(good time interval) data as well as bti (bad time interval) data. We can compare the two light curves obtained before and after applying gti cuts.
# Define timebins
binsize = 1
tbins = np.arange(time_bary.min(), time_bary.max(), binsize)
# Make histogoram
counts_time, t_bins = np.histogram(time_bary, bins=tbins)
t_bins = (t_bins[1:] + t_bins[:-1])/2
# Plot
plt.plot(t_bins, counts_time/binsize, drawstyle='steps-mid')
plt.xlabel('Time (s)')
plt.ylabel('$counts\ s^{-1}$')
plt.show()
Read and filter event file to keep in array only events within the good time interval. You can see from the plot above that there are two good time intervals in the above data set.
gtidata=fgti[1].data
gtidata
t_start=gtidata[0][0],gtidata[1][0]
t_stop=gtidata[0][1],gtidata[1][1]
Student Exercise: Using the GTIs defined above, plot the light curve in the second good time interval.
While it is not apparent in the plot above, the X-ray light curves includes a periodic modulation due to the Crab pulsar.
Student Exercise: Calculate the power spectrum of the light curve in the second good time interval, and search for periodicity due to the Crab pulsar (hint: the period is ~ 30 Hz)
With the period in hand, we can calculate a phase folded light curve to measure the pulse profile. We can use the period measured above (P = 29.655 Hz) to perform the phase folding.
Student exercise: Fold the light curve above at the period of the Crab pulsar to determine the pulse profile. Hint: use the following period: 29.6553306828504