0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

PyXspec/Xspecで始める:fakeitによる疑似スペクトル生成と活用

Last updated at Posted at 2025-09-22

はじめに

X線天文学では、観測時に「どのくらいのシグナルが期待できるか?」を見積もりたい場面がしばしばあります。
その際に役立つのが fakeit というコマンドです。これは XspecPyXspec から利用でき、既知の応答ファイル(RMF/ARF)を使って任意のモデルから擬似観測スペクトルを生成できます。

fakeit でできることの例:

  • 観測提案で「この露光時間なら検出できそうか?」を検討する
  • 複数の衛星や検出器の性能を比較する(例:XRISM Resolve と Xtend)
  • 解析コードのテスト用データを作成する

本記事では、Python から XSPEC を操作する PyXspec を用いて、シンプルな fakeit の実行例と、検出器ごとの応答の違いを比較します。
あわせて、XSPEC コマンドラインでの実行例についても簡単に紹介します。


本記事で使用した RMF/ARF ファイルは、ZIP でまとめたものと展開済みの個別ファイルの両方を用意しています。
ZIP は一括ダウンロードに便利ですが、サイズがやや大きいため、必要な方は個別ファイルから取得することもできます。

PyXspec とは

PyXspec は、スペクトル解析ソフト Xspec を Python から操作できるようにした公式ラッパーです。
モデル設定、フィッティング、fakeit による擬似スペクトル生成など、Xspec のほとんどの機能を Python スクリプトから実行できます。

PyXspec の使い方については、以下の記事でデモンストレーションしていますので、もしよろしければあわせてご覧ください:

また、公式サイトにもチュートリアルが用意されています:

fakeit の最小サンプルコード

簡単な例として、6.4 keV に中心を持つ単一ガウシアン線(σ = 0.01 keV)を仮定したスペクトルを用いて、擬似観測を行う最小コードを示します。

import xspec as xs

# --- モデル定義(例:6.4 keV の単一ガウシアン線) ---
m = xs.Model("gaussian")
m.gaussian.LineE = 6.4
m.gaussian.Sigma = 0.01
m.gaussian.norm  = 1e-2

# --- Fakeit 設定(最小限) ---
fs = xs.FakeitSettings(
    response="xtd_source.rmf",       # 応答ファイル (RMF)
    arf="xtd_source.arf",         # 有効面積ファイル (ARF)
    exposure=1000,                   # 露光時間 [s]
    fileName="fake_xtd_gauss.fak"    # 出力ファイル名
)

# --- fakeit 実行 ---
xs.AllData.fakeit(nSpectra=1, settings=[fs], applyStats=True)

これで fake_xtd_gauss.fak という疑似スペクトルが生成されます。

引数の簡単な説明

  • response : RMF (Redistribution Matrix File)。検出器のエネルギー分解能を反映するファイル。

  • arf : ARF (Auxiliary Response File)。有効面積や感度を反映するファイル。

  • exposure : 露光時間 [秒]。長くするとカウント数が増える。

  • fileName : 出力される擬似スペクトルのファイル名。

  • applyStats: Poisson統計を適用するかどうかを選択

    • True → Poisson統計に従ってカウントがランダムに揺らぐ(観測に近い疑似データ)
    • False → Poisson統計の揺らぎを入れず、期待値そのまま(理想的なスペクトル=Poisson分布の平均値)

詳細仕様について

ここで紹介したのは代表的な引数ですが、このほかにもいくつかオプションがあります。
細かい仕様やすべての引数について知りたい方は、PyXspec公式ドキュメントをご覧いただければと思います。

可視化例

生成したスペクトルは matplotlib で簡単に描画できます。

import xspec as xs
import matplotlib.pyplot as plt
import numpy as np

plt.rcParams["font.size"] = 15

# --- 生成したスペクトルを読み込み ---
spec = xs.Spectrum("fake_xtd_gauss.fak")

# --- XSPEC 側で一度計算して配列取得 ---
xs.Plot.xAxis = "keV"
xs.Plot.device = "/null"
xs.Plot("data")

x = np.array(xs.Plot.x(1,1))
y = np.array(xs.Plot.y(1,1))
xerr = np.array(xs.Plot.xErr(1,1))
yerr = np.array(xs.Plot.yErr(1,1))
label = xs.Plot.labels(1)

# --- Matplotlib で描画 ---
plt.figure(figsize=(7,5))
plt.errorbar(x, y, xerr=xerr, yerr=yerr, fmt='o', ms=0, elinewidth=0.8)
plt.xlabel(label[0])
plt.ylabel(label[1])
plt.xlim(6, 6.8)
plt.show()

