Yu-Ying Zhang

*


Computer Lab Course

Instructor: Dr. Yu-Ying Zhang


####################################################
##
## X-ray XMM-Newton EPIC data reduction manual (19/10/2012)
##
## Credit: Yu-Ying Zhang
## Argelander-Institut fuer Astronomie
## Rheinische Friedrich-Wilhelms-Universitaet Bonn
## Auf dem Huegel 71, D-53121 Bonn, Germany
## Phone: +49 228 733673
## Fax: +49 228 734022
##
## Relevant references:
## Each reference describes a different analyzing strategy concerning background treatment
## Zhang et al. 2004, A&A, 413, 49
## Zhang et al. 2005, A&A, 429, 85
## Zhang et al. 2006, A&A, 456, 55
## Zhang et al. 2007, A&A, 467, 437
## Zhang et al. 2009, ApJ, 699, 1178
##
####################################################


####################################################
##
## Why am I interested in XMM-Newton data?
##
## Large FOV; high spectral resolution; high effective area
##
####################################################


####################################################
##
## What can I do with XMM-Newton data?
##
## Intra-Cluster Medium
## hydrostatic mass: see eq. (2) in page 1182
##
####################################################


####################################################
##
## Where can I download XMM-Newton data?
##
## XMM-Newton Archive
##
####################################################


####################################################
##
## How can I set up the software?

bash

source setup.sh

## --------------- begin setup.sh
## Set up HEASOFT
## for light curve and spectral analysis

source /vol/software/software/astro/heasoft/heasoft-init.sh 6.11
source /vol/software/software/astro/heasarc_CALDB/software/tools/caldbinit.sh

## Set up CIAO
## for image analysis

#CIAO="/vol/software/software/astro/chandra/ciao/64/ciao-4.3"
CIAO="/vol/software/software/astro/chandra/ciao/64/ciao-4.4"
. ${CIAO}/bin/ciao.sh

## Set up SAS (XMM-Newton Science Analysis System)

export LHEASOFT=/vol/software/software/astro/heasoft/
export SAS_CCFPATH=/vol/software/software/astro/xmm/CCF/COMMON
#export SAS_CCF=./ccf.cif
. /vol/software/software/astro/xmm/xmmsas-init.sh
export SAS_VERBOSITY=4
PATH="$PATH:~/bin"

## --------------- end setup.sh
##
####################################################


####################################################
##
## Don't forget to setp up your path to my scripts.

##Download the tar:
##http://www.astro.uni-bonn.de/~yyzhang/WS/bin.tar.gz
##Otherwise: cd
tar zxvf /home/xray71/WS/bin.tar.gz
PATH="$PATH:~/bin"
##
#####################################################


####################################################
##
## How can I analyze the XMM-Newton data (e.g., galaxy cluster CL0024+17)?
##
## Results see Zhang et al. 2005, A&A, 429, 85
## I summarize the analysis into the following Steps:
##
####################################################


####################################################
##
## Steps:

(1) Untar the raw data.


## Create a working directory
cd
mkdir WS12

## Enter the working directory
cd WS12

## Unpack the raw data
tar xvf /home/xray71/WS/0050140201.tar
cd 0050140201/odf
tar zxvf 0050140201.tar.gz
rm -f 0050140201.tar.gz
tar xvf 0198_0050140201.TAR
rm -f 0198_0050140201.TAR
cd ..


(2) Build up the calibration file of my observation.


## Create a directory for primary scientific data
mkdir primary
cd primary

## Setup path to raw data
cd ../odf
odfdir=`pwd`

## Create data calibration
cd ../primary
export SAS_CCF=$SAS_CCFPATH
export SAS_ODF=$odfdir
cifbuild | tee ccf.log

## Point to newly created calibration
export SAS_CCF=./ccf.cif


(3) Create the calibrated data of my observation.


## Create a summary file
odfingest outdir=$SAS_ODF odfdir=$SAS_ODF
cd ../odf
file=`ls *SUM.SAS`
export SAS_ODF=${odfdir}/${file}

