Astro

Heasoft使用与CIAO软件包使用（一）

本文主要演示了Heasoft和CIAO部分功能的使用，涉及数据获取、彩图合成、光变曲线提取与光谱提取（X-ray方面）

Step1.数据获取

数据获取网址：https://heasarc.gsfc.nasa.gov/cgi-bin/W3Browse/w3browse.pl

记住这里的obsid编号，稍后需要使用，我们取4875为例子

数据下载

进入Linux环境，打开终端，我们先创建一个工作目录并进入

1
2
mkdir xray
cd xray/

先前下载好了heasoft软件包和CIAO软件包，CIAO软件包是专门为Chandra X-ray望远镜设计的一套软件包，建立好了进入环境的快捷指令，直接运行即可进入CIAO环境

1
ciao

我们安装Mrk1210的光谱数据，在先前我们获取了Mrk1210的编号4875

1
2
download_chandra_obsid 4875 #安装指令，下载到当前工作目录
chandra_repro 4875 outdir='' #解压缩

包含蟹状星云数据在内，工作目录应该有红框内的文件夹

Step2.数据成像

这一模块列举蟹状星云和Mrk1210的成像，截取了特定能量范围的光谱数据作为RGB三原色的FITS文件进行叠加，再通过色彩调整得到一幅伪彩色图。

Crab

1
2
cd 198/
cd primary/

我们可以看到有一个名称是acisf00198N005_evt2.fits.gz的压缩文件，无需解压，我们在终端命令行操作即可
提取不同能段的图像，由于X-ray能段并没有可见光那样的颜色，我们选取特定范围的X-ray作为不同的颜色，这样我们可以合成一张伪彩色的图片2

  • soft:0.2-1.5keV 红色
  • medium:1.5-2.5keV 绿色
  • hard: 2.5-8keV 蓝色
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    punlearn dmcopy
    # 提取三种波段的fits
    dmcopy infile='acisf00198N005_evt2.fits[energy=200:1500]' outfile='red_198.fits'

    dmcopy infile='acisf00198N005_evt2.fits[energy=1500:2500]' outfile='green_198.fits'

    dmcopy infile='acisf00198N005_evt2.fits[energy=2500:8000]' outfile='blue_198.fits'

    # RGB合并
    heasoft #进入heasoft环境,ds9是heasoft里的一个图像显示的软件
    ds9 -rgb -red red_198.fits -green green_198.fits -blue blue_198.fits
    效果如图：

调色

我们在左侧工具框中选择Lock，选中色阶图例和Smooth，然后按住鼠标右键在图上拖动，使颜色鲜艳，还可以在‘色阶’功能选择合适的颜色，或者在分析功能中选择Smooth Parameters进行调整

增加图例坐标轴

在分析工具中，我们选择坐标格线参数，Type一栏中，选中分析和两个Exterior Axes,在Grid中取消显示网格线

添加其他望远镜图样

回到主菜单,选择影像视窗的合重，这样可以显示多个图像，选择空白图像，我们在分析一栏的Archive选择一个望远镜，输入目标名称，选择其中一个图像载入即可

最终效果

Mrk1210

处理流程类似蟹状星云，仅展示最终效果图

1
2
3
4
5
6
7
8
9
10
11
12
13
cd
cd ./xray/4875/primary/
punlearn dmcopy
# 提取三种波段的fits
dmcopy infile='acisf04875N004_evt2.fits[energy=200:1500]' outfile='red_4875.fits'

dmcopy infile='acisf04875N004_evt2.fits[energy=1500:2500]' outfile='green_4875.fits'

dmcopy infile='acisf04875N004_evt2.fits[energy=2500:8000]' outfile='blue_4875.fits'

# RGB合并
heasoft #进入heasoft环境,ds9是heasoft里的一个图像显示的软件
ds9 -rgb -red red_4875.fits -green green_4875.fits -blue blue_4875.fits

Step3.光变曲线提取

能谱和光变曲线都能够从acisf04875_repro_evt2.fits提取出来。CIAO提供了一个很方便的教程Click_Here

选取目标源和背景区域范围

1
2
3
cd 
cd ./xray/4875/repro/
ds9 acisf04875_repro_evt2.fits

打开ds9后选择编辑的区域，选中源的位置，适当调整大小保证目标在圆圈内，然后选择区域功能的Save Selection，保存为src.reg，选择格式为CIAO。

在源周围再选定一个框，但不可覆盖源,重复上述步骤，保存为bkg.reg

