XSPEC光谱分析

在得到X波段图像之后,我们可以通过拟合光谱来测量天体温度、金属丰度等信息,就要用到XSPEC。作为X射线光谱分析的通用软件,已经被包含在NASA的高能物理套件HEASOFT之中。是X射线分析套件XANADU的主要组件。这个名字本来是元上都的英译,马可波罗曾在他的游记中极力描绘那里的奢华繁荣,后又经英国诗人柯勒律治题诗歌颂,而演变为东方仙境的代名词。

言归正传,下面我们采用上文得到的 MACS J0257.6-2209 定标文件来提取光谱。如果不想重复前面的步骤,只需提取压缩包primary目录下的evt2.fits 文件,执行下面的命令即可。

dmcopy "evt2.fits[events][energy=500:7000][bin x=2000:5500:2,y=2000:5500:2]" img_tot.fits

因为星系团是展源,为了消除背景AGN的影响,我们先用wavdetect做小波探测找出图中的所有点源。

punlearn wavdetect
pset wavdetect infile=img_tot.fits
pset wavdetect outfile=sources_tot.fits
pset wavdetect scellfile=cell_tot.fits
pset wavdetect imagefile=image_tot.fits
pset wavdetect defnbkgfile=background_tot.fits
pset wavdetect expfile=none
pset wavdetect scales="2 4 8 16"
pset wavdetect clobber=yes
pset wavdetect sigthresh = 1e-06
pset wavdetect regfile=sources_tot.reg
pset wavdetect mode=h
wavdetect

可以直接在事件文件中查看找到了多少点源

dmlist sources_tot.fits blocks

也可以直接在ds9中显示出来。

 ds9 img_tot.fits -region  sources_tot.reg& 

点源探测

检查是否有误判的区域,确认无误后,就可以剔除点源来提取光谱了。

dmcopy "evt2.fits[exclude sky=region(sources_tot.reg)]" field_nopt.fits
dmcopy "field_nopt.fits[events][energy=500:7000][bin x=2000:5500:2,y=2000:5500:2]" img_cluster.fits
ds9 img_cluster.fits&

在ds9中选出整个星系团。保存位置文件(cl.reg )。然后就可以提取光谱了。

dmcopy "field_nopt.fits[events][sky=region(cl.reg)]" evt2_cl.fits
dmextract "evt2_cl.fits[events][bin pi]" cl.pha clobber=yes
 

但是,我们选取的星系团区域中,并不全是来自星系团的光子。空间高能粒子,银河系X射线背景,无法分辨的遥远AGN都混杂其中。同可见光观测要做平场一样,我们还需要背景文件。获取背景文件有两种途径,一种是望远镜对空白天区的长时间曝光blank-sky background,但由于望远镜本身的状态变化,不同时期的背景文件并不通用,需要根据状态参数找到对应的修正文件。这个方法有些复杂,我们这里采用另外一种:利用CCD上的空白区域来提取背景。既然图像上的点源和星系团都已经挑出来了,剩下的就都是背景了。但选取时仍需要注意,尽量选取靠近星系团的区域,尽量使用目标所在的CCD,同时不要和图像边缘靠得太近。但如果星系团很大,占据了整个画面,就只有前一个方法可用了。

dmcopy "field_nopt.fits[events][sky=region(bck.reg)]" evt2_bck.fits
dmextract "evt2_bck.fits[events][bin pi]" bck_cl.pha
 

由于X射线望远镜对不同能段光子的有效探测面积和能量分辨率都不同,而且还随时间和位置变化。因此每次观测,我们都需要根据星系团所在位置计算响应文件。包括修正光子分布的Redistribution Matrix Function (RMF)和记录量子转换效率和有效面积的Ancillary Response File (ARF)。更具体的介绍可以参考MIT工作组的文档

dmcopy infile="evt2_cl.fits[energy=300:2000][bin det=8]" outfile=cl_sou.wmap clobber=yes

punlearn mkwarf
mkwarf infile=cl_sou.wmap outfile=arf_cl.fits weightfile=cl_sou.wfef spectrumfile="none" egridspec=0.3:11.0:0.01 pbkfile=pbk0.fits  clobber=yes

punlearn mkacisrmf
mkacisrmf infile=CALDB outfile=rmf_cl.fits wmap=cl_sou.wmap energy=0.3:11.0:0.01 channel=1:1024:1 chantype=PI gain=CALDB mode=h

生成的响应文件可以通过rmfimg转换为图像来检查。

pset rmfimg infile = rmf_cl.fits
pset rmfimg outfile = imgrmf_cl.fits
pset rmfimg arf = arf_cl.fits
pset rmfimg arfout = imgarf_cl.fits
pset rmfimg product = no
pset rmfimg verbose = 3
pset rmfimg clobber = yes
pset rmfimg mode = h
rmfimg

ds9 imgrmf_cl.fits&
ds9 imgarf_cl.fits&

响应文件

左侧为ARF右侧为RMF,横轴均为能道编号,纵轴为光子原始能量。对于ARF,由于硬X射线掠射角增加,望远镜有效面积减小,探测器量子效率也同时降低,因此呈现自下而上的递减趋势。而RMF在理想状态下应为狭窄的对角直线,每个光子都进相应的能道,但实现起来非常困难,只能根据所在能道反推原始能量分布。
现在我们已经准备好了全部文件,将所有这些文件同光谱关联到一起。

pset grppha clobber=yes
grppha infile=cl.pha outfile=cl_grp.pha chatter=0 comm=' group min 1 & chkey RESPFILE rmf_cl.fits & chkey ANCRFILE arf_cl.fits & chkey BACKFILE bck_cl.pha & exit' 