## Create event lists for all three EPIC instrument: MOS1+MOS2+pn
cd ../primary
emchain ## Create EPIC MOS event lists
ls P*M1*EV* P*M2*EV*
epchain odfaccess=all ## Create EPIC pn event list
ls P*PN*PIEV*
## Note the pn events has a bright column toward the readout node.
## For example,
ds9 flagP0112230301PNS003_eviw.fits.gz
## Therefore an EPIC pn out-of-time events list has to be created for corrections.
epchain odfaccess=all runbackground=N keepintermediate=raw withoutoftime=Y ## Create EPIC pn out-of-time events list
ls P*PN*OOEV*

## Pile-ups, pn mode, and out-of-time correction:
Pile-ups appear in the EPIC MOS and PN detectors when more than one photon
hits the same CCD pixel during the integration time of a frame (e.g. 68.702
ms for the EPIC PN in Full-Frame mode) and deposit their energies into the
detector. The result is that we are unable to tell if there were several
photons detected during this time or only one with the energy equivalent to
several photons. As a result of pile-ups the read out electronic registers
only one photon but with the energy of all involved photons. Therefore the
X-ray spectrum of a source looks harder than it actually is. The longer the
integration time and the brighter the source the larger is the chance to get
pile-ups in your observation. If you are just planning an observation, you
should use an observing mode that suits your source the best, e.g. it does
not make sense to observe a bright source with e.g. 20 cts/s with the PN in
Full-Frame mode. Use the Small-Window or even Timing Modes instead. The way
to check this before your observation is by estimating the pile-up fraction
by the mission Count Rate Simulator WebPIMMS at Goddard.


(4) Clean the calibrated data, and get rid of bad observing periods.


## Note evselect is a very important command in SAS for extracting information from the event lists, of which the parameters are generally defined here with some extra information here.
## In order to get a lightcurve with columns TIME and COUNTS, set maketimecolumn=true

## Create pn light curve in the hard band for pn to check particle background
evselect table=P0050140201PNS003PIEVLI0000.FIT withrateset=true \
maketimecolumn=true timebinsize=100 rateset=P0050140201PNS003_h_lc.fits \
expression="PATTERN.le.4 .and.FLAG.eq.0 .and. PI.ge.12000..and.PI.le.14000."

## Plot the light curve
fv P0050140201PNS003_h_lc.fits
## Click RATE-Plot, Supply X with TIME and y with COUNTS, click GO

## Export the light curve fits to an ASCII file
fdump infile=P0050140201PNS003_h_lc.fits outfile=P0050140201PNS003_h_lc.ascii columns=' ' rows='-' prhead- clobber=yes

## Make a histogram of the light curve
lc_qdp200 < P0050140201PNS003_h_lc.ascii
mv -f lc.qdp lc_pn_h.qdp

## Make a Gaussian fit to the historgram
qdp lc_pn_h.qdp
/xw
pl
re x 0 40
mo gaus
15
5
80
fit
fit
cpd /ps
pl
cpd /xw
rp
wmod lc_pn_h.mod
q

mv -f pgplot.ps lc_pn_h.ps

## Create light curve in the soft band for pn to check soft proton background
## In order to get a lightcurve with columns TIME and COUNTS, set maketimecolumn=true
evselect table=P0050140201PNS003PIEVLI0000.FIT withrateset=true \
maketimecolumn=true timebinsize=10 rateset=P0050140201PNS003_s_lc.fits \
expression="PATTERN.le.4 .and.FLAG.eq.0 .and.PI.ge.3000..and.PI.le.10000."

## Plot the light curve fits file
fv P0050140201PNS003_s_lc.fits
## Click RATE-Plot, Supply X with TIME and y with COUNTS, click GO

## Export the light curve fits to an ASCII file
fdump infile=P0050140201PNS003_s_lc.fits outfile=P0050140201PNS003_s_lc.ascii columns=' ' rows='-' prhead- clobber=yes