出力例
fakeit_xtd_gauss.png

入力したガウス線の真の広がりは σ = 0.01 keV ですが、観測されたスペクトルではそれよりもはるかに大きな幅になっています。
これは、検出器のエネルギー分解能の影響を強く受けて線が広がって見えていることを意味します。

輝線幅と観測の関係

観測される線幅は、ソース側の真の幅と検出器の分解能が二乗和で合成されます:

\sigma_{\text{observed}}^2 = \sigma_{\text{source}}^2 + \sigma_{\text{detector}}^2

一見すると「単純にソースと検出器の幅を重ね合わせただけ」と思いがちですが、実際にはこの性質が 直感とずれることが多い と筆者は感じています。

たとえば:

  • ソースの幅が検出器より小さい場合
    → 観測される幅は検出器に大きく引っ張られ、真の幅は見えにくくなります。
  • ソースの幅が検出器より大きい場合
    → 観測される幅はほとんどソースの幅で決まり、検出器の影響は目立たなくなります。

つまり「どんな場合でも検出器分解能が必ず効いている」ものの、寄与の仕方は相対的であり、思ったよりも検出器の影響が強く出たり、逆にほとんど消えてしまうこともあるのです。

このあたりは図で示すとより直感的に理解できるので、もしよろしければ こちらの記事 もご覧ください。

応用例

1. Poisson ノイズあり・なしの比較

通常、fakeit コマンドは Poisson ノイズを含む疑似スペクトルを生成します。
オプション applyStats=False を指定すると、統計的な揺らぎを入れない理想的なノイズなしスペクトルを得られます。

  • applyStats=True → 実際の観測に近い(ノイズあり)
  • applyStats=False → 応答関数を畳み込んだだけ(ノイズなし)

用途に応じて切り替えると便利です。

import xspec as xs
import matplotlib.pyplot as plt
import numpy as np

plt.rcParams["font.size"] = 15

# --- モデル定義 ---
m = xs.Model("gaussian")
m.gaussian.LineE = 6.4
m.gaussian.Sigma = 0.01
m.gaussian.norm  = 1e-2

# --- Fakeit 設定(2パターン) ---
fs_t = xs.FakeitSettings(
    response="xtd_source.rmf",
    arf="xtd_source.arf",
    exposure=1000,
    fileName="fake_xtd_true.fak"
)
fs_f = xs.FakeitSettings(
    response="xtd_source.rmf",
    arf="xtd_source.arf",
    exposure=1000,
    fileName="fake_xtd_false.fak"
)

# --- fakeit 実行 ---
xs.AllData.fakeit(nSpectra=1, settings=[fs_t], applyStats=True)
xs.AllData.clear()
xs.AllData.fakeit(nSpectra=1, settings=[fs_f], applyStats=False)

# --- 読み込み ---
xs.AllData.clear()
spec_t = xs.Spectrum("fake_xtd_true.fak")
spec_f = xs.Spectrum("fake_xtd_false.fak")

# --- XSPEC 側で配列を用意 ---
xs.Plot.xAxis = "keV"
xs.Plot.device = "/null"
xs.Plot("data")

xt, yt = np.array(xs.Plot.x(1,1)), np.array(xs.Plot.y(1,1))
xt_err, yt_err = np.array(xs.Plot.xErr(1,1)), np.array(xs.Plot.yErr(1,1))
xf, yf = np.array(xs.Plot.x(2,1)), np.array(xs.Plot.y(2,1))
xf_err, yf_err = np.array(xs.Plot.xErr(2,1)), np.array(xs.Plot.yErr(2,1))
label = xs.Plot.labels(1)

# --- プロット ---
plt.figure(figsize=(7,5))
plt.errorbar(xt, yt, xerr=xt_err, yerr=yt_err, fmt='o', ms=0, elinewidth=0.8, label="applyStats=True")
plt.errorbar(xf, yf, xerr=xf_err, yerr=yf_err, fmt='o', ms=0, elinewidth=0.8, label="applyStats=False")
plt.xlabel(label[0])
plt.ylabel(label[1])
plt.xlim(6, 6.8)
plt.legend()
plt.show()

出力例
fakeit_xtd_gauss_apply_true_false.png

図を見ると、applyStats=True では Poisson ノイズによってカウントがランダムに揺らいでいるのに対し、
applyStats=False では統計誤差のない滑らかなスペクトルになっていることがわかります。

2. 検出器(Resolve / Xtend / Chandra)の応答比較

fakeit が役立つ場面のひとつに、検出器の違いを比較できるという点があります。
単にモデルの分解能だけでなく、感度や有効面積なども込み込みでシミュレーションできるので、観測提案を考えるときにもとても重宝します。

