这是2008年Paolo Tozzi 来北京做的X射线数据处理讲座的笔记,一直拖到现在才开始整理。这里是第一部分,主要是些准备工作。高能数据处理要用HEASOFT,Chandra数据还要用到Ciao和CALDB,按照说明文件一步步安装并不困难,只是2G多文件下载起来要费些时间。附加脚本也经常用到,解压到安装路径下就可以了。安装好之后,先用 heainit 和ciao初始化运行环境,若没有错误提示就可以运行下面的命令了。
如果手头上没有数据。可以先到Chandra数据中心找个自己喜欢的目标,数据通常会在保密期过后自动公开,既避免了对同一目标的反复观测,也能最大限度地发挥已有数据的价值。 我这里选择了红移0.3的 MACSJ0257.6-2209 (即Abell 402, Obs ID 3267 )。输入到左上角的Target Name中,点Search,会得Chandra已有的观测数据列表。我们选个曝光时间短点的,数据文件小些,加到获取列表中(Add to Retrieval List)。点击retrieval之后,服务器会自动将文件打包放到匿名ftp的临时目录下。也有一个基于python的自动下载脚本。 官方详尽的下载介绍由此跳转。
压缩包里东西很多,图像处理只需要部分文件:
primary目录下,记录坏点的 bpix,确定指向的 asol
secondary目录下,观测事例 evt ,描述可用时间间隔的 flt,储存观测参数的 pbk 文件,
后面的数字代表处理级别,0为原始数据,1为基本处理数据,2为科学数据
完整命名规则由此去。
为防止误操作原始数据,简化文件名,先做符号链接
1 2 3 4 5 6 7 |
gunzip *.gz ln -s *evt1.fits evt1.fits ln -s *bpix1.fits bpix1.fits ln -s *asol1.fits asol1.fits ln -s *flt1.fits flt1.fits ln -s *pbk0.fits pbk0.fits ls -l |
启动ciao,用dmlist检查fits文件头中的基本信息。检查文件区块
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
$dmlist evt1.fits blocks -------------------------------------------------------------------------------- Dataset: evt1.fits -------------------------------------------------------------------------------- Block Name Type Dimensions -------------------------------------------------------------------------------- Block 1: PRIMARY Null Block 2: EVENTS Table 16 cols x 597013 rows Block 3: GTI3 Table 2 cols x 2 rows Block 4: GTI2 Table 2 cols x 3 rows Block 5: GTI1 Table 2 cols x 2 rows Block 6: GTI0 Table 2 cols x 2 rows Block 7: GTI6 Table 2 cols x 2 rows |
可以看到 fits 文件实际是按表格组织的,GTI是 可用时间间隔(Good Time Interval),记录起止时间每次两列。光子事件共有 597013 行,16列的具体含义可以用”dmlist evt1.fits cols”查看。
在ACIS的10片ccd中有两片是后照式BI(S1和S3),处理方式略有不同。详细介绍参考申请指南中的相关页面 下面以ACIS-I的四片FI式CCD为例简要介绍观测数据的定标处理过程。
1 |
dmcopy "evt1.fits[ccd_id=0,1,2,3]" evt1_acis_I.fits |
将CCD数据提取到新文件中,然后用acis_process_events作初步的处理。
1 2 3 4 |
$ dmlist evt1.fits header | grep DATAMODE | awk '{print $3}' VFAINT # 检查观测模式, 确定下面的 check_vf_pha 参数 $ dmlist evt1.fits header | grep Focal | grep FP_TEMP| awk '{print $3}' 153.44601440 # 检查焦面温度,153K需要做 CTI 修正, 对应下面的 apply_cti 参数 |
下面用 punlearn 更新命令参数,虽然也可以直接在命令后加参数,不过前者看起来比较清楚。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 |
punlearn acis_process_events pset acis_process_events infile=evt1_acis_I.fits pset acis_process_events outfile=evt1_updated.fits pset acis_process_events eventdef=")stdlev1" pset acis_process_events check_vf_pha=yes pset acis_process_events apply_cti=yes pset acis_process_events rand_pha=yes pset acis_process_events doevtgrade=yes pset acis_process_events apply_tgain=yes pset acis_process_events gradefile=CALDB pset acis_process_events gainfile=CALDB pset acis_process_events calculate_pi=yes pset acis_process_events badpixfile=bpix1.fits pset acis_process_events acaofffile=asol1.fits pset acis_process_events stop=none pset acis_process_events mode=h pset acis_process_events verbose=5 pset acis_process_events clobber = yes acis_process_events |
Chandra的正常光子事件是在可用时间间隔内 grade 为 [02346]之一,状态为0的那部分。我们用dmcopy来提取。
1 |
dmcopy "evt1_updated.fits[EVENTS][grade=0,2,3,4,6][status=0][@flt1.fits]" evt2.fits |
这就是从evt1到evt2的标准处理流程,如果没有什么特殊要求的话,直接用primary目录下的evt2文件也是一样的。在ds9中检查一下,应该和压缩包中的 full_img 图像是一样的。
为避免穿越辐射带,遭遇太阳风暴等情况带来的高能背景涨落,还要选取背景信号稳定的观测时段。下面统一处理四个CCD,谨慎一点的话可以对他们分别处理。
先在ds9中手动选出包括星系团在内的所有亮源,保存成souce.reg文件。然后扣除这些区域提取背景光谱,以免X射线源流量变化带来的干扰。接下来观察背景光子计数率变化曲线,时间间隔取100秒。
1 2 3 4 5 6 7 8 |
punlearn dmextract pset dmextract "evt2.fits[exclude sky=region(source.reg)][bin time=::100]" pset dmextract outfile=check.pha pset dmextract clobber=yes pset dmextract opt=ltc1 pset dmextract mode=h pset dmextract verbose=5 dmextract |
借助现成的 S-Lang 脚本,详细步骤可参考此页面。
1 2 3 4 5 6 7 8 9 10 11 12 13 |
slsh slsh> require ("lightcurves"); slsh> lc_sigma_clip("check.pha"); .... GTI limits calculated using a count-rate filter: (count_rate>1.2454 && count_rate<2.01135) The corresponding times are: ((time >= 121965879.253540) && (time < 121986779.253540)) ; 20.47 ksec, bin 1 .... slsh>quit |
程序默认剔除所有计数率偏离均值3σ以上的异常时段。
我们把脚本输出的正常时段保存成新的gti文件来更新源文件,注意要删除其中的等号,否则报错。
1 2 3 4 5 6 7 8 |
punlearn dmgti pset dmgti infile=evt2.fits pset dmgti outfile=evt2.gti pset dmgti userlimit = " ((time > 121965879.253540) && (time < 121986779.253540))" pset dmgti clobber=yes pset dmgti mode=h pset dmgti smooth=yes dmgti |
更新fits文件
1 |
dmcopy "evt2.fits[@evt2.gti]" evt2_clean.fits |
不放心的话可以再重新提取背景光谱检查一下。如果是对几片CCD分别处理,还要对其余三片重复上述过程,然后再把它们重新合并起来
1 2 3 4 5 6 7 |
punlearn dmmerge pset dmmerge infile=ccd0_clean.fits, ccd1_clean.fits, ccd2_clean.fits, ccd3_clean.fits pset dmmerge outfile=evt2_clean.fits pset dmmerge outBlock=events pset dmmerge mode=h pset dmmerge clobber=yes dmmerge |
如果有多次观测可以在分别处理后通过merge_all合并。定标就算完成了,接下来终于可以提取数据了。在命令dmcopy中,可以通过 energy 来选取波段,用bin来选取或合并像素
1 |
dmcopy "evt2_clean.fits[events][energy=500:5000][bin x=2000:5500:2,y=2000:5500:2]" img_tot.fits |
为了更好地展示图像,可将不同的能段赋予不同颜色,产生伪彩色图。这里我们选取0.3-1 keV 为红色,1-2keV 为绿色,2-5 keV 为蓝色
1 2 3 |
dmcopy "evt2_clean.fits[energy=300:1000][bin x=2000:5500:2,y=2000:5500:2]" img_cl_red.fits dmcopy "evt2_clean.fits[energy=1000:2000][bin x=2000:5500:2,y=2000:5500:2]" img_cl_green.fits dmcopy "evt2_clean.fits[energy=2000:5000][bin x=2000:5500:2,y=2000:5500:2]" img_cl_blue.fits |
在Chandra的高角分辨率下X射线光子分布十分弥散,我们可以用aconvolve命令做一下高斯平滑,方便查看。
1 2 3 4 5 6 7 8 9 10 |
punlearn aconvolve pset aconvolve clobber=yes pset aconvolve infile=img_cl_green.fits pset aconvolve outfile=img_cl_greensm.fits pset aconvolve kernelspec="lib:gaus(2,5,1,2.5,2.5)" pset aconvolve method=fft pset aconvolve verbose=5 pset aconvolve mode=h pset aconvolve normkernel = area aconvolve |
上面这步对三个波段都要做,我这里就不重复了。接下来用dmimg2jpg合并三色图像就可以得到伪彩色图了。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
punlearn dmimg2jpg pset dmimg2jpg infile=img_cl_redsm.fits pset dmimg2jpg greenfile=img_cl_greensm.fits pset dmimg2jpg bluefile=img_cl_bluesm.fits pset dmimg2jpg outfile=color.jpg pset dmimg2jpg maxred=0.7 maxblue=0.7 maxgreen=1.5 pset dmimg2jpg minred=0 minblue=0 mingreen=0 pset dmimg2jpg gridsize=60 pset dmimg2jpg fontsize=1 pset dmimg2jpg showaimpoint=no pset dmimg2jpg clobber=yes pset dmimg2jpg scalefunction=power pset dmimg2jpg scaleparam=0.6 pset dmimg2jpg mode=h dmimg2jpg |
当然,生成jpg格式主要是方便检查核对。图片确实不够漂亮,但包含了很多有用的信息。除了明亮的星系团,图中还有很多颜色各异的亮点,通常是AGN。蓝色的温度较高,由于低能段被介质吸收,因而缺少红色,不然便是白色。而红色的温度较低,高能光子较少。更具体的情况可以直接查看X射线光谱,将在下篇中介绍。