## Make a histogram of the light curve
lc_qdp200 < P0050140201PNS003_s_lc.ascii
mv -f lc.qdp lc_pn_s.qdp

## Make a Gaussian fit to the historgram
qdp lc_pn_s.qdp
/xw
pl
re x 0 40
mo gaus
15
5
80
fit
fit
cpd /ps
pl
cpd /xw
rp
wmod lc_pn_s.mod
q
mv -f pgplot.ps lc_pn_s.ps

## Create MOS1 light curve in the hard band for pn to check particle background
evselect table=P0050140201M1S001MIEVLI0000.FIT withrateset=true \
maketimecolumn=true timebinsize=100 rateset=P0050140201M1S001_h_lc.fits \
expression="PATTERN.le.12 .and.FLAG.eq.0 .and. PI.ge.10000..and.PI.le.12000."

## Plot the light curve
fv P0050140201M1S001_h_lc.fits
## Click RATE-Plot, Supply X with TIME and y with COUNTS, click GO

## Export the light curve fits to an ASCII file
fdump infile=P0050140201M1S001_h_lc.fits outfile=P0050140201M1S001_h_lc.ascii columns=' ' rows='-' prhead- clobber=yes

## Make a histogram of the light curve
lc_qdp50 < P0050140201M1S001_h_lc.ascii
mv -f lc.qdp lc_m1_h.qdp

## Make a Gaussian fit to the historgram
qdp lc_m1_h.qdp
/xw
pl
re x 0 20
mo gaus
10
5
100
fit
fit
cpd /ps
pl
cpd /xw
rp
wmod lc_m1_h.mod
q
mv -f pgplot.ps lc_m1_h.ps
## Create light curve in the soft band for M1 to check soft proton background
## In order to get a lightcurve with columns TIME and COUNTS, set maketimecolumn=true
evselect table=P0050140201M1S001MIEVLI0000.FIT withrateset=true \
maketimecolumn=true timebinsize=10 rateset=P0050140201M1S001_s_lc.fits \
expression="PATTERN.le.12 .and.FLAG.eq.0 .and.PI.ge.3000..and.PI.le.10000."

## Plot the light curve fits file
fv P0050140201M1S001_s_lc.fits
## Click RATE-Plot, Supply X with TIME and y with COUNTS, click GO

## Export the light curve fits to an ASCII file
fdump infile=P0050140201M1S001_s_lc.fits outfile=P0050140201M1S001_s_lc.ascii columns=' ' rows='-' prhead- clobber=yes

## Make a histogram of the light curve
lc_qdp100 < P0050140201M1S001_s_lc.ascii
mv -f lc.qdp lc_m1_s.qdp

## Make a Gaussian fit to the historgram
qdp lc_m1_s.qdp
/xw
pl
re x 0 10
pl
mo gaus
3.5
5
800
fit
fit
cpd /ps
pl
cpd /xw
rp
wmod lc_m1_s.mod
q
mv -f pgplot.ps lc_m1_s.ps

## Create MOS2 light curve in the hard band for pn to check particle background
evselect table=P0050140201M2S006MIEVLI0000.FIT withrateset=true \
maketimecolumn=true timebinsize=100 rateset=P0050140201M2S006_h_lc.fits \
expression="PATTERN.le.12 .and.FLAG.eq.0 .and. PI.ge.10000..and.PI.le.12000."

## Plot the light curve
fv P0050140201M2S006_h_lc.fits
## Click RATE-Plot, Supply X with TIME and y with COUNTS, click GO

## Export the light curve fits to an ASCII file
fdump infile=P0050140201M2S006_h_lc.fits outfile=P0050140201M2S006_h_lc.ascii columns=' ' rows='-' prhead- clobber=yes

## Make a histogram of the light curve
lc_qdp50 < P0050140201M2S006_h_lc.ascii
mv -f lc.qdp lc_m2_h.qdp

