宇宙X線スペクトル解析向けのxspecの使い方

この記事は、xspecの基礎はある程度は知っていて、bashスクリプトと連携する方法なども知りたい人向けに書いています。
xspecの基本については、榎戸さんの記事、

も合わせてご参照ください。

レスポンスやARFファイルについては、酒井さんの記事、

を参考にしてください。

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の間をコピペで入れればよいです．

mk_apec.sh

#!/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を例に実データでシミュレーションする例

データ取得

https://heasarc.gsfc.nasa.gov/docs/nicer/proposals/nicer_tools.html

から，

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$）は、プロットのエネルギービンの下限エネルギーと上限エネルギーの幾何平均が使われます。

単純に、有効面積A(E)という一次元のベクトルで

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) がだいたいガウシアンの Line Spread Function (LSF)を持つ仮定すると、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）には掛からない．

background_whim.mo

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する仕様である．

make_fig.sh

#!/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 をつけないと背景が透明になる．

実行結果

そのほか参考資料

新しいのはないですが，「すざく」ファーストステップガイドがよくまとまっている．(ただし，日本語版のすざくのマニュアルは途中で英語版オンリーになったので，すざくのデータ解析については英語版を参照すること)

見た目通りにビンまとめしたい場合はこちらの記事を参考にしてください。

XRISM衛星のレスポンス、ARFについては、下記を参考にしてください。