提取lc

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
# 检查相机id，确认我们选取的源和背景的范围是在那个CCD芯片上（芯片的序号）
punlearn dmstat
dmstat "acisf04875_evt2.fits[sky=region(src.reg)][cols ccd_id]"
### 检查id为7

dmlist acisf04875_evt2.fits #查看时间范围，选取

#提取减掉背景信号的源的光变曲线
punlearn dmextract
pset dmextract infile="acisf04875_evt2.fits[ccd_id=7,sky=region(src.reg)][bin time=194800731.1376200020:194801731.1376200020:20]" #起始时间:最终时间:分bin时间
pset dmextract outfile="src_sub_lc.fits" #输出文件名称
pset dmextract bkg="acisf04875_evt2.fits[ccd_id=7,sky=region(bkg.reg)]" #引入背景
pset dmextract opt="ltc1" #光变曲线文件格式
dmextract #执行

dmlist src_sub_lc.fits cols #检查生成文件

#画图，使用heasoft的lcurve
heasoft #进入环境
lcurve 1 src_sub_lc.fits

# 运行信息，rebinning和interval的乘积就是选定的时间区间
$ Name of the window file ('-' for default window)    ##默认回车可看到一些时间等信息
$ Newbin Time or negative rebinning  20   ##输入重新并合的时间bin大小
$ Number of Newbins/Interval  50   ##输入时间重新bin后可用的通道数
$ Name of output file     ##输入输出文件名或默认回车
$ Do you want to plot your results?  ##输入yes或回车
$ Enter PGPLOT device    ##默认回车或输入 /xw

#修改文件头
punlearn dmhedit
dmhedit src_sub_lc.fits filelist='' operation=add key=object value=JHCao_mrk1210

#保存
hardcopy

最终成图

Step4.能谱拟合

工作目录依然是repro，因此不需要其他操作

方法一：一步提取点光源

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
CIAO #进入环境

#提取点源光谱
#specextract 可以从事件文件的头文件中确定辅助文件（宽高比解决方案、坏像素、掩码），这意味着只需提供源文件、输出目录和可能的背景文件即可： 
punlearn specextract
pset specextract infile="acisf04875_repro_evt2.fits[sky=region(src.reg)]" #src.reg和bkg.reg都是刚才使用ds9提取的点源和背景文件
pset specextract bkgfile="acisf04875_repro_evt2.fits[sky=region(bkg.reg)]"
pset specextract outroot=mrk1210
pset specextract correctpsf=yes
pset specextract weight=no
# weight 参数设置为 “no ”是为了告诉脚本做一个点样的 ARF，而 correctpsf 参数设置为 “yes ”是为了运行 arfcorr 来修正 ARF。这些参数必须更改，因为默认行为是 weight=yes 和 correctpsf=no。
specextract #运行

dmstat "acisf04875_repro_evt2.fits[sky=region(src.reg)][bin sky=2]"
pget dmstat out_max_loc #确定中心点

# specextract步骤生成的文件应该有：
$ mrk1210.pi            ：ungrouped source spectrum
$ mrk1210.arf           ：source ARF, no PSF corrections
$ mrk1210.corr.arf      ：source ARF, with PSF corrections applied
$ mrk1210.rmf           ：source RMF
$ mrk1210_grp.pi        ：grouped source spectrum

$ mrk1210_bkg.pi        ：background spectrum
$ mrk1210_bkg.arf       ：background ARF
$ mrk1210_bkg.rmf       ：background RMF

mrk1210_grp.pi就是我们最终需要的光谱文件

方法二：step by step

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
#提取源光谱
#我们提取标准脉冲不变（PI）通道的光谱。这样就绘制出了计数数与 PI 信道的直方图： 
pset dmextract infile="acisf04875_repro_evt2.fits[sky=region(src.reg)][bin pi]"
pset dmextract outfile=step_by_step.pi  #生成文件
pset dmextract op=pha1 #格式
dmextract #运行

#设置 ardlib.par (acis_set_ardlib)
#坏像素和坏列信息存储在 bpix1.fits 文件中。该文件被多个工具使用，可通过 ardlib.par 参数文件访问。在处理单个数据集时，chandra_repro 默认会用正确的坏像素文件自动设置 ardlib.par 文件。但是，如果其他观测数据已经过重新校准，则 ardlib.par 文件可能无法指向正确的文件。
punlearn ardlib
acis_set_ardlib acisf04875_repro_bpix1.fits absolute=no


