费米LAT数据处理初体验·光谱分析篇
与其他波段相比，高能伽玛射线的光谱分析没有什么本质区别，无外乎得出一系列的谱参数，并且给出其统计意义上的置信度。不过高能光子观测有着自己的特点，在分析方法上，能谱处理自然也有与众不同之处，因此Fermi Science Tools为此准备了专门的工具，本文将相关的基本操作步骤讲解一遍。
说来这几个特点其实都不是什么优点。其中之一是由于光子数目过少（各能量bin中的光子数目基本遵从泊松分布），除非针对少量亮源，chi2检验一般不可用；其次是由于泊松噪声的存在，往往拟合结果不唯一；再有就是探测器性能所限，点扩散函数很宽，这样邻近的源非常容易彼此混淆，因此针对费米的谱处理，一般需要借助多变量分析才能完成。另外同样由于光子数目缺乏的原因，一般建议是不要再合并能段，否则结果起伏很大。限于篇幅，相关的统计理论本文不作详细介绍，有兴趣者请参阅相关文献和手册（如描述光子计数分析的Cash W., 1979, ApJ, 228, 939）。
进行谱分析时，首先需要确认光子的来源，这就要涉及探测器的响应了。响应函数的形式与LAT有效面积、点扩散函数以及一般可以忽略不计的能量色散有关，这后三者又都会涉及能量、辐射源对探测器轴向的倾角以及时间，所以分析过程中一个重要内容就是计算待分析天区中的响应。有效面积的形式已经在前文提过了，这里给出的是垂直入射时点扩散函数的形式：
在介绍基本步骤之前，照例先说一下需要用到的几个工具：
gtltcube：计算光子的livetime与辐射源的位置之间的关系。前面说过，LAT的响应函数与辐射源相对探测器轴向的夹角（所谓off-axis angle）有关。这样在巡天工作模式下，LAT记录下来自给定强度的辐射源的光子数目显然就取决于观测过程中辐射源处在LAT视场中的时间以及off-axis angle的变化情况。这一工具需要输入gtmktime筛选后的事件数据外加航天器数据。当然出于计算时间和效率的考虑，要对时间和空间合并，因此合并的步长是该工具的关键可调参数。
gtexpmap：生成似然分析所需要的exposure map，不过仅仅针对未合并能段的分析。这里的exposure map指的是响应函数针对目标天区的积分，需要提供事件数据、航天器数据以及gtltcube的结果。BTW，gtltcube和gtexpmap两步非常耗时，现在已经不太记得谁花费的时间更长了，反正完成这两个步骤需要的时间比一顿饭的工夫长得多……
gtlike：进行似然分析的工具。LAT的数据处理方法与本人之前对谱分析的印象不同，不是利用观测数据拟合结果，而是给定模型，给出获取观测数据的概率，这个概率被定义为似然参数L。gtlike的输入量包括gtbin生成的光子计数图、探测器数据、exposure map和gtltcube的结果，外加所有辐射源的光谱模型。
下面继续拿Terzan 5试刀。首先如前文所示，用gtselect、gtmktime筛选好数据，并用gtbin成图（注意输出文件参数与前文不同，要选成CMAP）。这些步骤在此不重复，接下来计算光子的livetime和exposure map：
[bzhang@Melipal terzan5]$ gtltcube
Event data file[ter5_filtered_5deg_gti.fits]
Spacecraft data file[L100618062133E0D2F37E56_SC00.fits]
Output file[expcube_me.fits]
Step size in cos(theta) (0.:1.) [0.025]
Pixel size (degrees)[1]
[bzhang@Melipal terzan5]$ gtexpmap
The exposure maps generated by this tool are meant
to be used for *unbinned* likelihood analysis only.
Do not use them for binned analyses.
Event data file[ter5_filtered_5deg_gti.fits]
Spacecraft data file[L100618062133E0D2F37E56_SC00.fits]
Exposure hypercube file[expcube_me.fits]
output file name[expmap_10.fits]
Response functions[P6_V3_DIFFUSE]
Radius of the source region (in degrees)[10]
Number of longitude points (2:1000) [120]
Number of latitude points (2:1000) [120]
Number of energies (2:100) [20]
在最后的分析之前，还要解决一个比较麻烦的事情：为视场中的各个辐射源定义。辐射源分为两类，一是弥漫的背景源，二是点源，其中背景源又分为两类，河内与河外。背景源要好办一些，可以直接在费米望远镜的网站上下载，点源就有摸索的意味了，依本人的理解是将人为定义的参数与观测比较，然后看效果如何。当然，这个过程需要根据计算结果一再的调整参数。
所有的辐射源的模型都要写入一个xml文件中，这种文件有固定格式，也可以利用Fermi Science Tools自带的modeleditor工具来生成文件，只要按提示输入对应参数就可以了。以下是模型文件中背景源的定义示例，至于各个参数的含义，烦请各位查阅手册吧，挨个解释实在太麻烦了……
<?xml version=”1.0″ ?>
<source_library title=”source library”>
<source name=”EG_v02″ type=”DiffuseSource”>
<spectrum file=”../isotropic_iem_v02.txt” type=”FileFunction”>
<parameter free=”1″ max=”1000″ min=”1e-05″ name=”Normalization” scale=”1″ value=”5.39258″ />
</spectrum>
<spatialModel type=”ConstantValue”>
<parameter free=”0″ max=”10.0″ min=”0.0″ name=”Value” scale=”1.0″ value=”1.0″/>
</spatialModel>
</source>
然后是试刀对象Terzan 5的辐射谱模型，定义为带指数cutoff的幂律谱：
<source name=”Terzan 5″ type=”PointSource”>
<!– point source units are cm^-2 s^-1 MeV^-1 –>
<spectrum type=”PLSuperExpCutoff”>
<parameter free=”1″ max=”1000″ min=”1e-05″ name=”Prefactor” scale=”1e-08″ value=”0.125358″/>
<parameter free=”1″ max=”0″ min=”-5″ name=”Index1″ scale=”1″ value=”-1.70573″/>
<parameter free=”0″ max=”1000″ min=”50″ name=”Scale” scale=”1″ value=”100″/>
<parameter free=”1″ max=”30000″ min=”500″ name=”Cutoff” scale=”1″ value=”3000.18″/>
<parameter free=”0″ max=”5″ min=”0″ name=”Index2″ scale=”1″ value=”1″/>
</spectrum>
<spatialModel type=”SkyDirFunction”>
<parameter free=”0″ max=”360.” min=”-360.” name=”RA” scale=”1.0″ value=”267.02″/>
<parameter free=”0″ max=”90.” min=”-90.” name=”DEC” scale=”1.0″ value=”-24.7792″/>
</spatialModel>
</source>
其他辐射源的定义依此类推，反正形式差不多。用自带的fv绘制Terzan 5的待分析光谱如下，嗯，光子计数谱，凑合看吧……
这样就可以进行似然分析了：
[bzhang@Melipal terzan5]$ gtlike
Statistic to use (BINNED|UNBINNED) [UNBINNED]
Spacecraft file[L100618062133E0D2F37E56_SC00.fits]
Event file[ter5_filtered_5deg_gti.fits]
Unbinned exposure map[expmap_10.fits]
Exposure hypercube file[expcube_me.fits]
Source model file[terzan5model.xml]
这一步的输出结果很长，包括一系列统计参数外加置信区间。针对Terzan 5得出的结果如下，其他源就不一一列举了：
Terzan 5:
Prefactor: 0.117792 +/- 0.065512
Index1: -1.7828 +/- 0.280544
Scale: 100
Cutoff: 4244.11 +/- 2149.74
Index2: 1
Npred: 506.68
ROI distance: 0
TS value: 144.814
这个TS value是置信度的表征，遵循康普顿伽玛射线天文台数据处理的定义，细节可以参考Mattox J. R. et al., 1996, ApJ, 461, 396等文。另外输出结果的最后还会输出理论与实测光子数的比较，这也可以看作是拟合好坏的反映：如果差异过大，肯定要调整参数（这次试刀的结果顶多只能算马马虎虎，模型预测光子数65341.2，实际62536，没有细调，知道大概过程就好）。更详细的谱分析可以借助HEAsoft中的XSPEC来完成，这就不是本文的范围所及了。如果是经常作数据处理，一个建议是将所有步骤写成脚本调用，操作起来要方便一些。谱处理之外，这次还接触了一下脉冲星的分析，不过似乎Fermi Science Tools在这此有些bug，建议还是使用官方推荐的其他小组开发的软件^_^
本人还没有接触到一个更加消耗时间的事情：搜索新源，一般是在光谱分析完成后进行的。这种活计需要生成TS map，据说那是可以轻易花掉好几天的事情，回来有机会再说吧……
嗯嗯，倒是希望今年下半年能把HEAsoft学会，自己先祝自己下月相关事宜顺利~~