ここでは XRISM に搭載されている 2 つの検出器──マイクロカロリーメータを搭載した Resolve と X 線 CCD を搭載した Xtend──および Chandra の HETG というグレーティング検出器を例に、同じモデルスペクトルを通したときの違いを比較してみます。

XRISM の RMF/ARF の詳細な特性については、以下の記事で紹介していますので、よろしければ合わせてご覧ください:

それでは、それぞれの検出器を通した疑似スペクトルを生成し、応答の違いを見てみましょう。

import xspec as xs
import matplotlib.pyplot as plt
import numpy as np

plt.rcParams["font.size"] = 15

# --- モデル定義 ---
m = xs.Model("gaussian")
m.gaussian.LineE = 6.4
m.gaussian.Sigma = 0.01
m.gaussian.norm  = 1e-2

# --- fakeit 設定(XRISM Resolve Xtend, Chandra HETG HEG±1) ---
fs_rsl = xs.FakeitSettings(
    response="rsl_source_Hp.rmf",
    arf="rsl_source_Hp.arf",
    exposure=1000,
    fileName="fake_xrism_rsl.fak"
)
fs_xtd = xs.FakeitSettings(
    response="xtd_source.rmf",
    arf="xtd_source.arf",
    exposure=1000,
    fileName="fake_xrism_xtd.fak"
)
fs_heg = xs.FakeitSettings(
    response="chandra_hetg_heg_abs1.rmf",
    arf="chandra_hetg_heg_abs1.arf",
    exposure=1000,
    fileName="fake_chandra_hetg_heg_abs1.fak"
)

# --- fakeit 実行 ---
xs.AllData.fakeit(nSpectra=3, settings=[fs_rsl, fs_xtd, fs_heg], applyStats=True)

# --- 読み込み ---
xs.AllData.clear()
spec_rsl = xs.Spectrum("fake_xrism_rsl.fak")
spec_xtd = xs.Spectrum("fake_xrism_xtd.fak")
spec_heg = xs.Spectrum("fake_chandra_hetg_heg_abs1.fak")

# --- XSPEC 側で配列を用意 ---
xs.Plot.xAxis = "keV"
xs.Plot.device = "/null"
xs.Plot("data")

xr, yr = np.array(xs.Plot.x(1,1)), np.array(xs.Plot.y(1,1))
xr_err, yr_err = np.array(xs.Plot.xErr(1,1)), np.array(xs.Plot.yErr(1,1))
xx, yx = np.array(xs.Plot.x(2,1)), np.array(xs.Plot.y(2,1))
xx_err, yx_err = np.array(xs.Plot.xErr(2,1)), np.array(xs.Plot.yErr(2,1))
xh, yh = np.array(xs.Plot.x(3,1)), np.array(xs.Plot.y(3,1))
xh_err, yh_err = np.array(xs.Plot.xErr(3,1)), np.array(xs.Plot.yErr(3,1))
label = xs.Plot.labels(1)

# --- プロット ---
plt.figure(figsize=(7,5))
plt.errorbar(xr, yr, xerr=xr_err, yerr=yr_err, fmt='o', ms=0, elinewidth=0.8, label="XRISM/Resolve")
plt.errorbar(xx, yx, xerr=xx_err, yerr=yx_err, fmt='o', ms=0, elinewidth=0.8, label="XRISM/Xtend")
plt.errorbar(xh, yh, xerr=xh_err, yerr=yh_err, fmt='o', ms=0, elinewidth=0.8, label="Chandra/HETG HEG±1")
plt.xlabel(label[0])
plt.ylabel(label[1])
plt.xlim(6, 6.8)
plt.legend(title="Detector", fontsize="small")
plt.show()

出力例
fakeit_xrism_vs_Chandra.png

  • 分解能の比較
    結果を見ると、Resolve は分解能が数 eV 程度と非常に高く、輝線をシャープに捉えているように見えます。
    一方、Xtend は分解能が百 eV 程度なので線が広がって見えます。
    Chandra/HETG HEG±1 はその中間で、Resolve よりは広がっていますが、Xtend よりはシャープに見える印象です。

  • 感度の比較
    感度については、ほぼ有効面積(ARF)の特性を反映していると考えてよさそうです。傾向としては Xtend > Resolve > Chandra の順に高いように見えます。特に、Chandra はグレーティングを使っているため有効面積が Resolve の一桁ほど小さく、光子を集めにくい構造になっています。もちろん Chandra/ACIS(グレーティングなし) では状況が異なり、ACIS は CCD のため分解能は劣るものの感度はそこそこ確保される、といった違いが出てきます。

  • bin の比較
    そして、最小 bin の違いも興味深い点だと思います。bin は細かいものをまとめることはできても、粗いものを後から細かく分けるのは難しく、不可逆的な操作になります。例えば Resolve は最小 bin が細かいため、統計誤差がやや大きく見える傾向がありますが、bin をまとめることで調整可能です。こうした最小 bin の違いも、検出器の性能の一つと言えると思います。
    XSPEC での bin まとめ方法については、山田さんの「X線解析ツール xspec 向けにスペクトルのビンまとめをする3つの方法」をご覧ください。

