owned this note
owned this note
Published
Linked with GitHub
# XspecによるＸ線スペクトル解析
最終更新: 2020-11-07 (榎戸輝揚)
###### tags: `lecture`,`xspec`
## Xspec とは？
X線のスペクトル解析を行うパッケージはいくつかあり、その代表的な一つが Xspec です。
- [Xspec Home Page](https://heasarc.gsfc.nasa.gov/xanadu/xspec/) (公式ウェブサイト)
- [Xspec Manual (pdf)](https://heasarc.gsfc.nasa.gov/xanadu/xspec/XspecManual.pdf)
この Xspec は、NASA を中心に開発されてきた高エネルギー天文学のデータ解析の統合パッケージ "[HEASoft](https://heasarc.gsfc.nasa.gov/docs/software/heasoft/)" に含まれており、HEASoft をインストールする必要があります。
- [Installing HEASoft](https://heasarc.gsfc.nasa.gov/lheasoft/install.html) (HEASoft のインストールガイド)
に従って、インストールできます。（昔はいろいろハマるところがありましたが、最近は、書かれている通りにやればインストールできます）
インストールして、Xspec を呼び出すことができるようになると、コマンドラインで対話的にスペクトルのフィットを行うことができます。もちろん、バッチ処理にすることもできます。
## Ｘ線スペクトル解析の基本
以下、コマンドラインを使える人は実際に自分でやってもらえると良いのですが、そうでない方は以下に動画を作っていますので、これを見てもらえればと思います。
- [YouTube: Xspecでのスペクトル解析](https://youtu.be/aPe2kL-tN50)
まず練習用のディレクトリを作ります。
```
mkdir practice
cd practice
```
必要なファイル一式は、[GitHub: lecture_xspec](https://github.com/tenoto/lecture_xspec) に置いています。以下のようにダウンロードできます。GitHub が使える前提ですが、使えない方は zip file でダウンロードできるようにしています。
```
git clone https://github.com/tenoto/lecture_xspec.git
```
- [zip file](TBD)
コマンドラインで以下のように入力してみましょう。EOF は[ヒアドキュメント](https://qiita.com/take4s5i/items/e207cee4fb04385a9952)と呼ばれていて、実際に自分でコマンドプロンプトから入力する際には、EOF は無視して順番に入れてみてください。
```
xspec <<EOF
data 1 sample/sim_ninja_crab_100s_bin100.fak
none
resp 1 response/ninja_2gmc_whole_191224.rmf
back 1 sample/sim_ninja_crab_100s_bkg.fak
setplot energy
cpd /xw
plot ldata
ignore 1-1:**-1.0 50.0-**
```
これで表示されるのは、検出器で測定された「かに星雲」(Crab Nebula)のスペクトルです（注: ここでは実習のために、シミュレーションで作成したものです）。かに星雲は、Ｘ線天文衛星が打ち上げ後に真っ先に観測する標準光源としての意味のある天体ですね。
データは、 検出器で計数されたＸ線光子のカウント数をエネルギーに対してプロットしたもので、縦軸は Counts s$^{-1}$ keV$^{−1}$ で表示される。XSPEC の plot コマンドでは、誤差棒は Poisson 統計から予想される 1σ が default ではつくようになっています。
### モデル関数でのフィット
このとき、検出器に入射するスペクトル情報を引き出したいわけです。この観測データを再現するモデル関数として、天体が「べき」関数で表せる X 線スペ クトルをもち、星間吸収を受けて検出器まで届いたと考える。これは「べき」関数 (power) と星間吸収 (phabs) の積 phabs * pow で与えられるとし、検出器に入射したエネルギー E (keV) を もつ X 線光子の個数は、
$A(E)=\exp (-N_{\rm H}\sigma(E)) \times K\left(\frac{E}{1~\rm{keV}}\right)^{-\Gamma}$ photons cm$^{-2}$ s$^{-1}$ keV$^{-1}$
と表せる。星間吸収において、$N_{\rm H}$ atoms cm$^{-2}$ は水素柱密度、$\sigma (E)$ cm$^{-2}$ は光電吸収の断面積であり、「べき」関数において、K Photons cm−2 s−1 keV−1 は 1 keV で規格化した強度 (Normalization; Norm)、Γ は光子指数である。初期パラメータを入れてみる。
```
model phabs * power
0.4
2.0
9.0
fit
```
まず初期パラメータである程度、実験データに近いものを見つけたあとで、実験データを最もよく再現するモデルパラメータの最良推定値を探してみよう。
$\chi^2$ の最小値を求めるアルゴリズムには複数あり、XSPEC のデフォルトは Bevington の教科書で
の CURFIT ルーチンに基づいた Levenberg-Marquardt 法 (leven) を用いている。Levenberg-Marquardt 法は$\chi^2$の関数面の傾きから最小値を探索する最急降下法 (Steepest descent method)とガウス-ニュートン法 (Gauss-Newton method)を合わせたようなものである。このほかにCERN ライブラリの MINUIT 関数も用いることができ、
```
method simplex
```
などとして切り替えることができる。method としては、migrad (MINUIT の MIGRAD), simplex (SIMPLEX)が用意されている。$\chi^2$を最小化するようにモデル関数のパラメータを最適化しデータに合わせ込むには、
```
fit 200
```
とする。これで、200回の試行を繰り返しを行い (iteration)、$\chi^2$の最小値$\chi^2_{\rm min}$を見つけ出す。収束した場合、下の図のように、データによく合うモデル関数（赤）が得られ、中図のように、両者の残差は 0 を中心として分布する。
図: X線スペクトルデータ(黒点)とスペクトルのモデル関数(赤)。(top)検出器で計数された生のスペクトルと最適なモデル関数の解。(middle) モデルに対するデータの比。(data-model)/error が表示されている。(bottom) 検出器の応答関数を考慮して算出される、検出器に入射する前の段階での天体のスペクトルを$\nu F_{\nu}$で表示している。星間吸収は含まれている。
試行回数を限定せずにフィットするには、
```
query yes
fit
```
などとする。フィットが収束すると、以下のような最適なパラメータと統計量が表示される。
```
========================================================================
Model phabs<1>*powerlaw<2> Source No.: 1 Active/On
Model Model Component Parameter Unit Value
par comp
1 1 phabs nH 10^22 0.614496 +/- 0.127168
2 2 powerlaw PhoIndex 2.14982 +/- 3.11995E-02
3 2 powerlaw norm 10.9687 +/- 0.666350
________________________________________________________________________
Fit statistic : Chi-Squared 52.66 using 66 bins.
Test statistic : Chi-Squared 52.66 using 66 bins.
Null hypothesis probability of 8.20e-01 with 63 degrees of freedom
```
これを読むと、66 ビンのデータを 3 パラメータのモデルでフィットしたので、自由度 $\nu$ は 66-3=63 であり、得られた $\chi^2=52.66$ であるから、フィットの妥当性の目安として用いられる reduced chi-square の値は $\chi^2_{\nu}=0.84$ 程度となったことがわかる。
$\chi^2$ 分布において、上記の自由度で$\chi^2$が得られる確率は、0.82 (Null hypothesis probability) であるから、もし 仮に 100 回同一の観測を行えたとすれば、このモデルに合致する観測結果が 82 回得られる、つまり最適なパラメータをもつモデル関数はデータをよく再現すると言える。このようにして妥当なモデル関数を得ることができた。
一般に、$0.05 < \chi^2_{\nu} < 1.20$ 程度の範囲であればモデル関数が妥当であると言える場合が多いが、最終的には null hypothesis probability に基づいた、つまり帰無仮説を背景にした考察に戻るべきである。なお、$\chi^2$が自由度 $n$ より非常に小さい場合、つまり $\chi^2_{\nu}$ が 0 に近 い場合、フィットがよく合っていると考えるよりも、誤差評価に誤りがあるなどを疑うべ きであろう。また、このフィットにおいて保証されるのは「用いたモデル関数がデータを表すものとして適切であったか」であり、「そのモデル関数が物理的に適切かどうか」に は言及していないことに注意しておいた方がよい。フィットを合わせられることと、デー タを正しく物理解釈することは等価ではない。
なお、検出器の応答関数を考慮し、天体から地球に届いた段階での検出器に入る前の $\nu F_{\nu}$ スペクトルは、
```
plot eeuf
```
とすれば表示され、上図の最後のパネルのような形になる。この天体は非熱的な「べき」成分の天体スペクトルが星間吸収を受けてきたものと解釈できる。
### パラメータ誤差の推定
次のステップとして、モデル関数パラメータの誤差の範囲 (信頼できる値の範囲) を error コマンドを用いて求めよう。XSPEC の error コマンドでは、fit コマンドで求まった $\chi^2_{\rm min}$ に対して、与えられた $\triangle \chi^2$ を足した $\chi^2_{\rm bound}=\chi^2_{\rm min}+\triangle \chi^2$ 以下の$\chi^2$を与えるようなパラメータ値の範囲を計算する。この際、誤差を求めたいパラメータ以外の、値を固定されていないフリーな全てのパラメータをふり$\chi^2$の最小化を行う。これは自由度 1 の$\chi^2$分布に対応し、XSPEC のデフォルト設定では 90% の信頼区間に対応する $\triangle \chi^2 =2.73$ の範囲から誤差を求める。星間吸収、「べき」の光子指数、強度はそれぞれパラメータ番号が 1, 2, 3 で、あれば、以下のように90% 信頼区間に対応する誤差を求めることができます。
```
XSPEC12>err 1,2,3
Parameter Confidence Range (2.706)
1 0.40639 0.830787 (-0.208098,0.216298)
2 2.09934 2.20165 (-0.0504843,0.0518309)
3 9.93301 12.1332 (-1.03559,1.16464)
```
つまり、以下のように 90% 信頼区間の誤差を得ることができます。
- 星間吸収: $N_{\rm H}=(6.1\pm 0.2)\times 10^{11}$ cm$^{-2}$
- 「べき」の光子指数: $\Gamma = 2.10\pm 0.05$
- 1 keV での強度: normalization = $11.0_{-1.0}^{+1.2}~$ photons sec$^{-1}$ keV$^{-1}$ cm$^{-2}$
実は、このデータは、以下のパラメータをもつモデルから擬似的に作成したもので、以下が本来の値です。
- $N_{\rm H}=3.8\times 10^{11}$ cm$^{-2}$
- $\Gamma = 2.1$
- Normalization = $10$ photons sec$^{-1}$ keV$^{-1}$ cm$^{-2}$ at 1 keV
- ちなみに、Crab Nebula + Pulsar のスペクトルのパラメータは歴代のＸ線検出器でばらつきがある。[Enoto Ph.D thesis Figure A.1](https://jaxa.repo.nii.ac.jp/?action=repository_action_common_download&item_id=30269&item_no=1&attribute_id=31&file_no=1)など参照。
光子指数と強度は誤差の範囲で正しく得られています。星間吸収は低エネルギー側で決まるので、あまりうまくパラメータを得られていないようです。
ちなみに、90% 信頼区間の誤差はＸ線天文学でよく使われていたが、より一般的な$1\sigma$誤差 (68% 信頼度) の場合には、$\triangle \chi^2=1.0$ として、
```
XSPEC12>err 1.0 1,2,3
Parameter Confidence Range (1)
1 0.487016 0.744919 (-0.127472,0.130431)
2 2.11897 2.18116 (-0.0308556,0.0313376)
3 10.3255 11.6591 (-0.6431,0.690492)
```
とすればよい。
これを確認するために、フリーなパラメータのうち 2 つについて $\chi^2$値の 変化を描くことができる。今回、星間吸収はあまりうまく決まっていなかったので、$N_{\rm H}$の値を最適値に固定してみる。つまり、
```
freeze 1
```
その上で、
```
steppar 2 2.08 2.22 100 3 9.8 12.2 100
```
などとすると、2 番目のパラメータ (今回の場合、べきの光子指数) を 2.08 から 2.22 まで 100 ステップ、3 番目のパラメータ (べきの強度) を 9.8 から 12.2 までの 100 ステップにそれぞれ分割し、100 × 100 のグリッド点について他のパラメータを最適化した上で計算する。この結果は
```
plot contour,,5,1.0,2.30,2.70,4.61,9.21
```
とすると、以下のように表示できる。
この図では、$\chi^2$ = 1.00, 2.30, 2.70, 4.61, 9.21 をもつ 5 つの confidence contour を描いてくれる。2 個のパラメータ空間において残り 1 個のパラメータを最適化して $\chi^2$ を求めるので、この場合は自由度 2 の $\chi^2$ 分布に従い、3 つの confidence contour はそれぞれ信頼度 68% ($\chi^2=2.30$)、90% ($\chi^2=4.61$)、99% ($\chi^2=9.21$) に対応することになる。一方で、XSPEC の標準的な解析で error コマンドを使用する場合には、パラメータ 1 個 ずつについて誤差を計算していので自由度 1 の $\chi^2$分布であり、信頼度 68%、90%での 誤差は $\triangle \chi^2=1.0, 2.7$ で決まる。自由度 2 の 場合よりも誤差範囲はわずかに小さくなる。一般的な解析では error コマンドを使用することが多いが、どうのようなデータ解析を行っているのかには常に注意が必要である。また、実験結果をまとめる際に、$1\sigma$ 誤差をつけることが多いが、XSPEC では 90% 信頼度 での誤差がデフォルトのことが多い。いずれも定義して使えば問題ないが、注意しよう。
### カウントレートとＸ線フラックス
検出器のカウントレートは、ignore コマンドで指定したエネルギー範囲について
```
XSPEC12>show rate
Spectral Data File: sample/sim_ninja_crab_100s_bin100.fak
Assigned to Data Group 1 and Plot Group 1
Net count rate (cts/s) for Spectrum:1 1.018e+02 +/- 1.015e+00 (98.9 % total)
Spectral data counts: 10296
Model predicted rate: 101.304
```
とすると表示され、実測値は 101 counts s$^{-1}$ (1-50 keV)とわかる。モデルの推定値や、合計のカウント数も表示される。
フラックスを求めるときは、
```
XSPEC12>flux 2.0 10.0
Model Flux 3.2971 photons (2.1224e-08 ergs/cm^2/s) range (2.0000 - 10.000 keV)
XSPEC12>flux 2.0 10.0 err 200,68.3
Parameter distribution is derived from fit covariance matrix.
Model Flux 3.2971 photons (2.1224e-08 ergs/cm^2/s) range (2.0000 - 10.000 keV)
Error range 3.259 - 3.328 (2.097e-08 - 2.139e-08) (68.30% confidence)
```
などとする。```flux 2.0 10.0``` としたときは、単に 2-10 keV のフラックスが求まり、$2.1\times 10^{-8}$ ergs cm$^{-2}$ s$^{-1}$ と求まる。```flux 2.0 10.0 err 200, 68.3``` などとすると、誤差も求めてくれる。この場合は、1$\sigma$誤差。詳しくは、[ヘルプ](https://heasarc.gsfc.nasa.gov/xanadu/xspec/manual/node98.html)を調べるとよい。
ちなみに、よく使われる換算値は、1 Crab = $2.3\times 10^{-8}$ ergs s$^{-1}$ cm$^{-2}$ である。
Ｘ線光度は
$L_{\rm x}=4\pi d^2 F_{\rm x}=1.2\times 10^{32}(d/1~{\rm kpc})^2(F_{\rm x}/10^{-12}~{\rm ergs/sec/cm}^2)$ erg/s
なので、Crab まで 2 kpc の距離を考えると、このＸ線フラックスから推定されるＸ線光度は、$L_{\rm x}=1.1\times 10^{37}$ ergs/s となる。
余談だが、Crab が銀河中心 (d=8 kpc)にいると、ちょうど $L_{\rm x}=2\times 10^{38}$ ergs/s なので、中性子星のエディントン光度くらいになる。覚えておくと便利なときもある。
### 検出器の応答関数とフィットの関係
最後に XSPEC の中での検出器応答の扱いを簡単に補足しておく。 簡単な物理観測では、測定値に対して検出器の検出効率の補正を施せば、元々の物理量、観測量に簡単に戻すことが可能である。一方でＸ線観測においては、検出器に光子の全エネルギーを与える光電吸収のほかにも、一部分のエネルギーを与えるコンプトン散乱などの過程もある確率で生じている。そのため、検出器の実際の計数から元の Ｘ線スペクトルを一意的に求めることは、一般にはできない。そこで、X 線のモデルスペクトルを仮定し、検出器応答を受けた場合に予想されるカウント数のスペクトルを、実際の観測結果と比較する。そして$\chi^2$最小化のアルゴリズムによって、再帰的に最適なモデルパラメータを求めるという方法を採用する。これは、
$M(PI_j)=\sum_i S\cdot A(E_i)\cdot R(E_i,PI_j) M(E_i) \triangle E_i$
と表せる。ここで、
- $M(PI_j)$: 検出されるイベント数 (counts cm$^{-2}$ s$^{-1}$)
- $S$: 有効面積 (cm$^2$)
- $A(E_i)$: 検出効率
- $R(E_i,PI_j)$: 検出器の応答関数
- $M(E_i)$: 天体Ｘ線のモデル関数 (photons cm$^{-2}$ s$^{-1}$ keV$^{-1}$)
- $\triangle E_i$: エネルギー幅 (keV)
である。エネルギー$E_i$ keVの単色X線で予想されるスペクトルは、$A(E_i)R(E_i,PI_j)$ であり、エネルギー$E_i$ keVに対応する検出効率は$A(E_i)\sum j=R(E_i,PI_j)$と表せる。 一般に$R(E_i,PI_j)$に逆行列が必ずしも見つかるわけではないため、$\chi^2$最小化アルゴリズ ムを採用しているとも考えられる。
XSPEC では $R(E_i,PI_j)$ を RESP (Response File)、$S\cot A(E_i)$ を ARF (Ancillary Response File) として読み込むが、検出器ごとに S をどこが担うかは異なる場合がある。
たとえば、「すざく」の X 線 CCD カメラ (XIS) では、ARF、RESP を両方用いているが、硬 X 線検 出器の PIN 型半導体 (HXD-PIN) では、RESP のみに全てを含めるようになっている。また、硬 X 線検出器の GSO 結晶シンチレータ (HXD-GSO) では、若干の補正として ARF を用いるようにしている。それぞれがどのような関数形をしているかは、FTOOLS の fv コマンドなどで中身を見てみるとよい。
### スペクトルのビンまとめと C-stat
統計の悪いデータを解析する際、ビンまとめという操作を行うことがある。```grppha``` などのコマンドで行うことができる。
なお、昔はビンまとめすることが多かったが、これはエッジやラインなどのシャープな構造がもっている情報を消してしまうので、C-stat を使う人も増えている。
## シミュレーションスペクトルを作る
観測を行う前に、どのようなスペクトルを持つかなどをシミュレーションしたいことも多い。こういった場合、fakeit という xspec のコマンドを使うことができる。
まず、仮定として、0.3 keV の温度をもつ中性子星が 1 kpc の距離にある場合にどう見えるかを考えてみよう。つまり、
- 黒体放射の温度 kT = 0.3 keV
- 放射の半径 R = 12 km
- 距離 d = 1 kpc
- 上記の仮定から、[bbodyrad model](https://heasarc.gsfc.nasa.gov/xanadu/xspec/manual/node139.html) の normalization $= R_{\rm km}^2/D_{\rm 10kpc}^2 = 1.44\times 10^4$
- 過去の観測例から、$N_{\rm H}=5\times 10^{21}$ cm$^{-2}$とする。
まず、xspec でソース・スペクトル（なんでもよい）とバックグラウンド、レスポンスを読み込んだ後、モデル関数を指定する。そして、fakeit コマンドで観測時間を指定する。これで、だいたいの統計誤差の大きさや、どれくらいのエネルギー帯域で検出ができそうかなどを把握することができる。
```
xspec
data 1 response/ninja_empty.pi
back 1 response/ninja_2gmc_whole_nxbcxb_200928.pi
resp 1 response/ninja_2gmc_whole_191224.rmf
model tbabs * bbodyrad
0.5
0.3
14400
fakeit
y
sim_ninja_nstar_5ks.fak
5000
```
すると、sim_ninja_nstar_5ks.fak と sim_ninja_nstar_5ks_bkg.fak という２つのファイルができる。これを読み込んで、同じモデルでフィットしてみると、以下のような結果を得る。
```
========================================================================
Model TBabs<1>*bbodyrad<2> Source No.: 1 Active/On
Model Model Component Parameter Unit Value
par comp
1 1 TBabs nH 10^22 1.16972 +/- 0.393162
2 2 bbodyrad kT keV 0.279850 +/- 1.31500E-02
3 2 bbodyrad norm 3.24012E+04 +/- 1.69984E+04
________________________________________________________________________
Fit statistic : Chi-Squared 67.82 using 71 bins.
Test statistic : Chi-Squared 67.82 using 71 bins.
Null hypothesis probability of 4.83e-01 with 68 degrees of freedom
XSPEC12>err 1.0 1,2,3
Parameter Confidence Range (1)
1 0.789101 1.58597 (-0.380622,0.416248)
2 0.266361 0.293468 (-0.0134894,0.0136178)
3 19324.3 56813.2 (-13076.9,24411.9)
```
### 参考サイト
- [初めての宇宙X線スペクトル解析ツール xspecの使い方](https://qiita.com/yamadasuzaku/items/b978a5cac6197b6f1090) (山田真也さん＠立教大学)
- [XMM のデータ解析](https://qiita.com/yamadasuzaku/items/0c5f28438f7342b63975) (山田真也さん＠立教大学)
- [すざくファーストステップガイド](http://cosmic.riken.jp/suzaku/help/guide/fstep_web/node2.html)
- [Xspec11ユーザのためのXspec12入門講座](http://ytkyk.info/wiki/xspec/Xspec11%E3%83%A6%E3%83%BC%E3%82%B6%E3%81%AE%E3%81%9F%E3%82%81%E3%81%AEXspec12%E5%85%A5%E9%96%80%E8%AC%9B%E5%BA%A7.html) (HongoWiki)
Google Drive
Gist
Clipboard
Sign in with Wallet
