宇宙X線スペクトル解析向けのxspecの使い方
この記事は、xspecの基礎はある程度は知っていて、bashスクリプトと連携する方法なども知りたい人向けに書いています。
xspecの基本については、榎戸さんの記事、
も合わせてご参照ください。
xspec manual を全部読むのは膨大ですが、
のページは一読しておくと、何をやるツールなのかがわかって良いです。
はじめに
HEASOFT をまだインストールしていない人は、
を参考にしてください。
ここではHEASOFTの設定ができていると想定します。まずは、xspec と打った時に、起動するxspecの本体を確認しましょう。
which というコマンドで、full path を教えてくれます。
$ which xspec
> /Users/syamada/work/software/heasoft-6.27.2/x86_64-apple-darwin18.7.0/bin/xspec
のように、自分が意図した xpsec が起動していることを確認しましょう。
それができたら、次に、
$ xspec
dummyrsp 0.1 20 100000
model apec
2
1.0
1.0
1
pl eemo
cpd /xs
を実行して、 apec というモデルが表示されることを確認してください。
実際に実行する風景は、こんな感じになります。
$ xspec
XSPEC version: 12.11.0
Build Date/Time: Tue Aug 4 22:49:16 2020
XSPEC12>dummyrsp 0.1 20 100000
XSPEC12>model apec
Input parameter value, delta, min, bot, top, and max values for ...
1 0.01( 0.01) 0.008 0.008 64 64
1:apec:kT>2
1 -0.001( 0.01) 0 0 5 5
2:apec:Abundanc>1.0
0 -0.01( 0.01) -0.999 -0.999 10 10
3:apec:Redshift>1.0
1 0.01( 0.01) 0 0 1e+20 1e+24
4:apec:norm>1
Reading APEC data from 3.0.9
========================================================================
Model apec<1> Source No.: 1 Active/Off
Model Model Component Parameter Unit Value
par comp
1 1 apec kT keV 2.00000 +/- 0.0
2 1 apec Abundanc 1.00000 frozen
3 1 apec Redshift 1.00000 frozen
4 1 apec norm 1.00000 +/- 0.0
________________________________________________________________________
XSPEC12>pl eemo
XSPEC12>cpd /xs
このように、手で一つ一つコマンドを打つことができますが、スペクトル解析においては、モデルやファイルが変わることが多いので、その度に手でコマンドを打っていては非効率です。自動化するためには、ヒアドキュメント bashのヒアドキュメントを活用する を使います。下記、その例を紹介します。
dummyrsp でモデルをプロットする簡単な例
dummyrsp という対角レスポンスを生成するコマンドで，疑似的な装置応答を生成して，モデルを表示する練習です．
見よう見まねで動かしてみるのが第一歩なので，まずは次のコマンドを打ってみよう．bashスクリプトとして動作させてもよいし，xspecと打ってからコマンドを入れる場合は，xspec << EOFの間をコピペで入れればよいです．
#!/bin/sh
rm -f test1.* # remove outfiles
xspec <<EOF
dummyrsp 0.1 20 100000
model apec
2
1.0
1.0
1
log test1.log
flux 0.5 2
log none
pl eemo
cpd /xs
iplot
r y 1e-5 1e3
r x 0.1 20
la t 0.5 keV
hard test1.ps/cps
we test1
exit
exit
EOF
ps2pdf test1.ps # convert ps file to pdf if needed.
と書かれたスクリプトを実行してみよう．
>chmod +x mk_apec.sh # chmod +x で実行権限を与える
>./mk_apec.sh
>bash ./mk_apec.sh # これでもよい
権限を与えて動かしてみよう．
ただし，mac or linux のターミナルで，xspec に path が通っていることが前提です．
これがうまく走ると，
このような図が生成される．
正しく動くと，下記のようなメッセージが出力されるはずである．
XSPEC version: 12.11.0
Build Date/Time: Tue Aug 4 22:49:16 2020
XSPEC12>dummyrsp 0.1 20 100000
XSPEC12>model apec
Input parameter value, delta, min, bot, top, and max values for ...
1 0.01( 0.01) 0.008 0.008 64 64
1:apec:kT>2
1 -0.001( 0.01) 0 0 5 5
2:apec:Abundanc>1.0
0 -0.01( 0.01) -0.999 -0.999 10 10
3:apec:Redshift>1.0
1 0.01( 0.01) 0 0 1e+20 1e+24
4:apec:norm>1
Reading APEC data from 3.0.9
========================================================================
Model apec<1> Source No.: 1 Active/Off
Model Model Component Parameter Unit Value
par comp
1 1 apec kT keV 2.00000 +/- 0.0
2 1 apec Abundanc 1.00000 frozen
3 1 apec Redshift 1.00000 frozen
4 1 apec norm 1.00000 +/- 0.0
________________________________________________________________________
XSPEC12>log test1.log
Logging to file:test1.log
XSPEC12>flux 0.5 2
Model Flux 0.18586 photons (2.3763e-10 ergs/cm^2/s) range (0.50000 - 2.0000 keV)
XSPEC12>log none
Log file closed
logging switched off
XSPEC12>pl eemo
XSPEC12>cpd /xs
XSPEC12>iplot
PLT>
PLT> r y 1e-5 1e3
PLT> r x 0.1 20
PLT> la t 0.5 keV
PLT> hard test1.ps/cps
PLT> we test1
PLT>
PLT> exit
XSPEC12>
XSPEC12>exit
XSPEC: quit
解説
ファイルの初期化
rm -f test1.* # remove outfiles
xspec に限らず，X線のデータ解析ツールは書き出す先のファイルがすでにある場合は，「消しますか？」と聞いてくるので，スクリプトなどでその処理次第では止まってしまうことがあるので，生成する予定のファイルがすでにある場合は消しておく．rm -f の -f をつけることで，可否を確認せずに強制的(forcely)に消す，という意味になる．
ヒアドキュメントの使い方
xspec の後のEOFは、ヒアドキュメントと呼ばれるもので，bashのヒアドキュメントを活用する がわかりやすい．
xspec <<EOF
...
EOF
と書くことで，EOFで囲まれた部分が，xspecに読み込まれるというイメージになります．
xspec 2>&1 > all.log <<EOF
とすると，全部のメッセージをall.logに書き出すことができる．
擬似的なレスポンスの生成
dummyrsp 0.1 20 100000
xspec の dummyrsp というとっても便利なコマンドを使います．意味は，0.1 keV から 20 keV まで 100000 ビンで，対角レスポンスを生成せよ，という意味で，xspec に登録されているモデルを表示したい時や，flux を計算したいときによく使います．
スペクトルモデルの指定
スペクトルモデルを指定します．ここでは，熱的プラズマからの放射モデルである apec を使います．
model apec
2
1.0
1.0
1
の意味は，実行時のログをみるとわかりますが，
1:apec:kT>2
2:apec:Abundanc>1.0
3:apec:Redshift>1.0
4:apec:norm>1
のように，パラメータの初期値を与えたことに対応しています．
fluxの計算とログの設定
flux というコマンドで，0.5keVから2keVまでのfluxを計算して，それを test1.log というテキストファイルに中身を書き出します．これは，部分的な情報だけをファイルに書き出したい場合によく使います．
log test1.log
flux 0.5 2
log none
プロットする
pl eemo
これは，plot eemodel の略で，そのほかにも plot には色んなオプションがあるので，自由に選んでよいです．
ちなみに，eemoやeeufの"e"は"energy"の意味で，"eemo" は，"energy" x "energy" x "光子スペクトル model "という意味です．宇宙X線の場合などは，エネルギーの -2乗などでスペクトルが下がることが多いので，見やすさのために光子スペクトルの絵エネルギーを２回掛けることが多いです．それをlog-log表示したときの積分値が放射エネルギーになる，という理由もあります．
表示
X window で表示させます．macでQuartsの設定などで動かない場合は，昔のX11が残っているとかそういう場合かと思います．
cpd /xs
QDPでプロット
xspecは，QDPという若い人は使わないほうがよい古いプロットツールがbuilt-inされていて，これを使うしかないです．
iplot
r y 1e-5 1e3 # Y軸の範囲の設定
r x 0.1 20 # X軸の範囲の設定
la t 0.5 keV # titile の設定
hard test1.ps/cps # psファイルに書き出す
we test1 # qdp, pco ファイルを生成する
exit
はxspecのルールではなくて，QDPのルールで書かれたものです．iplot という xspec のコマンドで，QDPに移行します．
横軸縦軸の変更は，
log x
log y
log x off
log y off
で変更する．
psファイルをpdfに変換
これはなくてもよいですが，pdf に変換した方が便利な人はやってください．
ps2pdf test1.ps # convert ps file to pdf if needed.
NICERを例に実データでシミュレーションする例
データ取得
から，
- nixtiback20190807.pi : バックグランドのスペクトル
- nixtiref20170601v002.rmf : 検出器のエネルギー分解能，Energy vs. PHA (パルスハイト)の関係，検出器の応答関数の入ったファイル
- nixtiaveonaxis20170601v004.arf : 望遠鏡の有効面積
の３つのファイルをダウンロードする．rmf と arf の住み分け方は，X線固有であるが，rmf には有効面積の情報は含めずに検出効率(入力したエネルギーによる出力のエネルギーの行列)を入れて，望遠鏡の有効面積はarfに入れるという使い分けをしている．
シミュレーション時には，バックグランドファイルが普通は必要となる．なくてもよいが，何かしらPHAファイルというスペクトルファイルが必要になる．
xspecのコマンドの使い方，およびfakeの仕方
xspecの使い方も含めて羅列した． # 以下はここでのコメントなので，実際には打たなくてよい．
> xspec
XSPEC12>data 1:1 nixtiback20190807.pi
Error: cannot read response file nixtiref20170601v001.rmf
# piファイルにrmfの名前が書き込まれているがこれがないと言われるだけ．
# rmfファイルが更新された場合になどによくあるエラー．次でrmfファイルを読み込むので無視してよい．
New filename ( "none" or "/*" to return to the XSPEC prompt): none
... skipped
response file read skipped
1 spectrum in use
Spectral Data File: nixtiback20190807.pi Spectrum 1
Net count rate (cts/s) for Spectrum:1 8.241e+00 +/- 2.871e-03
...
XSPEC12>resp 1 nixtiref20170601v002.rmf
Response successfully loaded.
XSPEC12>arf 1 nixtiaveonaxis20170601v004.arf
Arf successfully loaded.
# XSPEC12>back 1 nixtiback20190807.pi バックグランドを含めたシミュレーションの場合はここでバックグランドを入れる．
XSPEC12>cpd /xs
XSPEC12>pl eff # 有効面積の確認
XSPEC12>setp e # エネルギー軸に変更
XSPEC12>pl ld # バックグランドのスペクトルの確認
XSPEC12>model wabs * power # 適当にパワーローのモデルを入れる．
Input parameter value, delta, min, bot, top, and max values for ...
1 0.001( 0.01) 0 0 100000 1e+06
1:wabs:nH>10
1 0.01( 0.01) -3 -2 9 10
2:powerlaw:PhoIndex>2.1
1 0.01( 0.01) 0 0 1e+20 1e+24
3:powerlaw:norm>0.1
========================================================================
Model wabs<1>*powerlaw<2> Source No.: 1 Active/On
Model Model Component Parameter Unit Value
par comp
1 1 wabs nH 10^22 10.0000 +/- 0.0
2 2 powerlaw PhoIndex 2.10000 +/- 0.0
3 2 powerlaw norm 0.100000 +/- 0.0
________________________________________________________________________
XSPEC12>fakeit # シミュレーションを行う
Use counting statistics in creating fake data? (y):
Input optional fake file prefix: testpl_
Fake data file name (testpl_nixtiback20190807.fak):
Exposure time, correction norm, bkg exposure time (1.00000e+06, -1.00000, 1.00000): 1e5
No background will be applied to fake spectrum #1
1 spectrum in use
XSPEC12>pl ld delc # 残差を確認
XSPEC12>ignore **-0.5 10.0-** # 0.5-10 keV のみ取り出す． 10.0ではなく10とするとエネルギーではなく，チャンネルと判断されるので，エネルギーでカットする場合は実数で指定する．
50 channels (1-50) ignored in spectrum # 1
502 channels (1000-1501) ignored in spectrum # 1
Fit statistic : Chi-Squared 932.02 using 949 bins.
Test statistic : Chi-Squared 932.02 using 949 bins.
Null hypothesis probability of 6.21e-01 with 946 degrees of freedom
Current data and model not fit yet.
XSPEC12>pl eeuf # vFv でのプロットする．近似的に検出器応答を解く．
XSPEC12>save model pl.mo # モデルのみ保存
XSPEC12>save all pl.xcm # すべて保存
XSPEC12>show mo # モデルの確認
Current model list:
Model wabs<1>*powerlaw<2> Source No.: 1 Active/On
For Data Group(s): 1
Using energies from responses.
XSPEC12>editmo wabs * (powerlaw + gauss) # モデルの一つだけ修正，ここでは gauss を追加する．
Input parameter value, delta, min, bot, top, and max values for ...
6.5 0.05( 0.065) 0 0 1e+06 1e+06
4:gaussian:LineE>6.4
0.1 0.05( 0.001) 0 0 10 20
5:gaussian:Sigma>0.001
1 0.01( 0.01) 0 0 1e+20 1e+24
6:gaussian:norm>0.0001
Fit statistic : Chi-Squared 1451.00 using 949 bins.
Test statistic : Chi-Squared 1451.00 using 949 bins.
Null hypothesis probability of 2.89e-24 with 943 degrees of freedom
Current data and model not fit yet.
========================================================================
Model wabs<1>(powerlaw<2> + gaussian<3>) Source No.: 1 Active/On
Model Model Component Parameter Unit Value
par comp
1 1 wabs nH 10^22 10.0000 +/- 0.0
2 2 powerlaw PhoIndex 2.10000 +/- 0.0
3 2 powerlaw norm 0.100000 +/- 0.0
4 3 gaussian LineE keV 6.40000 +/- 0.0
5 3 gaussian Sigma keV 1.00000E-03 +/- 0.0
6 3 gaussian norm 1.00000E-04 +/- 0.0
________________________________________________________________________
XSPEC12>pl
XSPEC12>newp 6 1e-5 # パラメータを変更
XSPEC12>pl # 確認
XSPEC12>newp 6 1e-4 # 元に戻す．
Fit statistic : Chi-Squared 1451.00 using 949 bins.
Test statistic : Chi-Squared 1451.00 using 949 bins.
Null hypothesis probability of 2.89e-24 with 943 degrees of freedom
Current data and model not fit yet.
XSPEC12>iplot # QDPに入る．
***Warning: Fit is not current.
PLT> r x 6 7 # 6.0-7.0 keV まで拡大
PLT> log x off # X軸リニア
PLT> log y off # Y軸リニア
PLT> r y 0 0.3 #Y軸の範囲を 0 - 0.3 に指定
XSPEC12>show file # ファイルの確認
1 file 1 spectrum
Spectral Data File: testpl_nixtiback20190807.fak Spectrum 1
Net count rate (cts/s) for Spectrum:1 7.478e+00 +/- 8.647e-03
Assigned to Data Group 1 and Plot Group 1
Noticed Channels: 51-999
Telescope: NICER Instrument: XTI Channel Type: PI
Exposure Time: 1e+05 sec
Using fit statistic: chi
Using Response (RMF) File nixtiref20170601v002.rmf for Source 1
Using Auxiliary Response (ARF) File nixtiaveonaxis20170601v004.arf
観測シミュレーションで，観測時間を変更する場合には，
Exposure time, correction norm, bkg exposure time (1.00000e+06, -1.00000, 1.00000):
の時間を変更する． それで， pl ld delc eeuf ra など，xspecで残差を表示したり，モデルを変えて見たり，error や steppar コマンドで，エラーの範囲を見積るなど，普通の観測と同様の操作で可能となる．
vFvでの注意点
一つ注意点は，観測シミュレーションをするときに，エネルギー分解能よりも細いラインを入れた場合に vFv でみるとおかしなことになる．
この左は上記のケースで， wabs * power というパワーローに星間吸収の入ったモデルでシミュレーションした場合の例です． ここで， editmo により， ガウシアンを追加した場合にどうなるか？というと， vFv のモデルではなくて，データが輝線の周りで凹んで，尖っているような変な構造が出現してしまっている． どちらも同じデータで，モデルだけ変えて， plot eeuf ra しています．
鉄付近だけ拡大してみると，
このように，結構な違いが出てしまう． ちなみに， QDP で， log x off, log y off, r x 6 7, r y 0 0.3 として，このような拡大図を生成している．
これは， xspec の unfolded spectrum を作るときの仕様で，モデルが正しく見えるようにデータ点を操作してしまうために起こる．これを避けるためには，データに合わない尖った構造をもつモデルで vFv をしないとか，連続成分で合わせておいて， vFv をするなど，使うときには注意が必要となる．
「ufspec、eufspec、eeufspec」 の詳細な説明
より正確には、
の ufspec、eufspec、eeufspec の部分を読んでください。
簡単に説明すると、
plot eeuf (or plot euf, plot uf)などのコマンドは、モデルと unfolded スペクトルをプロットするもので、folded という言葉はジャーゴンで、検出器応答が畳み込まれたという意味で、fold という言い方をします。検出器応答が除かれた、という意味は、unfolded という言葉になります。
注意点は、このプロットはモデルに依存しており、モデルが変更されると unfolded model のポイントは移動します。unfolded model としてプロットされるデータポイントは、
f(E) = Data (E) * (unfolded~model(E) )/(folded~model (E) )
によって計算されます。
ここで、Data (E) は観測データ、unfolded~model(E)はプロットのエネルギービン幅で積分された理論モデルの値、folded~model (E)は応答を考慮した(標準プロットで表示される)モデルです。eufspecは、$Ef(E)$、または波長をプロットする場合は$\lambda f(\lambda)$でプロットします。eeufspecは、$E^{2}f(E)$、または波長をプロットする場合は$\lambda^{2}f(\lambda)$でプロットします。乗算因子で使用されるE（または$\lambda$）は、プロットのエネルギービンの下限エネルギーと上限エネルギーの幾何平均が使われます。
folded~model (E) = A(E) \times (unfolded~model(E))
であれば単純に、
f(E) = Data (E) / A(E)
となり、A(E)で決まるベクトル(いわゆる arf file)分だけ補正することになります。
問題となるのは、
folded~model (E) = \sum_i R(E, E_i) \times unfolded~model(E_i)
のように、単一のエネルギー E だけではなくて、周囲のエネルギー E_i に依存する場合です。
具体的に、unfolded~model(E) ~ δ(E-Ei) + 適当な連続成分、として Ei にピークを持つデルタ関数を考え、レスポンス R(E, E_i) がだいたいガウシアンを思うと、folded~model (E) は、Ei を中心として、エネルギー分解能程度で広がるガウシアンになります。
一方、unfolded~model (E) は、Eiを中心としたガウシアンを仮定しているので(unfolded~model(E) )/(folded~model (E) ) は、Ei で尖って、エネルギー分解能程度の幅までは小さい値になるような、上で示したような図になります。
エネルギー分解能程度に広げておけば、そのような変な構造が小さくなるので、「vFvの図を作るときには、ラインの幅を装置のエネルギー分解能程度を入れておきなさい」、というよく言われるルールの所以です。
シェルスクリプトを用いた簡単なプロット方法
xspecを使ったWHIMの簡単なシミュレーションを例として，シェルスクリプトとxspecの使い方を紹介する．
宇宙背景X線放射，LHB, MH
半径20'の円のfluxで足並みをそろえて， Hattori et al. 2017 のTable3の(i)のapec@z=0.3を入れる．
CXBについては，
https://heasarc.gsfc.nasa.gov/docs/suzaku/analysis/pin_cxb.html
に従うと，
CXB(E) = 9.412\times 10^{-3} \times (E/1keV)^{-1.29} \times \exp(-E/40keV) ~~~~ \rm{photons/cm2/s/HXDFOV/keV}
なので，これに，XISのFOV = (20 arcmin)^2 pi とHXDのFOV = 4deg^2 = (240 arcmin)^2 の割合をかけてあげて，
(20.203.14)/(240*240) * 9.4e-3 = 2.049e-4 photons/cm2/s/XISFOV/keV
と考える．
モデル
apec + TBabs(apec + powerlaw + apec) のモデルで，
apec[LWB] + TBabs(apec[MWH] + powerlaw[CXB] + apec[WHIM]) という意味で，ローカルホットバブル（LHB），天の川ハロー（MWH）とからの宇宙X線背景（CXB）であり，TBabsで表現される星間吸収は，近場にあるローカルホットバブル（LHB）には掛からない．
method leven 10 0.01
abund angr
xsect vern
cosmo 70 0 0.73
xset delta 0.01
systematic 0
model apec + TBabs(apec + powerlaw + apec)
0.34 0.01 0.008 0.008 64 64
1 -0.001 0 0 5 5
0 -0.01 -0.999 -0.999 10 10
0.000335 0.01 0 0 1e+20 1e+24
0.0139 0.001 0 0 100000 1e+06
0.13 0.01 0.008 0.008 64 64
1 -0.001 0 0 5 5
0 -0.01 -0.999 -0.999 10 10
0.00171 0.01 0 0 1e+20 1e+24
1.4 0.01 -3 -2 9 10
0.0002 0.01 0 0 1e+20 1e+24
0.17 0.01 0.008 0.008 64 64
0.2 -0.001 0 0 5 5
0.3 -0.01 -0.999 -0.999 10 10
0.00285 0.01 0 0 1e+20 1e+24
bayes off
QDPファイルの設定
col 0 on 1
col 2 on 2
col 3 on 3
col 4 on 5
col 5 on 4
lw 4 on 2
lw 4 on 3
lw 4 on 4
lw 4 on 5
LStyle 1 on 2
LStyle 1 on 3
LStyle 1 on 4
LStyle 1 on 5
fo ro
la ro
lw 4
la t "red:MWH green:LHB cyan:CXB blue:WHIM"
シェルスクリプト
seq 0 0.1 1 というコマンドで，0から1まで0.1刻みの数字列を生成し，for loop の変数 rshift に入れる．bash スクリプトは，for 変数 in 要素 で回して，do と done までの間をloopする仕様である．
#!/bin/sh
for rshift in `seq 0 0.1 1`
do
ftag=`echo $rshift | sed "s/\./p/g"`
fnamefig=whim_redshift${ftag}
echo $fnamefig
rm -f $fnamefig
xspec <<EOF
dummyrsp 0.1 3 3000
@background_whim.mo
newp 14 $rshift
plot eemo
iplot
@color.pco
r y 1e-8 1
r x 0.1 3
la t "red:MWH green:LHB cyan:CXB blue:WHIM redshift = ${rshift}"
hard ${fnamefig}.ps/cps
q
exit
EOF
ps2pdf ${fnamefig}.ps
done
実行方法
%chmod +x
./make_fig.sh
chmod で実行権限を与えて，実行する．
gif アニメ を生成するには，
convert -delay 50 -loop 5 *.pdf -alpha off checkwm.gif
とする．-alpha off をつけないと背景が透明になる．
実行結果
そのほか参考資料
新しいのはないですが，「すざく」ファーストステップガイドがよくまとまっている．(ただし，日本語版のすざくのマニュアルは途中で英語版オンリーになったので，すざくのデータ解析については英語版を参照すること)
見た目通りにビンまとめしたい場合はこちらの記事を参考にしてください。