## Make a Gaussian fit to the historgram
qdp lc_m2_h.qdp
/xw
pl
re x 0 20
pl
mo gaus
10
5
100
fit
fit
cpd /ps
pl
cpd /xw
rp
wmod lc_m2_h.mod
q
mv -f pgplot.ps lc_m2_h.ps

## Create light curve in the soft band for M2 to check soft proton background
## In order to get a lightcurve with columns TIME and COUNTS, set maketimecolumn=true
evselect table=P0050140201M2S006MIEVLI0000.FIT withrateset=true \
maketimecolumn=true timebinsize=10 rateset=P0050140201M2S006_s_lc.fits \
expression="PATTERN.le.12 .and.FLAG.eq.0 .and.PI.ge.3000..and.PI.le.10000."

## Plot the light curve fits file
fv P0050140201M2S006_s_lc.fits
## Click RATE-Plot, Supply X with TIME and y with COUNTS, click GO

## Export the light curve fits to an ASCII file
fdump infile=P0050140201M2S006_s_lc.fits outfile=P0050140201M2S006_s_lc.ascii columns=' ' rows='-' prhead- clobber=yes

## Make a histogram of the light curve
lc_qdp100 < P0050140201M2S006_s_lc.ascii
mv -f lc.qdp lc_m2_s.qdp

## Make a Gaussian fit to the historgram
qdp lc_m2_s.qdp
/xw
pl
re x 0 10
pl
mo gaus
3.5
5
800
fit
fit
cpd /ps
pl
cpd /xw
rp
wmod lc_m2_s.mod
q
mv -f pgplot.ps lc_m2_s.ps

## Create a table of the time intervals with tabgtigen, of which the counts is below
## 3.3 sigma above the mean of the light curve historgram
tabgtigen table=P0050140201M1S001_h_lc.fits gtiset=P0050140201M1S001_h_gti.fits expression="COUNTS.lt.18.4476"
tabgtigen table=P0050140201M1S001_s_lc.fits gtiset=P0050140201M1S001_s_gti.fits expression="COUNTS.lt.11.1307"
tabgtigen table=P0050140201M2S006_h_lc.fits gtiset=P0050140201M2S006_h_gti.fits expression="COUNTS.lt.17.5311"
tabgtigen table=P0050140201M2S006_s_lc.fits gtiset=P0050140201M2S006_s_gti.fits expression="COUNTS.lt.11.5491"
tabgtigen table=P0050140201PNS003_h_lc.fits gtiset=P0050140201PNS003_h_gti.fits expression="COUNTS.lt.29.1474"
tabgtigen table=P0050140201PNS003_s_lc.fits gtiset=P0050140201PNS003_s_gti.fits expression="COUNTS.lt.23.4005"

## Create event lists in good time intervals
## M1
evselect table=P0050140201M1S001MIEVLI0000.FIT withfilteredset=yes \
destruct=yes keepfilteroutput=yes updateexposure=yes writedss=yes \
expression='GTI('P0050140201M1S001_h_gti.fits',TIME).and.GTI('P0050140201M1S001_s_gti.fits',TIME).and.FLAG.eq.0.' \
filteredset=flagP0050140201M1S001_clean.fits
## M2
evselect table=P0050140201M2S006MIEVLI0000.FIT withfilteredset=yes \
destruct=yes keepfilteroutput=yes updateexposure=yes writedss=yes \
expression='GTI('P0050140201M2S006_h_gti.fits',TIME).and.GTI('P0050140201M2S006_s_gti.fits',TIME).and.FLAG.eq.0.' \
filteredset=flagP0050140201M2S006_clean.fits
## pn
evselect table=P0050140201PNS003PIEVLI0000.FIT withfilteredset=yes \
destruct=yes keepfilteroutput=yes updateexposure=yes writedss=yes \
expression='GTI('P0050140201PNS003_h_gti.fits',TIME).and.GTI('P0050140201PNS003_s_gti.fits',TIME).and.FLAG.eq.0.' \
filteredset=flagP0050140201PNS003_clean.fits
## pn out-of-time events
evselect table=P0050140201PNS003OOEVLI0000.FIT withfilteredset=yes \
destruct=yes keepfilteroutput=yes updateexposure=yes writedss=yes \
expression='GTI('P0050140201PNS003_h_gti.fits',TIME).and.GTI('P0050140201PNS003_s_gti.fits',TIME).and.FLAG.eq.0.' \
filteredset=flagP0050140201PNS003OO_clean.fits


