或者: http://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/gbm_grb_analysis.html
这里是想把NaI和BGO两个仪器的数据拿来联合分析. 以GRB 080810为例, 步骤如下:
1. 获取数据
位置:
http://heasarc.gsfc.nasa.gov/FTP/fermi/data/gbm/
也可以用ftp下载.
注意实际上Fermi提供更方便的方式, 登录它的服务器, 这样就不用下载数据, 但还是要安装软件. 通过软件链接服务器的. (貌似还可以用另外一种方式登录, 连软件都不用安装.) 另外, 数据的获取还可以通过这里: http://heasarc.gsfc.nasa.gov/cgi-bin/W3Browse/w3table.pl?Action=Detailed%20Mission&Observatory=fermi
这个应该是更官方的.
2. 进入目录
~/data/triggers/2008/bn080810549/current$ ls
glg_cspec_b0_bn080810549_v02.pha glg_ctime_n2_bn080810549_v02.pha
glg_cspec_b0_bn080810549_v02.rsp glg_ctime_n2_bn080810549_v02.rsp
glg_cspec_b1_bn080810549_v02.pha glg_ctime_n3_bn080810549_v02.pha
glg_cspec_b1_bn080810549_v02.rsp glg_ctime_n3_bn080810549_v02.rsp
glg_cspec_n0_bn080810549_v02.pha glg_ctime_n4_bn080810549_v02.pha
glg_cspec_n0_bn080810549_v02.rsp glg_ctime_n4_bn080810549_v02.rsp
glg_cspec_n1_bn080810549_v02.pha glg_ctime_n5_bn080810549_v02.pha
glg_cspec_n1_bn080810549_v02.rsp glg_ctime_n5_bn080810549_v02.rsp
glg_cspec_n2_bn080810549_v02.pha glg_ctime_n6_bn080810549_v02.pha
glg_cspec_n2_bn080810549_v02.rsp glg_ctime_n6_bn080810549_v02.rsp
glg_cspec_n3_bn080810549_v02.pha glg_ctime_n7_bn080810549_v02.pha
glg_cspec_n3_bn080810549_v02.rsp glg_ctime_n7_bn080810549_v02.rsp
glg_cspec_n4_bn080810549_v02.pha glg_ctime_n8_bn080810549_v02.pha
glg_cspec_n4_bn080810549_v02.rsp glg_ctime_n8_bn080810549_v02.rsp
glg_cspec_n5_bn080810549_v02.pha glg_ctime_n9_bn080810549_v02.pha
glg_cspec_n5_bn080810549_v02.rsp glg_ctime_n9_bn080810549_v02.rsp
glg_cspec_n6_bn080810549_v02.pha glg_ctime_na_bn080810549_v02.pha
glg_cspec_n6_bn080810549_v02.rsp glg_ctime_na_bn080810549_v02.rsp
glg_cspec_n7_bn080810549_v02.pha glg_ctime_nb_bn080810549_v02.pha
glg_cspec_n7_bn080810549_v02.rsp glg_ctime_nb_bn080810549_v02.rsp
glg_cspec_n8_bn080810549_v02.pha glg_tcat_all_bn080810549_v09.fit
glg_cspec_n8_bn080810549_v02.rsp glg_trigdat_all_bn080810549_v02.fit
glg_cspec_n9_bn080810549_v02.pha glg_tte_b0_bn080810549_v02.fit
glg_cspec_n9_bn080810549_v02.rsp glg_tte_b1_bn080810549_v02.fit
glg_cspec_na_bn080810549_v02.pha glg_tte_n0_bn080810549_v02.fit
glg_cspec_na_bn080810549_v02.rsp glg_tte_n1_bn080810549_v02.fit
glg_cspec_nb_bn080810549_v02.pha glg_tte_n2_bn080810549_v02.fit
glg_cspec_nb_bn080810549_v02.rsp glg_tte_n3_bn080810549_v02.fit
glg_ctime_b0_bn080810549_v02.pha glg_tte_n4_bn080810549_v02.fit
glg_ctime_b0_bn080810549_v02.rsp glg_tte_n5_bn080810549_v02.fit
glg_ctime_b1_bn080810549_v02.pha glg_tte_n6_bn080810549_v02.fit
glg_ctime_b1_bn080810549_v02.rsp glg_tte_n7_bn080810549_v02.fit
glg_ctime_n0_bn080810549_v02.pha glg_tte_n8_bn080810549_v02.fit
glg_ctime_n0_bn080810549_v02.rsp glg_tte_n9_bn080810549_v02.fit
glg_ctime_n1_bn080810549_v02.pha glg_tte_na_bn080810549_v02.fit
glg_ctime_n1_bn080810549_v02.rsp glg_tte_nb_bn080810549_v02.fit
3. 查看哪些NaI的仪器触发了
~/data/triggers/2008/bn080810549/current$ fdump glg_tcat_all_bn080810549_v09.fit
(然后一路回车, 就有下面的信息了. 注意红字体的, 1表示触发了, 我们就需要那个对应的fits文件. 注意从0开始数数: 0123456789ab. 另外Trigtime也是有用的.)
Name of optional output file[STDOUT]
Names of columns[]
Lists of rows[-]
SIMPLE = T / file does conform to FITS standard
BITPIX = 8 / number of bits per data pixel
NAXIS = 0 / number of data axes
COMMENT FITS (Flexible Image Transport System) format is defined in 'Astronomy
COMMENT and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H
CREATOR = 'GBM_TCAT_Regen.pl v1.2' / Software and version creating file
FILETYPE= 'TRIGGER ENTRY' / Name for this type of FITS file
TELESCOP= 'GLAST ' / Name of mission/satellite
INSTRUME= 'GBM ' / Specific instrument used for observation
OBSERVER= 'Meegan ' / GLAST Burst Monitor P.I.
ORIGIN = 'GIOC ' / Name of organization making file
DATE = '2009-06-24T17:12:14' / file creation date (YYYY-MM-DDThh:mm:ss UT)
DATE-OBS= '2008-08-10T13:08:02' / Date of start of observation
DATE-END= '2008-08-10T13:18:08' / Date of end of observation
TIMESYS = 'TT ' / Time system used in time keywords
TIMEUNIT= 's ' / Time since MJDREF, used in TSTART and TSTOP
MJDREFI = 51910 / MJD of GLAST reference epoch, integer part
MJDREFF = 7.428703703703703D-4 / MJD of GLAST reference epoch, fractional part
TSTART = 240066482.379032 / Time of first TRIGDAT INFO packet
TSTOP = 240067088.594364 / Time of last TRIGDAT INFO packet
FILENAME= 'glg_tcat_all_bn080810549_v09.fit' / Name of this file
TRIGTIME= 240066613.580576 / Trigger time relative to MJDREF, double precisi
OBJECT = 'GRB080810549' / Burst name in standard format, yymmddfff
More?[Yes]
RADECSYS= 'FK5 ' / Stellar reference frame
EQUINOX = 2000.0 / Equinox for RA and Dec
RA_OBJ = 356.8 / Calculated RA of burst
DEC_OBJ = 0.32 / Calculated Dec of burst
ERR_RAD = 0. / Calculated Location Error Radius
CLASS = 'GRB ' / Classification of trigger
RELIABLT= 0.5412 / Reliability of classification.
TRIGSCAL= 256 / [ms] Triggered timescale
TRIG_ALG= 9 / Triggered algorithm number
CHAN_LO = 3 / Trigger channel: low
CHAN_HI = 4 / Trigger channel: high
ADC_LO = 259 / Trigger channel: low (ADC: 0 - 4095)
ADC_HI = 1352 / Trigger channel: high (ADC: 0 - 4095)
DET_MASK= '000000011000' / Triggered detectors: (0-11)
CHECKSUM= 'PcogSaoZPaofPaoZ' / HDU checksum updated 2009-06-24T17:43:56
DATASUM = ' 0' / data unit checksum updated 2008-08-10T15:31:07
OBJ_CLAS= 'GRB ' / Classification of trigger
LOC_SRC = 'Swift ' / Mission/Instrument providing the localization
GEO_LONG= 999. / [deg] Spacecraft geographical east longitude
GEO_LAT = 999. / [deg] Spacecraft geographical north longitude
RA_SCX = 999. / [deg] Pointing of spacecraft x-axis: RA
DEC_SCX = 999. / [deg] Pointing of spacecraft x-axis: Dec
RA_SCZ = 999. / [deg] Pointing of spacecraft z-axis: RA
More?[Yes]
DEC_SCZ = 999. / [deg] Pointing of spacecraft z-axis: Dec
END
从上面看, 仪器7和8有比较好的数据, 那我们就选7吧.
4. 对数据进行bin, 产生一个有光变的文件
~/data/triggers/2008/bn080810549/current$ gtbin
This is gtbin version ScienceTools-v9r15p2-fssc-20090808
Type of output file (CCUBE|CMAP|LC|PHA1|PHA2) [PHA2] lc
Event data file name[glg_tte_b1_bn090323002_v01.fit] glg_tte_n7_bn080810549_v02.fit
Output file name[glg_tte_b1_bn090323002_v01.pha] glg_tte_n7_bn080810549_v02.lc
Spacecraft data file name[none]
Algorithm for defining time bins (FILE|LIN|SNR) [FILE] lin
Start value for first time bin in MET[259459340.241652] indef
Stop value for last time bin in MET[259459665.038936] indef
Width of linearly uniform time bins in MET[0.25]
注意后面是输入的内容. 软件总是记住前面输入过的, 如果保持不变的话, 就直接回车.
5. 换算成容易看的trigger time
~/data/triggers/2008/bn080810549/current$ fcalc
Name of FITS file and [ext#][glg_tte_b1_bn090323002_v01.lc] glg_tte_n7_bn080810549_v02.lc
Name of output FITS file[glg_tte_b1_bn090323002_v01_trig.lc] glg_tte_n7_bn080810549_v02_trig.lc
Resultant column name[TRIGTIME]
Arithmetic expression[Time-259459364.627422] Time-240066613.580576
注意前面fdump查看到的trigtime这里就用上了.
6. 可以查看一下光变
知道从哪里到哪里是信号, 大致哪里是背景, 便于选择背景.(这一步可以用fv, 图形话的界面, 来完成. 其实很多查看类的都可以用fv.)
~/data/triggers/2008/bn080810549/current$ fplot
Name of FITS file and [ext#][glg_tte_n7_bn080810549_v02_trig.lc]
Name of X Axis Parameter[error][TrigTime]
Name of Y Axis Parameter[error] up to 8 allowed[COUNTS]
Lists of rows[-]
Device: /XWindow, /XTerm, /TK, /PS, etc[/ps] /xw
Any legal PLT command[/XWindow]
Illegal command. At token = /XWINDOW
Type HELP to get command list.
PLT> quit
于是可以看到下面这个图. 大致上信号到70s的样子.
7. 制作time_bin文件
生成一个文件, 比如 glg_tte_n7_bn080810549_v02_time_bin.txt
内容如下
240066613.580576 240066683.580576
240066713.580576 240066753.580576
然后把它变成fit格式:
~/data/triggers/2008/bn080810549/current$ gtbindef
This is gtbindef version ScienceTools-v9r15p2-fssc-20090808
Type of bins, E -- energy or T -- time[T]
File containing the ascii bin definitions[n9_b1_time_bin.txt] glg_tte_n7_bn080810549_v02_time_bin.txt
Output file name[n9_b1_time_bin.fit] glg_tte_n7_bn080810549_v02_time_bin.fit
8. 制作最终拟合谱需要的pha
~/data/triggers/2008/bn080810549/current$ gtbin
This is gtbin version ScienceTools-v9r15p2-fssc-20090808
Type of output file (CCUBE|CMAP|LC|PHA1|PHA2) [LC] pha2
Event data file name[glg_tte_n7_bn080810549_v02.fit] glg_tte_n7_bn080810549_v02.fit
Output file name[glg_tte_n7_bn080810549_v02.lc] glg_tte_n7_bn080810549_v02.pha
Spacecraft data file name[none]
Algorithm for defining time bins (FILE|LIN|SNR) [LIN] file
Name of the file containing the time bin definition[n9_b1_time_bin.fit] glg_tte_n7_bn080810549_v02_time_bin.fit
如果只用NaI的数据拟合谱的话, 现在就可以了. 可以进入xpsec阶段. 不过我们要BGO的数据. 现在BGO有两个文件, 如何选择那个需要呢? 可以看看它们的光变(个人认为的), 看那个比较好比较强就用哪个. 那么就相当于重复上面的步骤.
9. 检查b0和b1的光变
~/data/triggers/2008/bn080810549/current$ gtbin
This is gtbin version ScienceTools-v9r15p2-fssc-20090808
Type of output file (CCUBE|CMAP|LC|PHA1|PHA2) [PHA2] lc
Event data file name[glg_tte_n7_bn080810549_v02.fit] glg_tte_b0_bn080810549_v02.fit
Output file name[glg_tte_n7_bn080810549_v02.pha] glg_tte_b0_bn080810549_v02.lc
Spacecraft data file name[none]
Algorithm for defining time bins (FILE|LIN|SNR) [FILE] lin
Start value for first time bin in MET[240066584.909792] indef
Stop value for last time bin in MET[240066914.41369] indef
Width of linearly uniform time bins in MET[0.25]
~/data/triggers/2008/bn080810549/current$ gtbin
This is gtbin version ScienceTools-v9r15p2-fssc-20090808
Type of output file (CCUBE|CMAP|LC|PHA1|PHA2) [LC]
Event data file name[glg_tte_b0_bn080810549_v02.fit] glg_tte_b1_bn080810549_v02.fit
Output file name[glg_tte_b0_bn080810549_v02.lc] glg_tte_b1_bn080810549_v02.lc
Spacecraft data file name[none]
Algorithm for defining time bins (FILE|LIN|SNR) [LIN]
Start value for first time bin in MET[240066584.9089] indef
Stop value for last time bin in MET[240066914.415662] indef
Width of linearly uniform time bins in MET[0.25]
然后可以用fv查看
fv glg_tte_b0_bn080810549_v02.lc
如上图, 首先是出现左边和上面两个窗口, 点第一个plot, 出现右边的窗口, 将time设为X, counts设为Y, 点plot, 就得到中间这个图. 还可以具体调解.另外图片还可以打印下来, 或者保存成ps(也是点打印).
得到这个. 选择"file"可以保存成文件.注意有个小技巧, 就是要将窗口拖动一下, 拖到屏幕的下面一点点, 然后用鼠标将窗口上面的边网上拉, 这样那些文字才能完整地显示出来. 点print就可以了. 默认打印出来的文件名是:powGraph.ps .
再看第二个文件b1的光变:
要稍微好点, 那就用这个吧.从光变上看, b0和b1都没有什么很强的光变. 意味着它们对谱的限制也不会很好.
10. 把b1转换到trig时间 (这一步不必要)
~/data/triggers/2008/bn080810549/current$ fcalc
Name of FITS file and [ext#][glg_tte_n7_bn080810549_v02.lc] glg_tte_b1_bn080810549_v02.lc
Name of output FITS file[glg_tte_n7_bn080810549_v02_trig.lc] glg_tte_b1_bn080810549_v02_trig.lc
Resultant column name[TRIGTIME]
Arithmetic expression[Time-240066613.580576]
11. 用同样的time bin生成谱
~/data/triggers/2008/bn080810549/current$ gtbin
This is gtbin version ScienceTools-v9r15p2-fssc-20090808
Type of output file (CCUBE|CMAP|LC|PHA1|PHA2) [LC] pha2
Event data file name[glg_tte_b1_bn080810549_v02.fit]
Output file name[glg_tte_b1_bn080810549_v02.lc] glg_tte_b1_bn080810549_v02.pha
Spacecraft data file name[none]
Algorithm for defining time bins (FILE|LIN|SNR) [LIN] file
Name of the file containing the time bin definition[glg_tte_n7_bn080810549_v02_time_bin.fit]
12. 用xspec拟合
这里就不一步一步的做了, 每次都这样做比较麻烦. 可以写一个相当于批处理的文件, 然后一次做完就可以了. 文件内容如下(文件名n7_b1_spectrum.xcm, 注意#表示说明, xspec就不理它们了.):
#This is the batch file to excute the fitting commands
data glg_tte_n7_bn080810549_v02.pha{1}, glg_tte_b1_bn080810549_v02.pha{1}
back glg_tte_n7_bn080810549_v02.pha{2}, glg_tte_b1_bn080810549_v02.pha{2}
resp glg_cspec_n7_bn080810549_v02.rsp, glg_cspec_b1_bn080810549_v02.rsp
ignore 1:1-3,125-**
ignore 2:1-5,121-**
model grbm & /*
#cutoffpl, grbm
#indicate all the prompted queries using the default values
fit 1000
setplot energy
#cpd spectrum_n9_na_res.eps /cps
#plot data res
cpd spectrum_n7_b1.eps /cps
plot ufspec
最终结果如下, 注意红体字部分是主要结果.
~/data/triggers/2008/bn080810549/current$ xspec n7_b1_spectrum.xcm
XSPEC version: 12.5.0ac
Build Date/Time: Thu Apr 2 03:38:12 2009
2 spectra in use
Spectral Data File: glg_tte_n7_bn080810549_v02.pha{1} Spectrum 1
Net count rate (cts/s) for Spectrum:1 1.285e+03 +/- 4.292e+00
Assigned to Data Group 1 and Plot Group 1
Noticed Channels: 1-128
Telescope: GLAST Instrument: GBM Channel Type: PHA
Exposure Time: 69.73 sec
No response loaded.
Spectral Data File: glg_tte_b1_bn080810549_v02.pha{1} Spectrum 2
Net count rate (cts/s) for Spectrum:2 1.469e+03 +/- 4.591e+00
Assigned to Data Group 1 and Plot Group 2
Noticed Channels: 1-128
Telescope: GLAST Instrument: GBM Channel Type: PHA
Exposure Time: 69.7 sec
No response loaded.
***Warning! One or more spectra are missing responses,
and are not suitable for fit.
Net count rate (cts/s) for Spectrum:1 1.779e+02 +/- 6.797e+00 (13.9 % total)
Net count rate (cts/s) for Spectrum:2 1.061e+02 +/- 7.436e+00 (7.2 % total)
***Warning! One or more spectra are missing responses,
and are not suitable for fit.
***Warning: No TLMIN keyword value for response matrix FCHAN column.
Will assume TLMIN = 1.
Response successfully loaded.
***Warning: No TLMIN keyword value for response matrix FCHAN column.
Will assume TLMIN = 1.
Response successfully loaded.
3 channels (1,3) ignored in spectrum # 1
4 channels (125,128) ignored in spectrum # 1
5 channels (1,5) ignored in spectrum # 2
8 channels (121,128) ignored in spectrum # 2
========================================================================
Model grbm<1> Source No.: 1 Active/On
Model Model Component Parameter Unit Value
par comp
1 1 grbm alpha -1.00000 +/- 0.0
2 1 grbm beta -2.00000 +/- 0.0
3 1 grbm tem keV 300.000 +/- 0.0
4 1 grbm norm 1.00000 +/- 0.0
________________________________________________________________________
Chi-Squared = 2.669440e+07 using 236 PHA bins.
Reduced chi-squared = 115062.1 for 232 degrees of freedom
Null hypothesis probability = 0.000000e+00
Current data and model not fit yet.
Chi-Squared Lvl Par # 1 2 3 4
273.146 -1 -0.990225 -2.72208 380.86 0.00320712
272.375 -1 -1.03279 -5.57841 447.841 0.00302134
272.367 1 -1.03304 -4.9986 447.832 0.00302167
==================================================
Variances and Principal Axes
1 2 3 4
6.7881E-08| -0.0001 -0.0000 -0.0000 -1.0000
4.2610E-03| 1.0000 0.0002 0.0007 -0.0001
3.6926E+04| -0.0006 -0.0507 0.9987 -0.0000
2.5564E+03| 0.0002 -0.9987 -0.0507 0.0000
--------------------------------------------------
================================================
Covariance Matrix
1 2 3 4
1.951e-02 5.800e-01 -2.361e+01 8.493e-05
5.800e-01 2.645e+03 -1.741e+03 4.364e-03
-2.361e+01 -1.741e+03 3.684e+04 -1.326e-01
8.493e-05 4.364e-03 -1.326e-01 5.466e-07
------------------------------------------------
========================================================================
Model grbm<1> Source No.: 1 Active/On Model
Model Component Parameter Unit Value
par comp
1 1 grbm alpha -1.03304 +/- 0.139667
2 1 grbm beta -4.99860 +/- 51.4278
3 1 grbm tem keV 447.832 +/- 191.931
4 1 grbm norm 3.02167E-03 +/- 7.39321E-04 ________________________________________________________________________
Chi-Squared = 272.37 using 236 PHA bins.
Reduced chi-squared = 1.1740 for 232 degrees of freedom
Null hypothesis probability = 3.548480e-02
Do you really want to exit? (y)
XSPEC: quit
拟合的谱如图:
正如 http://glast-ground.slac.stanford.edu/workbook/pages/sciTools_latGbmJointSpectralAnalysis/sciTools_latGbmJointSpectralAnalysis.htm 最后所说, 联合拟合的谱的图并不直观. (看起来好像拟合的不好样.)
也可以换个模型拟合试试看. 比如把上面的grbm(实际上就是band谱)换成cutoffpl试试.
主要结果如下:
========================================================================
Model cutoffpl<1> Source No.: 1 Active/On
Model Model Component Parameter Unit Value
par comp
1 1 cutoffpl PhoIndex 1.10457 +/- 0.134432
2 1 cutoffpl HighECut keV 499.988 +/- 231.849
3 1 cutoffpl norm 0.466366 +/- 0.231280
________________________________________________________________________
Chi-Squared = 272.23 using 236 PHA bins.
Reduced chi-squared = 1.1684 for 233 degrees of freedom
Null hypothesis probability = 3.965676e-02
误差Reduced chi-squared = 1.1684, 看起来也挺好的. 说明高能部分限制太松了, band谱和cutoff power law都可以拟合.
