Chandra Analysis

Let's log in to calgary and start the software we need

Log in:

ssh -CXY username@calgary.phy.bris.ac.uk

Make a directory to work in

cd /data/cluster2
mkdir username
cd username

So now I am working in /data/cluster2/bjm and this is where I will keep all my project work. From now on when I log in, I'll just go right there using

cd /data/cluster2/bjm

Now start the HEASoft and CIAO software

setenv HEADAS /usr/local/heasoft-6.26/x86_64-pc-linux-gnu-libc2.12
source $HEADAS/headas-init.csh
source /usr/local/ciao64/ciao-4.12/bin/ciao.csh

download_chandra_obsid 520

The data are now stored in a directory called 520.

Reprocess the data to match the latest calibration files

cd 520
chandra_repro indir="./" outdir=repro

This takes about 10 minutes, and creates a directory called repro within which is the new calibrate "level 2" events file called (in this case) acisf00520_repro_evt2.fits. Remind yourself of what an events file is by reviewing my X-ray astronomy lectures

Let's take a quick look at this even file using ds9

cd repro
ds9 acisf00520_repro_evt2.fits &

We should set the preferences in the ds9 edit menu to make it easier to see the data. Go to edit -> preferences -> Menus and Buttons and set Scale -> Menu -> log to set the colour scale to be logarithmic with image intensity, and Scale -> Menu -> Color -> bb for a nice colour scale. Now if you close and reopen ds9 you'll see the effect of the changes.

Setting the binning factor to 4 (bin -> bin 4) I see this, with the cluster clearly visible

You need to do this each time you start a new terminal session or if you switch to analysing a different cluster (more info).

The bad pixel file is called *bpix1* in the repro directory:

cd repro
ls *bpix1*
acisf00520_repro_bpix1.fits

The following command tells CIAO where to find this file

punlearn ardlib
acis_set_ardlib acisf00520_repro_bpix1.fits

Let's now clean the lightcurve of the observation.

First extract a lightcurve, working in the repro directory

dmextract "acisf00520_repro_evt2.fits[energy=300:12000][bin time=::260]" outfile=520_lc.fits opt=ltc1

Now identify good times in the observation

deflare 520_lc.fits 520_gti.fits method=clean plot=no save=520_lc.png

This creates a good time file called 520_gti.fits and a plot called 520_lc.png. We can look at the latter with

eog 520_lc.png

Apply the good time interval file to the events list to remove any bad times

dmcopy "acisf00520_repro_evt2.fits[@520_gti.fits]" 520_evt2.fits

This creates our final cleaned events list that we will use for all our analysis, 520_evt2.fits

We can make a blank sky background events file to match our observation, which can be helpful for some analysis

blanksky 520_evt2.fits 520_bg_evt2.fits

Look at this is ds9 and you should see a very empty looking observation!

We can make an image and exposure map using the fluximage command. Refer to the X-ray astronomy lectures to remind yourself what exposure maps are.

Let's organise ourselves by making a new directory for our image analysis. Starting in our 520 directory

mkdir work
cd work
mkdir img
cd img

Now make the images, in this example using the 0.5-2 keV band with a binning factor of 4 (this might not be optimal for your science)

fluximage ../../repro/520_evt2.fits outroot=520b4 binsize=4 bands="0.5:2:1.5"

Use ds9 to look at the output images.

You can make a simple smoothed image using the Analysis -> Smooth option in ds9 and use the Analysis -> Smooth Parameters to adjust the parameters. These images are just for display though, to make saved versions of the smoothed images we can use fgauss.

This command smooths the exposure-corrected image with a 4 pixel Gaussian:

fgauss 520b4_0.5-2_flux.img 520b4_0.5-2_flux_g4.img 4

Open the smoothed image 520b4_0.5-2_flux_g4.img in ds9

Make a circular region in ds9 centred roughly on the bright part of the cluster and save it in CIAO format with physical pixels. My region looks like this

and I saved it as 520_cent.reg.

Now use dmstat to find the centroid

dmstat 520b4_0.5-2_thresh.img"[sky=region(520_cent.reg)]"

EVENTS_IMAGE(x, y)
    min:	0             @:	( 4038.5 3438.5 )
    max:	23            @:	( 4074.5 3710.5 )
cntrd[log] :	( 53.38769389 51.07027166 )
cntrd[phys]:	( 4076.0507756 3634.7810866 )
sigma_cntrd:	( 65.755769518 63.018643912 )
   good:	8194
   null:	2415