(5) Create a weighting column in the calibrated data, which accounts for vignetting effect.


Want to know more about the vignetting effect, click here.

## Weight EPIC event lists with inverse effective area with evigweight.
evigweight ineventset=flagP0050140201M1S001_clean.fits newoutput=yes outeventset=flagP0050140201M1S001_eviw.fits
evigweight ineventset=flagP0050140201M2S006_clean.fits newoutput=yes outeventset=flagP0050140201M2S006_eviw.fits
evigweight ineventset=flagP0050140201PNS003_clean.fits newoutput=yes outeventset=flagP0050140201PNS003_eviw.fits
evigweight ineventset=flagP0050140201PNS003OO_clean.fits newoutput=yes outeventset=flagP0050140201PNS003OO_eviw.fits


(6) Extract the images.


## For example, wth a binsize of 80 pixels.
evselect table=flagP0050140201M1S001_eviw.fits xcolumn=X ycolumn=Y imagebinning=binSize ximagebinsize=80 yimagebinsize=80 \
withimageset=true imageset=flagP0050140201M1S001_0520.fits expression="PATTERN.le.12..and.FLAG.eq.0..and.PI.ge.500..and.PI.le.2000."

evselect table=flagP0050140201M2S006_eviw.fits xcolumn=X ycolumn=Y imagebinning=binSize ximagebinsize=80 yimagebinsize=80 \
withimageset=true imageset=flagP0050140201M2S006_0520.fits expression="PATTERN.le.12..and.FLAG.eq.0..and.PI.ge.500..and.PI.le.2000."

evselect table=flagP0050140201PNS003_eviw.fits xcolumn=X ycolumn=Y imagebinning=binSize ximagebinsize=80 yimagebinsize=80 \
withimageset=true imageset=flagP0050140201PNS003_0520.fits expression="PATTERN.le.4..and.FLAG.eq.0..and.PI.ge.500..and.PI.le.2000."

(7) Extract the spectra.


## Check the data and decide the region to extract the spectrum
## Note 1'=1200 pixels in physical coordinates
ds9 flagP0050140201M1S001_eviw.fits &

## Extract global spectrum
## M1
evselect table=flagP0050140201M1S001_eviw.fits withspectrumset=true withspecranges=true \
energycolumn=PI specchannelmin=0 specchannelmax=11999 spectralbinsize=15 \
expression="circle(27764.502,26813.536,1200,X,Y) .and. PATTERN.le.12 .and. FLAG.eq.0" \
spectrumset=detsrc1200_1200m1.pha withzcolumn=Y withzerrorcolumn=N
## M2
evselect table=flagP0050140201M2S006_eviw.fits withspectrumset=true withspecranges=true \
energycolumn=PI specchannelmin=0 specchannelmax=11999 spectralbinsize=15 \
expression="circle(27764.502,26813.536,1200,X,Y) .and. PATTERN.le.12 .and. FLAG.eq.0" \
spectrumset=detsrc1200_1200m2.pha withzcolumn=Y withzerrorcolumn=N
## pn
evselect table=flagP0050140201PNS003_eviw.fits withspectrumset=true withspecranges=true \
energycolumn=PI specchannelmin=0 specchannelmax=20479 spectralbinsize=5 \
expression="circle(27764.502,26813.536,1200,X,Y) .and. PATTERN.le.4 .and. FLAG.eq.0" \
spectrumset=detsrc1200_1200pn.pha withzcolumn=Y withzerrorcolumn=N
## pn out-of-time events
evselect table=flagP0050140201PNS003OO_eviw.fits withspectrumset=true withspecranges=true \
energycolumn=PI specchannelmin=0 specchannelmax=20479 spectralbinsize=5 \
expression="circle(27764.502,26813.536,1200,X,Y) .and. PATTERN.le.4 .and. FLAG.eq.0" \
spectrumset=detsrc1200_1200pnoo.pha withzcolumn=Y withzerrorcolumn=N