在拟合光谱之前,还要考虑银河系晕中的中性氢消光,这与天体位置有关。

$ dmlist evt2_cl.fits header | grep NOM
0055 RA_NOM                      44.4640975717           Real8        Nominal RA
0056 DEC_NOM                    -22.1583582964           Real8        Nominal Dec
0057 ROLL_NOM                   355.2881562242           Real8        Nominal Roll
$ nh 2000   44.4640975717 -22.1583582964
  >> Leiden/Argentine/Bonn (LAB) Survey of Galactic HI
  LII , BII 210.101242  -60.925470
  Requested position at X and Y pixel   346.67     21.64
  Search nH in  4  X  4 box
  Each pixel is  0.675 deg 0.675 deg
  nH calculated using all points within
   1.0000 deg from input position
     RA       DEC      Dist        nH
   43.8807 -22.2531   0.5485    2.12E+20
   43.8991 -21.7153   0.6856    2.41E+20
   44.7581 -22.9792   0.8649    1.82E+20
   44.7833 -22.4343   0.4044    1.97E+20
   44.7998 -21.8885   0.4117    1.97E+20
  LAB >> Average nH (cm**-2)  2.06E+20
  LAB >> Weighted average nH (cm**-2)  2.06E+20

我们采用最后一个加权平均值。终于可以运行XSPEC了。

cpd /xw
#设置输出到屏幕,如果保存为文件则是 cpd file.ps /cps
data cl_grp.pha
# 载入光谱
show rate
# 显示计数率,曝光时间等信息
query yes
setplot energy
# 設定显示模式为能谱
abund angr
# 設定太阳金属丰度参考值,这里取Anders和Grevesse于1989发表的数据。详细列表
notice **-**
ignore bad
ignore 0.0-0.5
ignore 8.0-**
# 显示光谱 0.5-8.0 keV 的部分
setplot rebin 3 20
# 显示时合并部分通道,方便查看。后面两个值分别是最小信噪比和最大bin size。
plot ldata
# 对数显示,就可以看到文章开头的光谱了。

statistic cstat
# 采用 C-statistic 统计
model tbabs*mekal
# 定义模型,下面按提示输入参数。注意这里的中性氢密度nH 单位是10^22 atoms/cm^–2。详见Xspec用户手册。温度(2)、金属丰度(4) 以及归一化参数(7)是我们要拟合的,随便给个初值。红移在BAX星系团数据库中可以查到: 0.322(Simbad 和 NED中还没有)。

Input parameter value, delta, min, bot, top, and max values for ...
              1      0.001          0          0     100000      1e+06
1:TBabs:nH>0.0206
              1       0.01     0.0808     0.0808       79.9       79.9
2:mekal:kT>5
              1      -0.01      1e-06      1e-05      1e+19      1e+20
3:mekal:nH>0.001
              1      -0.01          0          0       1000       1000
4:mekal:Abundanc>0.3
              0      -0.01          0          0         10         10
5:mekal:redshift>0.322
              1
6:mekal:switch>1
              1       0.01          0          0      1e+24      1e+24
7:mekal:norm>0.0001

freeze 1 3 5 6
thaw 2 4 7
#設定自由参数为2 4 7
renorm
fit 20 1e-2
#默认采用 Levenberg-Marquardt 算法拟合,需指定迭代次数,精度可选。

Model TBabs <1> * mekal <2> Source No.: 1   Active/On
Model Model Component  Parameter  Unit     Value
 par  comp
   1    1   TBabs      nH         10^22    2.06000E-02  frozen
   2    2   mekal      kT         keV      10.7532      +/-  0.920716     
   3    2   mekal      nH         cm-3     1.00000E-03  frozen
   4    2   mekal      Abundanc            0.438082     +/-  0.129011     
   5    2   mekal      redshift            0.322000     frozen
   6    2   mekal      switch              1            frozen
   7    2   mekal      norm                2.24587E-03  +/-  9.69028E-05  

steppar 2 0.1 10 20 4 0.1 0.6 20
# 如果为单参则直接给出一维分布
二参概率分布
plot contour
# 显示结果如右图所示。可以看出温度在10keV左右误差较小,金属丰度不确定性较大。下面计算这两个参数的 $$1\sigma$$误差,虽然前面的拟合结果中也有,但并不可靠。如果没有steppar一步,算误差时还可能发现新的最优值,所以还是再算一下比较放心。

error 1.0 2 4
 Parameter   Confidence Range (1)
     2      9.84033      11.6813    (-0.912888,0.928126)
     4     0.312238     0.573214    (-0.125844,0.135133)

setplot device spec.eps /cps
# 输出到指定文件
setplot rebin 3 20
#setplot 只改变显示效果,不影响实际数据,因此每次绘图都要重设
plot ldata res
# 增加残差部分
quit

这样我们就完成了X射线光谱拟合,得到了星系团MACS J0257.6-2209 的温度 T= 10.75 +/- 0.92keV 以及金属丰度 Z=0.44 +/- 0.13 。

拟合结果

参考资料:《Spectral fitting: XSPEC》2009

订阅评论
提醒
5 评论
最旧
最新 最多投票
内联反馈
查看所有评论
hbcai
2010 年 8 月 13 日 11:22

呵呵,我也做过chandra的数据处理。这个处理过程写的比较详细。方便后学。盼继续!

赞!!

hbcai
2010 年 8 月 13 日 17:02

还差分析Timing的。呵呵。

GH
2011 年 3 月 19 日 02:52

写的真好!当年你教了我们那么多东西,我都没好好学,现在要用了,上网一搜,还得是你教得好!谢谢拉!