This tells me that the centroid is at coordinates (4076.0507756,3634.7810866) - check this by inspecting the image in ds9.

The dmstat output also tells me the location of the peak (4074.5,3710.5), but it is usually better to find the peak in a slightly smoothed image to avoid noise fluctuations influencing the value.

You can make a radial profile of the emission by following this thread or for a quick analysis, you can make one directly in ds9 by creating annulus regions (as in the thread linked above) and then in the region editing box, choose analysis -> radial profile.

You should usually centre your profile on the peak or centroid.

You can detect sources in your image following this thread - a minimal example would be:

wavdetect infile=520b4_0.5-2_thresh.img expfile=520b4_0.5-2_thresh.expmap outfile=520_src.fits regfile=520_src.reg scellfile=520_scell.fits imagefile=520_imgfile.fits defnbkgfile=520_nbgd.fits scales="1.0 2.0 4.0 8.0 16.0 32.0" psffile=""

Load the 520_src.reg file into ds9 to see how it did.

First decide on a source and background region to use and save them to separate files. For this example I'll use 520_cent.reg that I created earlier, and the following background region

Note that here I have excluded part of one of the background regions using Property -> exclude in the region property box. I also set that region to be at the back of the set of regions by selecting it and choosing Region -> back. I'll save this in a file called 520_bg.reg

I'll make a new directory in my work directory to work in and copy over the files I need

mkdir spec
cd spec
cp ../img/520_cent.reg .
cp ../img/520_bg.reg .

Now I can extract the spectra using specextract

specextract infile="../../repro/520_evt2.fits[sky=region(520_cent.reg)]" outroot=520 bkgfile="../../repro/520_evt2.fits[sky=region(520_bg.reg)]" grouptype=NUM_CTS binspec=5 weight="no"

This extracts a source and background spectrum and makes the appropriate RMF and ARF response files. Review the X-ray astronomy lectures to remind yourself what these are.

Note that I could have used the blank-sky background here if I wanted. I have also grouped the spectrum to 5 counts per bin.

We can measure the properties of the cluster by fitting a model to the spectrum using XSpec.

First of all we need to look up the absorbing column density towards the cluster. For this we need the RA and DEC. We can read these from the header of the events file

dmkeypar ../../repro/520_evt2.fits RA_TARG echo+
dmkeypar ../../repro/520_evt2.fits DEC_TARG echo+

which give me the RA and DEC in degrees

4.641088
16.435486

Now I can look up the nH value

nh 2000 4.641088 16.435486

and we want the Weighted average value, which is 3.77E+20 in this case.

Now start Xspec

xspec

and then we type the following commands at the xspec (ignore the comment lines starting with #)

# set the nH (in units of 1e22 cm^-3 and redshift
set nh 0.0377
set z 0.541

# set up plot window
cpd /xw

# use the c statistic instead of the default chi squared
statistic cstat

#Keep going until fits converge
query yes


#make sure auto renorm is off
renorm NONE

#read data - automatically reads in background and response files
data 520_grp.pi

#ignore unwanted energies
setplot energy
ignore 0.0-0.5
ignore 7.0-**

# plot data
plot ldata

#Set up the model with rough starting values
model phabs*mekal & /*
newpar 1 $nh
newpar 2 5
newpar 4 0.3
newpar 5 $z

# unfreeze the abundance
thaw 4
# show model
newpar 0

# fit model to data
fit

# plot fit and residuals
plot ldata delchi

#make plot and remove labels and teak appearance
setplot command lab top
setplot command time off
setplot command lab f
setplot command lw 3
setplot command cs 1.25

#rebin up to 5 adjacent bins to get min sn=3 in plot ONLY - doesn't affect fit
setplot rebin 3 10

# adjust plotting range on y axis
setplot command rescale y1 1e-3 0.3

# update plot
plot ldata delchi

# save plot
cpd 520_spec_fit.ps/cps
plot ldata delchi
# to view the plot again use cpd /xw

# print model to screen
newpar 0

# get 1 sigma errors on parameters
error 1.0 2 4 7

# get unabsorbed flux by turning off nh
newpar 1 0
dummyrsp
# flux in 0.5-2.0 keV band
flux 0.5 2.0
# bolometric flux
flux 0.001 100

My spectrum and best fitting model looks like this