## The generation of the effective area ARF (Ancillary Response File) by the task arfgen.
## The generation of the detector response file RMF (Redistribution Matrix File) by the task rmfgen.
## M1 rmf
rmfgen spectrumset=detsrc1200_1200m1.pha rmfset=m1on.rmf
## M1 arf
arfgen spectrumset=detsrc1200_1200m1.pha arfset=m1on.arf \
withrmfset=yes rmfset=m1on.rmf withsourcepos=yes sourcecoords="tel" \
sourcex=0 sourcey=0 withdetbounds=yes withbadpixcorr=yes \
badpixlocation=flagP0050140201M1S001_eviw.fits \
extendedsource=yes detmaptype=flat detxbins=10 detybins=10
## M2 rmf
rmfgen spectrumset=detsrc1200_1200m2.pha rmfset=m2on.rmf
## M2 arf
arfgen spectrumset=detsrc1200_1200m2.pha arfset=m2on.arf \
withrmfset=yes rmfset=m2on.rmf withsourcepos=yes sourcecoords="tel" \
sourcex=0 sourcey=0 withdetbounds=yes withbadpixcorr=yes \
badpixlocation=flagP0050140201M2S006_eviw.fits \
extendedsource=yes detmaptype=flat detxbins=10 detybins=10
## pn rmf
rmfgen spectrumset=detsrc1200_1200pn.pha rmfset=pnon.rmf
## pn arf
arfgen spectrumset=detsrc1200_1200pn.pha arfset=pnon.arf \
withrmfset=yes rmfset=pnon.rmf withsourcepos=yes sourcecoords="tel" \
sourcex=0 sourcey=0 withdetbounds=yes withbadpixcorr=yes \
badpixlocation=flagP0050140201PNS003_eviw.fits \
extendedsource=yes detmaptype=flat detxbins=10 detybins=10

## Prepare local background
## As shown in Zhang et al. 2005, A&A, 429, 85, the r200 for CL0024+17 is at 1.05Mpc which is 3.28' at redshift z=0.395.
## The region beyond r200 has little cluster X-ray emission, which can be used as local background (i.e., 3.5' - 3.64')
## M1
evselect table=flagP0050140201M1S001_eviw.fits withspectrumset=true withspecranges=true \
energycolumn=PI specchannelmin=0 specchannelmax=11999 spectralbinsize=15 \
expression="ring(27764.502,26813.536,4200,4368,X,Y) .and. PATTERN.le.12 .and. FLAG.eq.0" \
spectrumset=detbkg1200_1200m1.pha withzcolumn=Y withzerrorcolumn=N
## M2
evselect table=flagP0050140201M2S006_eviw.fits withspectrumset=true withspecranges=true \
energycolumn=PI specchannelmin=0 specchannelmax=11999 spectralbinsize=15 \
expression="ring(27764.502,26813.536,4200,4368,X,Y) .and. PATTERN.le.12 .and. FLAG.eq.0" \
spectrumset=detbkg1200_1200m2.pha withzcolumn=Y withzerrorcolumn=N
## pn
evselect table=flagP0050140201PNS003_eviw.fits withspectrumset=true withspecranges=true \
energycolumn=PI specchannelmin=0 specchannelmax=20479 spectralbinsize=5 \
expression="ring(27764.502,26813.536,4200,4368,X,Y) .and. PATTERN.le.4 .and. FLAG.eq.0" \
spectrumset=detbkg1200_1200pn.pha withzcolumn=Y withzerrorcolumn=N
## pn out-of-time events
evselect table=flagP0050140201PNS003OO_eviw.fits withspectrumset=true withspecranges=true \
energycolumn=PI specchannelmin=0 specchannelmax=20479 spectralbinsize=5 \
expression="ring(27764.502,26813.536,4200,4368,X,Y) .and. PATTERN.le.4 .and. FLAG.eq.0" \
spectrumset=detbkg1200_1200pnoo.pha withzcolumn=Y withzerrorcolumn=N

