Lecturer: Brad Cenko
Jupyter Notebook Author: Brad Cenko, Dipankar Bhattacharya & Cameron Hummels
This is a Jupyter notebook lesson taken from the GROWTH Summer School 2019. For other lessons and their accompanying lectures, please see: http://growth.caltech.edu/growth-school-2019.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
XTENSION= 'BINTABLE' / binary table extension BITPIX = 8 / 8-bit bytes NAXIS = 2 / 2-dimensional binary table NAXIS1 = 33 / width of table in bytes NAXIS2 = 477590 / number of rows in table PCOUNT = 0 / size of special data area GCOUNT = 1 / one data group (required keyword) TFIELDS = 12 / number of fields in each row TTYPE1 = 'Time ' / label for field 1 TFORM1 = 'D ' / data format of field: 8-byte DOUBLE TTYPE2 = 'CZTSECCNT' / label for field 2 TFORM2 = 'D ' / data format of field: 8-byte DOUBLE TTYPE3 = 'CZTNTICK' / label for field 3 TFORM3 = 'I ' / data format of field: 2-byte INTEGER TZERO3 = 32768 / offset for unsigned integers TSCAL3 = 1 / data are not scaled TTYPE4 = 'PHA ' / label for field 4 TFORM4 = 'I ' / data format of field: 2-byte INTEGER TZERO4 = 32768 / offset for unsigned integers TSCAL4 = 1 / data are not scaled TTYPE5 = 'DetID ' / label for field 5 TFORM5 = 'B ' / data format of field: BYTE TTYPE6 = 'pixID ' / label for field 6 TFORM6 = 'B ' / data format of field: BYTE TTYPE7 = 'DETX ' / label for field 7 TFORM7 = 'B ' / data format of field: BYTE TTYPE8 = 'DETY ' / label for field 8 TFORM8 = 'B ' / data format of field: BYTE TTYPE9 = 'veto ' / label for field 9 TFORM9 = 'I ' / data format of field: SHORT INTEGER TZERO9 = 32768 / offset for unsigned integers TSCAL9 = 1 / data are not scaled TTYPE10 = 'alpha ' / label for field 10 TFORM10 = 'B ' / data format of field: BYTE TUNIT1 = 's ' / physical unit of field 1 TUNIT2 = 's ' / physical unit of field 2 TUNIT3 = 'micro_sec' / physical unit of field 3 TUNIT4 = 'counts ' / physical unit of field 4 TUNIT9 = 'counts ' / physical unit of field TUNIT10 = 'count ' / physical unit of field EXTNAME = 'Q0 ' / name of this binary table extension QUADID = 0 / Quadrant Number TSTARTI = 221288272 / Start time of observation Integer part TSTARTF = 0. / Start time of observation Fractional part TSTOPI = 221295247 / Stop time of observation Integer part TSTOPF = 0. / Stop time of observation Fractional part EXTVER = 1 / auto assigned by template parser EXPOSURE= 3942.84576797485 / Exposure time MISSION = 'ASTROSAT' / Name of the mission/satellite TELESCOP= 'ASTROSAT' / Name of the mission/satellite INSTRUME= 'CZTI ' / Name of the instrument/detector ORIGIN = 'CZTI POC' / Source of FITS FILE TIMESYS = 'UTC ' / Time is UTC TIMEUNIT= 's ' / Time is in seconds MJDREFI = 55197 / MJDREF Integer part MJDREFF = 0 / MJDREF Fractional part TSTART = 221288272. / Start time of observation TSTOP = 221295247. / Stop time of observation OBJECT = 'Mrk421 ' / Target name RADECSYS= 'ICRS ' / Reference frame EQUINOX = 2000 / J2000 RA_PNT = 166.1138 / Nominal Pointing RA DEC_PNT = 38.20883 / Nominal Pointing DEC OBS_ID = 'A02_005T01_9000000948' / Observation ID OBS_MODE= 'POINTING' DATE-OBS= '2017-01-05' / Start date of observation TIME-OBS= '04:57:51.739900000' / Start time of observation DATE-END= '2017-01-05' / End date of observation TIME-END= '06:54:09.740600000' / End time of observation TIMEDEL = 2.00000027772809E-05 / Time resolution TELAPSE = 6975. / Elapsed time HISTORY Module run by czti HISTORY Parameter List START for cztscience2event_v1.0 HISTORY P1 infile=/data2/czti/level1/20170103_A02_005T01_9000000948_level1/czti/ HISTORY orbit/06884/modeM0/AS1A02_005T01_9000000948cztM0_level1.fits HISTORY P2 TCTfile=/data2/czti/level1/20170103_A02_005T01_9000000948_level1/czti HISTORY /orbit/06884/aux/AS1A02_005T01_9000000948czt_level1.tct HISTORY P3 outfile=/data2/czti/level2/20170103_A02_005T01_9000000948_level2/czti HISTORY /orbit/06884/modeM0/AS1A02_005T01_9000000948_06884cztM0_level2.fits HISTORY P4 hdrInfoFile=/data2/czti/level2/20170103_A02_005T01_9000000948_level2/ HISTORY czti/orbit/06884/modeM0/AS1A02_005T01_9000000948_06884cztM0_level2.hdr HISTORY P7 BigEndian=yes HISTORY P8 clobber=yes HISTORY P9 history=yes HISTORY P10 debug=yes HISTORY Parameter List END DATE = '2017-01-05T08:08:21' / file creation date (YYYY-MM-DDThh:mm:ss UT) CREATOR = 'cztevtclean_v1.0' / Module that created this file FILENAME= '/data2/czti/level2/20170103_A02_005T01_9000000948_level2/czti/orbit/' CHECKSUM= '6AKA63K369K969K9' / HDU checksum updated 2017-01-05T08:08:21 DATASUM = '2829133320' / data unit checksum updated 2017-01-05T08:08:21 TTYPE11 = 'PI ' / label for field TFORM11 = 'I ' / format of field TTYPE12 = 'ENERGY ' / label for field TFORM12 = 'E ' / format of field HISTORY Module run by czti HISTORY Parameter List START for cztpha2energy_v1.0 HISTORY P1 infile=/data2/czti/level2/20170103_A02_005T01_9000000948_level2/czti/ HISTORY orbit/06884/modeM0/AS1A02_005T01_9000000948_06884cztM0_level2.fits HISTORY P2 Ebounds file=/home/czti/new_CALDB/data/as1/czti/bcf/AS1cztebounds2015 HISTORY 0707v01.fits HISTORY P3 Gain file=/home/czti/new_CALDB/data/as1/czti/bcf/AS1cztgain20150707v0 HISTORY 1.fits HISTORY P5 outfile file=/data2/czti/level2/20170103_A02_005T01_9000000948_level2 HISTORY /czti/orbit/06884/modeM0/AS1A02_005T01_9000000948_06884cztM0_level2.evt HISTORY P5 clobber=yes HISTORY P6 history=yes HISTORY Parameter List END HISTORY Module run by czti HISTORY Parameter List START for cztbunchclean HISTORY P1 infile=/data2/czti/level2/20170103_A02_005T01_9000000948_level2/czti/ HISTORY orbit/06884/modeM0/AS1A02_005T01_9000000948_06884cztM0_level2.evt HISTORY P2 outfile=/data2/czti/level2/20170103_A02_005T01_9000000948_level2/czti HISTORY /orbit/06884/modeM0/AS1A02_005T01_9000000948_06884cztM0_level2_bc.evt HISTORY P3 bunchdeftime=30 HISTORY P4 bunch_length_thresh=20 HISTORY P5 skipT1=0.000000 HISTORY P6 skipT2=0.001000 HISTORY P7 skipT3=0.001000 HISTORY P8 clobber=yes HISTORY P9 history=yes HISTORY Parameter List END GTITYPE = 'COMMON ' HISTORY Module run by czti HISTORY Parameter List START for cztdatasel_v1.0 HISTORY P1 infile=/data2/czti/level2/20170103_A02_005T01_9000000948_level2/czti/ HISTORY orbit/06884/modeM0/AS1A02_005T01_9000000948_06884cztM0_level2_bc.evt HISTORY P2 gtifile=/data2/czti/level2/20170103_A02_005T01_9000000948_level2/czti HISTORY /orbit/06884/modeM0/AS1A02_005T01_9000000948_06884cztM0_level2.gti HISTORY P3 outfile=/data2/czti/level2/20170103_A02_005T01_9000000948_level2/czti HISTORY /orbit/06884/modeM0/AS1A02_005T01_9000000948_06884cztM0_level2_common_bc HISTORY _ds.evt HISTORY P4 gtitype=COMMON HISTORY P5 clobber=yes HISTORY P6 history=yes HISTORY Parameter List END HISTORY Module run by czti HISTORY Parameter List START for cztevtclean_v1.0 HISTORY P1 infile=/data2/czti/level2/20170103_A02_005T01_9000000948_level2/czti/ HISTORY orbit/06884/modeM0/AS1A02_005T01_9000000948_06884cztM0_level2_common_bc_ HISTORY ds_pc.evt HISTORY P2 outfile=/data2/czti/level2/20170103_A02_005T01_9000000948_level2/czti HISTORY /orbit/06884/modeM0/AS1A02_005T01_9000000948_06884cztM0_level2_common_cl HISTORY ean.evt HISTORY P3 alphaval=0 HISTORY P4 vetorange=0 HISTORY P5 clobber=yes HISTORY P6 history=yes HISTORY Parameter List END
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
dtype((numpy.record, [('Time', '>f8'), ('CZTSECCNT', '>f8'), ('CZTNTICK', '>i2'), ('PHA', '>i2'), ('DetID', 'u1'), ('pixID', 'u1'), ('DETX', 'u1'), ('DETY', 'u1'), ('veto', '>i2'), ('alpha', 'u1'), ('PI', '>i2'), ('ENERGY', '>f4')]))
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, ls='steps-mid')
plt.xlabel('Time (s)')
plt.ylabel('$counts\ s^{-1}$')
plt.show()
/home/growth/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:10: MatplotlibDeprecationWarning: Passing the drawstyle with the linestyle as a single string is deprecated since Matplotlib 3.1 and support will be removed in 3.3; please pass the drawstyle separately using the drawstyle keyword argument to Line2D or set_drawstyle() method (or ds/set_ds()). # Remove the CWD from sys.path while we load stuff.
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.
binsize =0.5
tbins = np.arange(221292820.0, 221292880.0, binsize)
counts, bins = np.histogram(times, bins=tbins)
bins = (bins[1:] + bins[:-1])/2
plt.plot(bins, counts/binsize, ls='steps-mid')
plt.xlabel('Time (s)')
plt.ylabel('$counts\ s^{-1}$')
plt.show()
/home/growth/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:8: MatplotlibDeprecationWarning: Passing the drawstyle with the linestyle as a single string is deprecated since Matplotlib 3.1 and support will be removed in 3.3; please pass the drawstyle separately using the drawstyle keyword argument to Line2D or set_drawstyle() method (or ds/set_ds()).
Part 2: Estimate the duration of GRB170105A in the CZTI bandpass (i.e., how long was there signal above the background level)?
binsize =1.0
tbins = np.arange(221292840.0, 221292860.0, binsize)
counts, bins = np.histogram(times, bins=tbins)
bins = (bins[1:] + bins[:-1])/2
plt.plot(bins, counts/binsize, ls='steps-mid')
plt.plot([221292840.0, 221292860.0], [550.0, 550.0], "k")
plt.xlabel('Time (s)')
plt.ylabel('$counts\ s^{-1}$')
plt.show()
# Approximately 4 bins above background, so ~ 4 s.
/home/growth/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:8: MatplotlibDeprecationWarning: Passing the drawstyle with the linestyle as a single string is deprecated since Matplotlib 3.1 and support will be removed in 3.3; please pass the drawstyle separately using the drawstyle keyword argument to Line2D or set_drawstyle() method (or ds/set_ds()).
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
XTENSION= 'BINTABLE' / binary table extension BITPIX = 8 / 8-bit bytes NAXIS = 2 / 2-dimensional binary table NAXIS1 = 18 / width of table in bytes NAXIS2 = 13614384 / number of rows in table PCOUNT = 0 / size of special data area GCOUNT = 1 / one data group (required keyword) TFIELDS = 5 / number of fields in each row TTYPE1 = 'TIME ' / Time elapsed since MJDREF TFORM1 = '1D ' / data format of field: 8-byte DOUBLE TUNIT1 = 's ' / physical unit of field TTYPE2 = 'Channel ' / label for field 2 TFORM2 = 'I2 ' / data format of field: 2-byte INTEGER TTYPE3 = 'Layer ' / label for field 3 TFORM3 = 'I2 ' / data format of field: 2-byte INTEGER TTYPE4 = 'LAXPC_No.' / label for field 4 TFORM4 = 'I2 ' / data format of field: 2-byte INTEGER TTYPE5 = 'Energy ' / label for field 5 TFORM5 = '1E3.2 ' / data format of field: 4-byte REAL TUNIT5 = 'keV ' / physical unit of field EXTNAME = 'event file' / name of this binary table extension RESPFIL1= 'lx10_17aug16.rmf' / Response file for LAXPC 10 RESPFIL2= 'lx20_17aug16.rmf' / Response file for LAXPC 20 MJDREFI = 55197 / TDB time reference; Modified Julian Day (int) MJDREFF = 0 / TDB time reference; Modified Julian Day (frac) TSTOP = 1.971066916643481E+08 / Elapsed seconds since MJDREF at end of file TSTART = 1.970990932256668E+08 / Elapsed seconds since MJDREF at start of file RESPFIL3= 'lx30_17aug16.rmf' / Response file for LAXPC 30 TIMEDEL = 1.00E-05 MISSION = 'ASTROSAT' TELESCOP= 'ASTROSAT' INSTRUME= 'LAXPC1 ' RADECSYS= 'FK5 ' / Coordinate Reference System TIMESYS = 'TDB ' / All times in this file are TDB HISTORY File modified by user 'dipankar' with fv on 2017-01-17T13:30:48 HISTORY File modified by user 'dipankar' with fv on 2017-01-17T13:44:59 DATE-OBS= '2016-03-31T05:45:59' / Date and time (TIMESYS) at start of file DATE-END= '2016-03-31T07:52:37' / Date and time (TIMESYS) at end of file TIMEZERO= 0.000000E+00 / & RA_OBJ = 8.36330800E+01 / Right Ascension used for barycenter corrections DEC_OBJ = 2.20144600E+01 / Declination used for barycenter corrections TIMEREF = 'SOLARSYSTEM' / Times are pathlength-corrected to barycenter TREFPOS = 'BARYCENTER' / Time reference position TREFDIR = 'RA_OBJ,DEC_OBJ' / Keywords of reference direction PLEPHEM = 'JPL-DE200' / Solar system ephemeris used for baryctr corr. TIERRELA= 1.000000E-09 / Short-term clock stability TIERABSO= 4.000000E+00 / Absolute precision of clock correction CREATOR = 'axBary - $Revision: 1.10 $ $Date: 2013/12/02 20:28:26 $' / & DATE = '2017-01-28T13:41:47' / Date and time (UTC) of file creation CHECKSUM= '31Ge509e30Ee309e' / HDU checksum updated 2017-01-28T13:41:47 DATASUM = '965414762' / data unit checksum updated 2017-01-28T13:41:47
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
73.24125623703003
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, ls='steps-mid')
plt.xlabel('Time (s)')
plt.ylabel('$counts\ s^{-1}$')
plt.show()
/home/growth/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:10: MatplotlibDeprecationWarning: Passing the drawstyle with the linestyle as a single string is deprecated since Matplotlib 3.1 and support will be removed in 3.3; please pass the drawstyle separately using the drawstyle keyword argument to Line2D or set_drawstyle() method (or ds/set_ds()). # Remove the CWD from sys.path while we load stuff.
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
FITS_rec([(1.97099975e+08, 1.97103300e+08), (1.97105810e+08, 1.97106645e+08)], dtype=(numpy.record, [('START', '>f8'), ('STOP', '>f8')]))
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.
eventsdata=fevents_bary[1].data
time_gti = time_bary[np.where((time_bary>=t_start[1]) & (time_bary<=t_stop[1]))]
# Define timebins
binsize = 1
tbins = np.arange(time_gti.min(), time_gti.max(), binsize)
# Make histogoram
counts_time, t_bins = np.histogram(time_gti, bins=tbins)
t_bins = (t_bins[1:] + t_bins[:-1])/2
# Plot
plt.plot(t_bins, counts_time/binsize, ls='steps-mid')
plt.xlabel('Time (s)')
plt.ylabel('$counts\ s^{-1}$')
plt.show()
/home/growth/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:13: MatplotlibDeprecationWarning: Passing the drawstyle with the linestyle as a single string is deprecated since Matplotlib 3.1 and support will be removed in 3.3; please pass the drawstyle separately using the drawstyle keyword argument to Line2D or set_drawstyle() method (or ds/set_ds()). del sys.path[0]
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)
binsize = 0.001
tbins = np.arange(time_gti.min(), time_gti.max(), binsize)
counts_time, t_bins = np.histogram(time_gti, bins=tbins)
t_bins = (t_bins[1:] + t_bins[:-1])/2
import numpy.fft as fft
n = len(counts_time)
T=time_gti.max()-time_gti.min()
frq = np.arange(n)/T
power = fft.fft(counts_time/binsize)
plt.plot(frq, abs(power))
plt.xlim(29.0, 31.0)
plt.ylim(0.0, 0.1e9)
plt.show()
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
tref=fevents_bary[1].header['MJDREFI']
tref=tref*86400
time_current =time_gti + tref
T0=57480.5635028439*86400 # Reference Epoch for which nu is known
t_phase=time_current-T0
# freq
nu = 29.6553306828504
#calculating phase
phi=t_phase*nu
phi_int=np.zeros(len(phi))
for i in range(len(phi)):
phi_int[i]= int(phi[i])
phase=phi - phi_int
phase=phase+1
min(phase),max(phase)
#folding phase
no_bins=250
pbins=np.linspace(0.,1.,no_bins+1)
phase_counts, pbin_edges = np.histogram(phase,bins=pbins)
pbin_edges = (pbin_edges[1:] + pbin_edges[:-1])/2
phase_counts=[float(i) for i in phase_counts]
phase_counts=np.array(phase_counts)
pcounts=phase_counts/max(phase_counts) # normalising
plt.plot(pbin_edges, pcounts) # plotting pulse profile
plt.xlabel('Phase')
plt.ylabel('norm ')
plt.show()