#定位中心点（dmstat 和 dmcoords）
#由于各芯片的校准情况不尽相同，我们需要找到源的中心点（以芯片坐标为单位）。创建 ARFs 和选择使用 mkrmf 计算 RMFs 的 FEF（FITS 嵌入函数）时都需要这一信息。
dmstat "acisf04875_repro_evt2.fits[sky=region(src.reg)][bin sky=1]" centroid=yes
# 输出得到中心点坐标x=4064.3547562，y=4150.6423499 重要！！！

punlearn dmcoords #然后将其转换为 CHIP 坐标
dmcoords acisf13858_repro_evt2.fits op=sky x=4064.3547562 y=4150.6423499 #使用刚才得到的坐标
pget dmcoords chip_id chipx chipy
#输出：7 272.6542441495169 496.3159845780768 近似为272和496


#计算 RMFs
#望远镜使用的观测数据是在 -120 摄氏度焦平面温度下拍摄的，光源位于 ACIS-7 上，是一个背照式芯片。因此，可以使用 mkacisrmf 创建 RMF 文件。
pset mkacisrmf infile=CALDB
pset mkacisrmf wmap=none
pset mkacisrmf obsfile=acisf04875_repro_evt2.fits
pset mkacisrmf ccd_id=7 chipx=272 chipy=496
pset mkacisrmf chantype=PI channel=1:1024:1 energy=0.3:11.0:0.01
pset mkacisrmf outfile=step_by_step.rmf 
mkacisrmf #运行


#使用 mkrmf（acis_fef_lookup、mkrmf）首先需要 acis_fef_lookup 来确定正确的输入校准文件，即 FITS 嵌入函数或 FEF。
pset acis_fef_lookup infile=acisf04875_repro_evt2.fits 
pset acis_fef_lookup chipid=7 chipx=272 chipy=496
acis_fef_lookup
pset mkrmf infile=")acis_fef_lookup.outfile" 
pset mkrmf axis1=energy=0.3:11.0:0.01
pset mkrmf axis2=pi=1:1024:1
pset mkrmf outfile=step_by_step_fef.rmf
mkrmf # 运行


#计算ARFs
#计算高宽比直方图（asphist）现在我们需要创建高宽比直方图，它是观测过程中高宽比运动的分档表示。
#高宽比直方图使用事件文件中的 “良好时间间隔”（GTI）信息。由于有多个 GTI，每个 CCD_ID 都有一组，因此我们需要用正确的 CCD_ID 过滤器过滤事件文件，以选择正确的 GTI。
pset asphist infile=pcadf04875_001N001_asol1.fits
pset asphist evtfile="acisf04875_repro_evt2.fits[ccd_id=7]" 
pset asphist outfile=asp7.hist
asphist #运行


#制作ARFs
punlearn mkarf
pset mkarf asphistfile=asp7.hist
pset mkarf sourcepixelx=4064.3547562 sourcepixely=4150.6423499 #中心点坐标
pset mkarf engrid=0.3:11.0:0.01
pset mkarf obsfile=acisf04875_repro_evt2.fits
pset mkarf detsubsys="ACIS-7"
pset mkarf maskfile=acisf04875_001N004_msk1.fits
pset mkarf outfile=step_by_step_mkarf.arf
mkarf #运行


#修正源 ARF (arfcorr)
#arfcorr 计算源的能量相关点源孔径校正（PSF 分数），并将校正结果应用于 ARF 文件。运行 mkarf 时使用的 (x,y) 位置会输入到 arfcorr，源光谱使用的提取区域也是如此。需要输入图像来设置将要生成的 PSF 图像的大小和比例。
pset arfcorr infile="acisf04875_repro_evt2.fits[sky=region(src.reg)][bin sky]"
pset arfcorr region="region(src.reg)"
pset arfcorr x=4064.3547562 y=4150.6423499 #中心点坐标
pset arfcorr arf=step_by_step_mkarf.arf
pset arfcorr outfile=step_by_step_corrected.arf #输出
arfcorr #运行


#打包点源光谱
punlearn dmgroup
pset dmgroup infile=step_by_step.pi
pset dmgroup outfile=step_by_step_grp.pi
pset dmgroup grouptype=NUM_CTS 
pset dmgroup grouptypeval=15
pset dmgroup xcolumn=channel
pset dmgroup ycolumn=counts
dmgroup
#更新头文件
punlearn dmhedit
dmhedit step_by_step.pi filelist="" operation=add key=ANCRFILE value=step_by_step_corrected.arf 
dmhedit step_by_step.pi filelist="" operation=add key=RESPFILE value=step_by_step.rmf
dmhedit step_by_step_grp.pi filelist="" operation=add key=ANCRFILE value=step_by_step_corrected.arf
dmhedit step_by_step_grp.pi filelist="" operation=add key=RESPFILE value=step_by_step.rmf