## The fraction for the out-of-time events correction for pn
## For source spectrum
mathpha expr=""detsrc1200_1200pn.pha"-0.023*"detsrc1200_1200pnoo.pha"" \
outfil="detsrc1200_1200pn_pure.pha" exposure="detsrc1200_1200pnoo.pha" \
ncomments=0 areascal="detsrc1200_1200pnoo.pha" units=R clobber=yes properr=yes ERRMETH="POISS-1"
## For local background
mathpha expr=""detbkg1200_1200pn.pha"-0.023*"detbkg1200_1200pnoo.pha"" \
outfil="detbkg1200_1200pn_pure.pha" exposure="detbkg1200_1200pnoo.pha" \
ncomments=0 areascal="detbkg1200_1200pnoo.pha" units=R clobber=yes properr=yes ERRMETH="POISS-1"

## Rebin the data with grppha
## The grouping is set such that each new grouping contains a minimum of certain counts (i.e., 100) in each bin.
## Channels that are defined as BAD are not included.
## Any spare channels at the end of the data are defined BAD by software.
## M1
grppha detsrc1200_1200m1.pha detsrc1200_1200m1_100g.pha
group min 100
exit
grppha detsrc1200_1200m1_100g.pha detsrc1200_1200m1_100g.pha \
comm="chkey respfile m1on.rmf & chkey ANCRFILE m1on.arf & exit" \
clobber=yes

## M2
grppha detsrc1200_1200m2.pha detsrc1200_1200m2_100g.pha
group min 100
exit
grppha detsrc1200_1200m2_100g.pha detsrc1200_1200m2_100g.pha \
comm="chkey respfile m2on.rmf & chkey ANCRFILE m2on.arf & exit" \
clobber=yes

## pn
grppha detsrc1200_1200pn_pure.pha detsrc1200_1200pn_pure_100g.pha
group min 100
exit
grppha detsrc1200_1200pn_pure_100g.pha detsrc1200_1200pn_pure_100g.pha \
comm="chkey respfile pnon.rmf & chkey ANCRFILE pnon.arf & exit" \
clobber=yes


## Prepare the hydrogen column density for the spectral fit
## I often take the weighted average nH from the LAB survey, which is 4.03E+20 /cm/cm for CL0024+17
nh
2000
6.6486614
17.159701


(8) Fit the spectra in XSPEC.


## Spectral fitting in XSPEC
xspec
cpd /xw
data 1:1 detsrc1200_1200m1_100g.pha
backgrnd 1 detbkg1200_1200m1.pha
data 2:2 detsrc1200_1200m2_100g.pha
backgrnd 2 detbkg1200_1200m2.pha
data 3:3 detsrc1200_1200pn_pure_100g.pha
backgrnd 3 detbkg1200_1200pn_pure.pha
setplot energy
ignore 1:**-0.7 7.8-**
ignore 2:**-0.7 7.8-**
ignore 3:**-0.7 7.8-**
ignore bad
plot ldata
statistic chi
query yes
model wabs*mekal
0.0403
5.0
1.e-6
0.3
0.395
1
0.1
0.0403
5.0
1.e-6
0.3
0.395
1
0.1
0.0403
5.0
1.e-6
0.3
0.395
1
0.1
untie 14 21
freeze 1 3 5 6 8 10 12 13 15 17 19 20
thaw 2 4 7 9 11 14 16 18 21
fit
newpar 9 = 2
newpar 16 = 2
newpar 11 = 4
newpar 18 = 4
fit
show all
pl
error 1. 2 4
cpd /cps
plot l chi
cpd /xw
gv pgplot.ps
log xspec.log
show all
error 1.0 2 4
log none
emacs xspec.log

:)