このように fakeit を使うと、検出器ごとの特性を簡単に比較できて面白いので、性能を理解する際にぜひ活用してみてください。

検出器の性能は時々刻々変化する

検出器の性能は、経年劣化や観測条件など様々な要因によって変わります。 例えば、どこまでをソース領域とするかによっても ARF が変化します。そのため、ここで紹介している結果はあくまで標準的なサンプルの一例として捉えてください。

fakeit で複数スペクトルを生成する際の注意点

AllData.fakeit() では、一度に複数のスペクトルを生成することができます。

まとめて作る場合の例

xs.AllData.fakeit(nSpectra=3, settings=[fs_rsl, fs_xtd, fs_heg], applyStats=True)

ただし applyStats は全体で一つの設定しか選べず、True / False を混在させることはできません
そのため、異なる設定で生成したい場合は個々に呼び出す必要があります。

個々に分けて書く場合、先ほどの「まとめて作る例」を展開すると以下のようになります。

個々に作る場合の例

xs.AllData.fakeit(nSpectra=1, settings=[fs_rsl], applyStats=True)
xs.AllData.clear()
xs.AllData.fakeit(nSpectra=1, settings=[fs_xtd], applyStats=True)
xs.AllData.clear()
xs.AllData.fakeit(nSpectra=1, settings=[fs_heg], applyStats=True)

ここでの注意点は xs.AllData.clear() を必ず挟むことです。
これを入れないと前に読み込んだ RMF/ARF がそのまま引き継がれてしまいます。
同じレスポンスを使う場合は一見問題なく動きますが、異なるレスポンスを指定したつもりでも前の設定で走ってしまうので注意してください。

3. Xspecを使用した例

これまで PyXspec の方法を紹介しましたが、コマンドラインの Xspec でも同様に fakeit が実行できます。
例えば、単一ガウス線モデルを Xspec で再現すると以下のようになります。

XSPEC12>model gaussian

Input parameter value, delta, min, bot, top, and max values for ...
            6.5       0.05(     0.065)          0          0      1e+06      1e+06
1:gaussian:LineE>6.4
            0.1       0.05(     0.001)          0          0         10         20
2:gaussian:Sigma>0.01
              1       0.01(      0.01)          0          0      1e+20      1e+24
3:gaussian:norm>1e-2

========================================================================
Model gaussian<1> Source No.: 1   Active/Off
Model Model Component  Parameter  Unit     Value
 par  comp
   1    1   gaussian   LineE      keV      6.40000      +/-  0.0          
   2    1   gaussian   Sigma      keV      1.00000E-02  +/-  0.0          
   3    1   gaussian   norm                1.00000E-02  +/-  0.0          
________________________________________________________________________

XSPEC12>fakeit none
For fake spectrum #1 response file is needed: xtd_source.rmf
   ...and ancillary file: xtd_source.arf
 Use counting statistics in creating fake data? (y): y
 Input optional fake file prefix: 
 Fake data file name (xtd_source.fak): fake_xtd_gauss_xspec.fak
 Exposure time, correction norm, bkg exposure time (1.00000, 1.00000, 1.00000): 1000, 1, 1

ここで fakeit none とすることで、実データをロードしていない状態からでも RMF/ARF を指定して擬似スペクトルを生成できます。

最後に尋ねられる correction normbkg exposure time は、特別な事情がなければデフォルト値(1.0)で問題ないと思います
(筆者自身はこれらを変更して使ったことがなく、説明に自信がありません。恐縮ですが、詳しい方にご相談いただければと思います。)

まとめ

今回は、PyXspec の fakeit コマンドを使って擬似スペクトルを生成し、その可視化や応用例を紹介しました。
さらに、Xspec コマンドラインからの fakeit 実行方法についても簡単に触れました。

  • applyStats=True/False を切り替えて観測らしさや理想スペクトルを再現できる
  • 異なる検出器応答を通したシミュレーションで性能を比較できる
  • Xspec コマンドラインでも同様の手順で再現可能

といった点が特に便利です。

この記事がみなさんの研究や解析に少しでも役立てば嬉しいです。

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?