#提取背景光谱
dmextract "acisf04875_repro_evt2.fits[sky=region(bkg.reg)][bin pi]" step_by_step_bkg.pi op=pha1


#计算背景ARFs，创建一个ARF的权重图
pset sky2tdet asphist=asp7.hist 
pset sky2tdet infile="acisf04875_repro_evt2.fits[sky=region(bkg.reg)][energy=300:2000][bin sky]" 
pset sky2tdet outfile="step_by_step_tdet_bkg.wmap[wmap]" 
sky2tdet #运行


#制作ARF
#输出文件中附加的 [wmap] 指定输出块名称应为 wmap。这样做是为了避免 mkwarf 发出一些非致命警告。
pset mkwarf egrid=0.3:11.0:0.01
pset mkwarf mskfile=acisf04875_001N004_msk1.fits
pset mkwarf infile=step_by_step_tdet_bkg.wmap
pset mkwarf outfile=step_by_step_bkg.warf
pset mkwarf feffile=CALDB
pset mkwarf weight=step_by_step_bkg.wghts
mkwarf


#制作背景RMFs
#使用 mkacisrmf 创建加权 RMF 也需要 WMAP。不过，由于 RMF 并不是在与 ARF 相同的空间尺度上进行校准的，因此用户可以通过 mkacisrmf 使用分档更粗的 wmap。wmap 可以在使用 dmextract 提取光谱的同时创建，也可以像这里显示的那样即时计算：
pset mkacisrmf infile=CALDB
pset mkacisrmf wmap="acisf04875_repro_evt2.fits[sky=region(bkg.reg)][energy=300:2000][bin tdet=8]"
pset mkacisrmf obsfile=acisf04875_repro_evt2.fits
pset mkacisrmf ccd_id=7 chipx=272 chipy=496   # these parameter value are required by not used
pset mkacisrmf chantype=PI channel=1:1024:1 energy=0.3:11.0:0.01
pset mkacisrmf outfile=step_by_step_bkg.rmf 
mkacisrmf #运行


#打包背景光谱
pset mkrmf infile=CALDB
pset mkrmf weight=step_by_step_bkg.wghts
pset mkrmf axis1=energy=0.3:11.0:0.01
pset mkrmf axis2=pi=1:1024:1
pset mkrmf outfile=step_by_step_fef_bkg.rmf
mkrmf

pset dmgroup infile=step_by_step_bkg.pi
pset dmgroup outfile=step_by_step_bkg_grp.pi
pset dmgroup grouptype=BIN
pset dmgroup binspec=1:1024:20
pset dmgroup xcolumn=channel
pset dmgroup ycolumn=counts
dmgroup #运行

绘制能谱

使用heasoft的Xspec对刚才输出的光谱进行能谱拟合

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
heasoft #进入环境
xspec #快捷指令进入xspec模块
cpd /xw #打开画板
data step_by_step_grp.pi #读取光谱文件
plot ldata #画出光谱
setplot energy   #设置很坐标为energy，默认为channel
plot ldata       ##重新画出光谱
setplot command rescale y 0.001 0.2    #后两位为纵坐标的上下范围
setplot command rescale x 0.4 9    #后两位为横坐标的上下范围
setplot command lab top caojh_mrk1210    #设置光谱图上方的显示字样，自定义名称
setplot add         #设置显示模型成分，默认为noadd不显示
model phabs(mekal+pow+zphabs(pow+gau))   ##利用模型拟合，回车后，要求输入红移的，输入红移值，第一个phabs的参数值，可以在高能网站查询

# 查询得到红移量为0.013496，nH为3.85E+20，此处nH单位是10^22因此换算为3.85E-02

fit    #利用上面定义的模型去拟合
plot ld res   #画出拟合后的光谱和残差，如设置了显示模型成分，会同时显示模型成分，这里只显示可以相加的模型的成分，相乘的成分不显示
hardcopy    #保存画板中显示的内容

红移量查询网站 https://ned.ipac.caltech.edu/
nH查询网站 https://heasarc.gsfc.nasa.gov/cgi-bin/Tools/w3nh/w3nh.pl

最终成图