はじめに

X線天文学では、観測装置の特性を理解し、正確なデータ解析を行うことが重要です。特に、X線スペクトルの解析では、観測データと物理モデルを比較するための応答関数（RMF・ARF） が不可欠であり、XSPECなどの多くの解析ソフトウェアで使用されます。

本記事では、XRISM X線天文衛星の観測装置（Resolve / Xtend）を例に、Python を用いた RMF・ARF の可視化手法を紹介します。実際のデータ解析に役立つよう、可視化を通じて各装置の特性や違いを直感的に理解することを目指します。

本記事のARF・RMF・ARFの作成は下記の公式ページに基づいて作成したものを用いています。

RMF・ARF・観測されるスペクトルの説明

RMF（Response Matrix File）

RMFは、X線望遠鏡の観測において、入射したX線のエネルギー $E_j$ が、どの観測エネルギー $E_i$ として記録されるかを表す変換行列です。これは、観測装置のエネルギー応答を記述したものであり、各入射エネルギーに対して観測エネルギーがどのように分布するかを2次元の行列（$R_{ij}$）として持ちます。そのため、RMFは「入射エネルギー × 観測エネルギー」の2次元データを持つ情報として扱われます。

R_{ij} = P(E_i | E_j)

ここで、ある特定の入射エネルギー $E_j$ に対して、X 線が検出器によってどのエネルギー $E_i$ に分布するかを記述するため、次の正規化条件を満たします。

\sum_{i} R_{ij} = 1

これは、ある入射エネルギー $E_j$ の X 線が検出された場合、その測定値がどこかのエネルギービン $E_i$ に割り当てられることを意味しています。

ARF（Auxiliary Response File）

ARFは、X線望遠鏡と検出器のシステム全体の有効面積（effective area）をエネルギーごとに表したものです。これは、各エネルギーにおける感度を示し、望遠鏡の集光効率や検出器の感度などを考慮しています。ARFは「入射エネルギー $E_{j}$ ごとの有効面積」という1次元の情報（$A_{j}$）を持ち、入射エネルギーに対する関数として表されます。

観測されるスペクトル

実際に観測されるスペクトル $D_i$ は、天体の真のスペクトル $S_j$ に対して、有効面積 $A_j$（ARF）と応答行列 $R_{ij}$（RMF）を適用した総和として表されます。

D_i = \sum_j R_{ij} A_j S_j

ここで：

$D_i$ ：観測されたエネルギー $E_i$ におけるスペクトル（観測データ）
$S_j$ ：真の入射エネルギー $E_j$ における天体のスペクトル（物理モデル）
$A_j$ ：有効面積（ARF）
$R_{ij}$ ：入射エネルギー $E_j$ の X 線が、観測エネルギー $E_i$ にどう分布するかを示す応答行列（RMF）

これは、「入射エネルギーごとの天体のスペクトル $S_j$」が、「望遠鏡と検出器の物理的応答（ARFとRMF）」によって変換された結果、観測されるスペクトル $D_i$ が得られることを意味しています。
数学的には $D_i$ は実数の分布を持ちますが、実際の観測では観測時間などの影響を受け、ポアソン統計に従うカウント値として記録されます。

RMFの可視化

RMFデータの特性

RMFデータは、その多くの成分がゼロであるため、スパース（疎行列） として効率的に管理されています（非ゼロ成分のみが格納されています）。しかし、可視化する際には、これらのスパース行列を2次元配列形式のフルマトリックスに戻す必要があります。
もちろん、スパース形式のまま可視化する方法も存在します（例えば、plt.scatter() を用いる方法など）が、ここでは取り扱いません。

このようなスパース行列を可視化する際、まずRMFデータがどのような構造になっているかを把握することが重要です。本節では、実データと似せて用意したトイデータを使い、スパース行列の取り扱い方法について説明します。

RMFファイルの構造の確認

スパース行列といってもイメージがつきにくいと思いますので、まずは実データを確認します。fstruct は、HEASoft（High Energy Astrophysics Software）のパッケージの一部として提供されているツールで、FITS形式のファイルを解析し、その構造を詳細に表示するための便利なコマンドです。

例えば、XtendのRMFのデータ構造を確認する際、以下のような出力が得られます。

XtendのRMF

$ fstruct xtd.rmf  
  No. Type     EXTNAME      BITPIX Dimensions(columns)      PCOUNT  GCOUNT
 
   0  PRIMARY                 16     0                           0    1
   1  BINTABLE MATRIX          8     34(6) 6400           31909344    1
 
      Column Name                Format     Dims       Units     TLMIN  TLMAX
      1 ENERG_LO                   E                   keV
      2 ENERG_HI                   E                   keV
      3 N_GRP                      I
      4 F_CHAN                     PJ(1)                             0     4095
      5 N_CHAN                     PJ(1)
      6 MATRIX                     PE(4096)
 
   2  BINTABLE EBOUNDS         8     12(3) 4096                  0    1
 
      Column Name                Format     Dims       Units     TLMIN  TLMAX
      1 CHANNEL                    J                                 0     4095
      2 E_MIN                      E                   keV
      3 E_MAX                      E                   keV

ResolveのRMF

ResolveのRMF

$ fstruct rsl_Hp_L.rmf

  No. Type     EXTNAME      BITPIX Dimensions(columns)      PCOUNT  GCOUNT
 
   0  PRIMARY                 16     0                           0    1
   1  BINTABLE MATRIX          8     34(6) 60000         368716996    1
 
      Column Name                Format     Dims       Units     TLMIN  TLMAX
      1 ENERG_LO                   E                   keV
      2 ENERG_HI                   E                   keV
      3 N_GRP                      I
      4 F_CHAN                     PJ(41)                            0    59999
      5 N_CHAN                     PJ(41)
      6 MATRIX                     PE(2199)
 
   2  BINTABLE EBOUNDS         8     12(3) 60000                 0    1
 
      Column Name                Format     Dims       Units     TLMIN  TLMAX
      1 CHANNEL                    J                                 0    59999
      2 E_MIN                      E                   keV
      3 E_MAX                      E                   keV

RMFにはさまざまな情報が含まれていますが、ここではその中から主要な要素について説明します。

チャンネル番号に対応するエネルギー配列
- ENERG_LO: 入力エネルギーのビンごとの下限値（keV単位）
- ENERG_HI: 入力エネルギーのビンごとの上限値（keV単位）
- E_MIN: 出力エネルギーのビンごとの下限値（keV単位）
- E_MAX: 出力エネルギーのビンごとの上限値（keV単位）
RMFの読み込み用の情報
- F_CHAN: 基準となるチャネネル番号（出力チャンネルの開始位置）
- N_CHAN: F_CHANから何チャンネル分が含まれているか（出力チャンネル範囲の幅）
- MATRIX: スパース行列（観測器の応答を格納した非ゼロ部分のみを抽出したデータ）
  - スパース形式からフルサイズの行列を再構築するには、F_CHAN（開始位置）と N_CHAN（幅）を使って、非ゼロデータ（MATRIX の一部）を指定された範囲に埋め込む処理を行います。

トイデータでスパース行列の可視化

RMFのスパース行列を言葉で説明するだけではイメージしにくい場合もあるため、まずはトイデータを用いて可視化します。

トイデータでスパース行列の可視化

import numpy as np
import matplotlib.pyplot as plt

# サンプルデータを縮小
num_in_ene_bins = 3  # 出力チャンネルの数
num_out_ene_bins = 10  # 出力チャンネルの数

# 簡略化したスパースデータ
f_chan = [[2, 6], [1, 5], [3]]  # 開始チャンネル
n_chan = [[2, 3], [3, 2], [4]]  # チャンネル数
compressed_matrix = [
    [0.2, 0.3, 0.2, 0.1, 0.2],  # チャンネル0（出力の合計が1）
    [0.5, 0.1, 0.1, 0.2, 0.1],  # チャンネル1（出力の合計が1）
    [0.1, 0.3, 0.4, 0.2],       # チャンネル2（出力の合計が1）
]

# フルマトリックスを初期化
matrix = np.zeros((num_in_ene_bins, num_out_ene_bins))

# スパースデータをフルマトリックスに展開
annotations = []
for j in range(len(f_chan)):
    tot_n = 0
    for f, n in zip(f_chan[j], n_chan[j]):
        matrix[j, f:f + n] = compressed_matrix[j][tot_n:tot_n + n]
        annotations.append((j, f, f"F={f}\n       → N={n}"))
        tot_n += n

# 描画
fig, ax = plt.subplots(figsize=(8, 3))
cax = ax.imshow(matrix, cmap="viridis", aspect="auto", origin='lower')
plt.colorbar(cax)

# 軸ラベルを設定
ax.set_xlabel("Output Channels")
ax.set_ylabel("Input Channels")
ax.set_title("Sparse to Full Response Matrix")
plt.xticks(range(num_out_ene_bins))
plt.yticks(range(num_in_ene_bins))

# f_chanの値を赤字で表示
for j, f, text in annotations:
    ax.text(f, j, text, color="red", ha="center", va="center", fontsize=15)

fig.tight_layout()
plt.savefig('sparse_to_full_rmf.png', bbox_inches='tight')
plt.show()

出力結果

F: F_CHAN（開始位置）
N: N_CHAN（チャンネル数）

上図では、スパース行列がどのようにフルマトリックスとして展開されるかを視覚的に示しています。

コードの説明

j（ループのインデックス）
- 入力チャンネルを指します（例: 0番目のチャンネル、1番目のチャンネルなど）。
- f_chan[j] と n_chan[j] は、このビンに対応する出力チャンネル情報を持っています。
f_chan[j]
- この入力チャンネルにおける出力チャンネルの開始位置を示します。
- 例: f_chan[j] = [2, 6] なら、このビンの非ゼロ応答データはチャンネル2から開始し、次にチャンネル6から開始する応答が続いています。
n_chan[j]
- f_chan[j] の各開始位置から、何出力チャンネル分の応答を含んでいるかを示しています。
- 例: f_chan[j] = [2, 6]で、n_chan[j] = [2, 3] なら、出力チャンネル2から始まり2チャンネル分、チャンネル6から始まり、3チャンネル分の応答を含むことを意味します。
compressed_matrix[j]
- スパース形式で格納された応答データ（非ゼロ部分のみ）。
- tot_n を利用して、現在どの部分のスパースデータを展開すべきかを管理しています。

RMFスパース行列の要点整理

compressed_matrix は、スパース行列の中で実際に値が存在する（非ゼロの）部分だけが詰め込まれた配列です。
f_chan は開始位置のチャンネル番号を示し、n_chan でチャンネル数を表します。
この両者を組み合わせることで、どの位置に非ゼロ値が割り当てられるか を特定し、それ以外をゼロ埋めすることで完全なマトリックスを再構築することが可能になります。

チャンネル情報とエネルギー情報を組み合わせることで、RMFの完全なレスポンスマトリックスが完成します。次節では、実際にXRISM衛星のRMFを用いて可視化していきます。

RMFの読み込みと可視化

XtendとResolveのRMFを astropy.io.fits を用いて読み込み、スパース形式のレスポンスマトリックスを完全なフルマトリックスに再構築します。さらに、これをエネルギー軸に対応させて可視化することで、両検出器のエネルギー応答特性を比較します。

なお、XtendとResolveではチャンネル数が異なるため、事前にそれぞれの構造を以下のように整理します：

Xtend
- 入力チャンネル数: 6,400
- 出力チャンネル数: 4,096
Resolve
- 入力チャンネル数: 60,000
- 出力チャンネル数: 60,000

Resolve の RMF には、データ量と解析の目的に応じた 4種類 があります。詳細は「補足 Resolve RMFの種類」で説明しますが、ここでは rslmkrmf の RMF 生成コマンドで作成される whichrmf=L を使用します。

Resolveは非常に多くのチャンネル数を持つため、可視化の効率性を向上させるために、出力チャンネルをビンまとめして可視化します。このビンまとめ処理では、指定された数のチャンネルを単純に足し込むことで、元のレスポンス確率が1を保たれるように設計されています。

このビンまとめ処理は、次に示す binning_out_ene() 関数を利用してbin_size引数で指定できるようになっています。

RMFの読み込みと可視化

import os
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from matplotlib.ticker import AutoLocator, AutoMinorLocator
import astropy.io.fits as pf

# matplotlibのグローバル設定（プロットの見た目を統一）
plt.rcParams['xtick.direction'] = 'in'  # x軸の目盛りを内向きに設定
plt.rcParams['ytick.direction'] = 'in'  # y軸の目盛りを内向きに設定


def construct_rmf_matrix(filename):
    """
    RMFからレスポンスマトリックスを構築する関数。

    Parameters:
        filename (str): FITSファイルのパス。

    Returns:
        dict: 構築されたレスポンスマトリックスと関連データを含む辞書。
    """
    with pf.open(filename) as f:
        # 必要なデータを抽出
        f_chan = f[1].data['F_CHAN']  # チャンネルの開始位置
        n_chan = f[1].data['N_CHAN']  # チャンネルの数
        compressed_matrix = f[1].data['MATRIX']  # 密な（非ゼロのみの）マトリックスデータ
        in_ene_lo = f[1].data['ENERG_LO']  # 入力エネルギー下限
        in_ene_hi = f[1].data['ENERG_HI']  # 入力エネルギー上限
        out_ene_lo = f[2].data['E_MIN']  # 出力エネルギー下限
        out_ene_hi = f[2].data['E_MAX']  # 出力エネルギー上限

        # フルレスポンスマトリックスの次元を決定
        num_in_ene_bins = len(in_ene_lo)  # 入力エネルギーのビン数（行数）
        num_out_ene_bins = len(out_ene_lo)  # 出力エネルギーのビン数（列数）

        # フルレスポンスマトリックスを構築
        matrix = np.zeros(shape=(len(in_ene_lo), len(out_ene_lo)), dtype='float32')
        for j in range(f_chan.shape[0]):
            tot_n = 0
            for f, n in zip(f_chan[j], n_chan[j]):
                # 密なデータをフルマトリックスに埋め込む
                matrix[j, f:f + n] = compressed_matrix[j][tot_n:tot_n + n]
                tot_n += n

    return matrix, in_ene_lo, in_ene_hi, out_ene_lo, out_ene_hi

def binning_out_ene(matrix, bin_size, ene_lo, ene_hi):
    """
    2D配列のx方向（列方向）でビニングを実行し、出力エネルギー境界もビニングする関数。

    Parameters:
        matrix (numpy.ndarray): 2D入力配列。
        bin_size (int): 列方向でのビンサイズ（足し込む列数）。
        ene_lo (numpy.ndarray): 出力エネルギーの下限配列。
        ene_lo (numpy.ndarray): 出力エネルギーの上限配列。

    Returns:
        tuple:
            numpy.ndarray: ビニングされた2D配列。
            numpy.ndarray: ビニングされた出力エネルギー境界。
    """
    # 必要な列数を計算し、ゼロ埋めで配列を拡張
    num_cols = matrix.shape[1]
    padding = (bin_size - (num_cols % bin_size)) % bin_size  # 不足分の列数を計算
    padded_matrix = np.pad(matrix, ((0, 0), (0, padding)), mode='constant', constant_values=0)

    # 配列をリシェイプしてビンごとに合計
    new_shape = (matrix.shape[0], padded_matrix.shape[1] // bin_size, bin_size)
    binned_matrix = padded_matrix.reshape(new_shape).sum(axis=2)

    # 出力エネルギー境界をビンサイズでサンプリング
    binned_out_ene_edges = np.append(ene_lo[::bin_size], ene_hi[-1])

    return binned_matrix, binned_out_ene_edges


bin_size = 200  # resolveの場合はデータ量の関係で bin_size = 200 程度を推奨
rmf_filename = 'rsl_Hp_L.rmf'
out_filename = f'resolve_rmf_matrix_bin_{bin_size}_X.png'
# rmf_filename = 'xtd_source.rmf'
# out_filename = f'xtend_rmf_matrix_bin_{bin_size}.png'

# RMFデータの読み込みと構築
matrix, in_ene_lo, in_ene_hi, out_ene_lo, out_ene_hi = construct_rmf_matrix(rmf_filename)

# エネルギー境界を計算
in_ene_edges = np.append(in_ene_lo, in_ene_hi[-1])

# 出力エネルギー方向をビニング
binned_matrix, binned_out_ene_edges = binning_out_ene(matrix, bin_size, out_ene_lo, out_ene_hi)

# プロット作成
fig, ax = plt.subplots(figsize=(8, 6))
c = ax.pcolormesh(
    binned_out_ene_edges,
    in_ene_edges,
    binned_matrix,
    cmap='turbo',
    norm=LogNorm()
)

ax.set_xlabel('Output (keV)')
ax.set_ylabel('Input (keV)')
ax.set_aspect('equal')
fig.colorbar(c, ax=ax, label='Probability')

# 軸の設定とグリッド
ax.xaxis.set_major_locator(AutoLocator())
ax.xaxis.set_minor_locator(AutoMinorLocator())
ax.yaxis.set_major_locator(AutoLocator())
ax.yaxis.set_minor_locator(AutoMinorLocator())
ax.grid(which='major', color='gray', linestyle='--', linewidth=0.5)
ax.grid(which='minor', color='lightgray', linestyle=':', linewidth=0.3)

# ファイル名をタイトルに設定し保存
base_filename = os.path.splitext(out_filename)[0]
ax.set_title(base_filename)
fig.savefig(out_filename, bbox_inches='tight')
plt.show()

出力結果
以下のプロットは、縦軸に入力エネルギー、横軸に出力エネルギーを取り、色は確率を示しています。それぞれの検出器（XtendとResolve）のレスポンス特性が視覚的に比較できます。
Xtend（出力：bin_size=1）

Resolve（RMF: rslmkrmf (whichrmf=L), 出力: bin_size=200）
（binsizeが小さいと描画に時間がかかるので注意）

RMFマトリックスの説明

対角線の強度: 入射エネルギーと出力エネルギーが一致する対角線上には強いピークが見られます。これらのピークは、検出器が入力エネルギーを正確に検出した際のスペクトルのエネルギー分解能を反映しており、LSF（Line spread function¹）を示しています。特に、Xtend に比べて Resolve はエネルギー分解能が桁で高く、よりシャープなピークを持つことがわかります。LSFは、は、検出器の分解能や材質の特性に依存します。特に、エネルギーが高くなるにつれて検出効率が低下し、分布が広がる傾向が確認できます。
エスケープピーク: Xtend・Resolve の両方で、対角成分よりも少し低いエネルギー帯に斜め方向の分布が見られます。これは、Si 検出器特有のエスケープピークによるものです²。さらに、Resolve では、検出器の HgのM線やTeのL線などに起因するエスケープピークも確認できます³。
エネルギー分解能: 分解能が低い場合、対角線から外れる成分が広がる傾向があります。Resolveのような高分解能検出器ではピークがシャープに集中しますが、Xtendでは分解能が低いため、広がったピークが現れます。この違いにより、両者のエネルギー応答特性を可視化で比較することが可能です。
非ゼロ要素の分布: 非ゼロ要素の分布は、検出器が観測可能なエネルギー範囲を示します。Siを用いた検出器のため、1.7 keV付近のSiのピークも顕著に見られます。

単色エネルギーでXtend・ResolveのRMFの比較

RMFの違いを比較するため、入力エネルギー 6.4 keV（Fe Iに相当）に対する出力エネルギー分布をプロットします。

1.7 keV付近のSiのピークと4.7 keVのSiのエスケープピークが見えます。そのほか、Resolveは、Hg, Te由来のエスケープピークも見えているのがわかります。

RMFのエネルギー依存性の確認

一つ前の節では単色のレスポンスを比較しましたが、ここではエネルギー依存性に着目します。Xtend と Resolve の RMF を 2D ヒストグラムで可視化することで、全体の分布を把握できます。

Xtend と Resolve はエネルギー分解能が 1桁違い、さらにエスケープ成分の寄与も異なるため、それぞれの特性が適切に見えるようにグラフの表示方法を調整しています。

可視化のポイント

Xtend → log スケールを使用
Resolve → symlog スケールで、小さい成分（$10^{-5}$以下）はリニア、大きい成分は対数で表示
カラースケール → 入力エネルギーごとに色分けし、分布の違いを視覚化

RMFエネルギー依存性の可視化コード

Xtend: RMFのエネルギー依存性の確認

import os
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from matplotlib.colors import LogNorm
from matplotlib.ticker import AutoLocator, AutoMinorLocator
import astropy.io.fits as pf

# matplotlibのグローバル設定（プロットの見た目を統一）
plt.rcParams['xtick.direction'] = 'in'  # x軸の目盛りを内向きに設定
plt.rcParams['ytick.direction'] = 'in'  # y軸の目盛りを内向きに設定
plt.rcParams['font.size'] = 16

# ~~~中略(前節のコード)~~~

rmf_filename = 'xtd_source.rmf'
out_filename = 'xtend_rmf_specific_ene_output.png'

# RMFデータの読み込みと構築
matrix, in_ene_lo, in_ene_hi, out_ene_lo, out_ene_hi = construct_rmf_matrix(rmf_filename)

# 各チャンネルの平均エネルギー値を計算
in_ene_mean = (in_ene_lo + in_ene_hi) / 2
out_ene_mean = (out_ene_lo + out_ene_hi) / 2

# プロット作成
fig, ax = plt.subplots(figsize=(18, 9))

# カラーマップの設定
vmin, vmax = 1, 20
cmap = plt.get_cmap("plasma")
norm = plt.Normalize(vmin=vmin, vmax=vmax)

for target_ene in range(vmin, vmax+1, 1):
    closest_index = np.abs(in_ene_mean - target_ene).argmin()
    closest_ene = in_ene_mean[closest_index]
    color = cmap(norm(closest_ene))
    ax.plot(out_ene_lo, matrix[closest_index], label=f'{closest_ene:.0f}', color=color)

# Xtendは普通にlogでも見やすい
ax.set_yscale('log')

ax.set_xlim(-0.5, 21.5)

ax.xaxis.set_major_locator(AutoLocator())
ax.xaxis.set_minor_locator(AutoMinorLocator())

ax.grid(which='major', color='gray', linestyle='--', linewidth=0.5)
ax.grid(which='minor', color='lightgray', linestyle=':', linewidth=0.3)
ax.set_xlabel('Output (keV)')
ax.set_ylabel('Probability')
fig.legend(title='Input (keV)', loc='upper right', bbox_to_anchor=(0.99, 0.91))

# ファイル名をタイトルに設定し保存
ax.set_title('Xtend: Energy Dependence of RMF')
fig.savefig(out_filename, bbox_inches='tight')

plt.show()

Resolve: RMFのエネルギー依存性の確認

import os
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from matplotlib.colors import LogNorm
from matplotlib.ticker import SymmetricalLogLocator, LogLocator, AutoLocator, AutoMinorLocator
import astropy.io.fits as pf

# matplotlibのグローバル設定（プロットの見た目を統一）
plt.rcParams['xtick.direction'] = 'in'  # x軸の目盛りを内向きに設定
plt.rcParams['ytick.direction'] = 'in'  # y軸の目盛りを内向きに設定
plt.rcParams['font.size'] = 16

# ~~~中略(前節のコード)~~~

rmf_filename = 'rsl_Hp_L.rmf'
out_filename = f'resolve_rmf_specific_ene_output.png'

# RMFデータの読み込みと構築
matrix, in_ene_lo, in_ene_hi, out_ene_lo, out_ene_hi = construct_rmf_matrix(rmf_filename)

# 各チャンネルの平均エネルギー値を計算
in_ene_mean = (in_ene_lo + in_ene_hi) / 2
out_ene_mean = (out_ene_lo + out_ene_hi) / 2

# プロット作成
fig, ax = plt.subplots(figsize=(18, 9))

# カラーマップの設定
vmin, vmax = 1, 20
cmap = plt.get_cmap("plasma")
norm = plt.Normalize(vmin=vmin, vmax=vmax)

for target_ene in [1, 5, 10, 15, 20]:
    closest_index = np.abs(in_ene_mean - target_ene).argmin()
    closest_ene = in_ene_mean[closest_index]
    color = cmap(norm(closest_ene))
    ax.plot(out_ene_lo, matrix[closest_index], label=f'{closest_ene:.0f}', color=color)

# Resolveは小さい値をlogスケールで表示すると見にくいためsymlogを活用
min_log_thresh = 1e-5
ax.set_yscale('symlog', linthresh=min_log_thresh)

ax.set_xlim(-0.5, 21.5)

ax.xaxis.set_major_locator(AutoLocator())
ax.xaxis.set_minor_locator(AutoMinorLocator())
ax.yaxis.set_major_locator(SymmetricalLogLocator(base=10.0, linthresh=min_log_thresh))
ax.yaxis.set_minor_locator(LogLocator(base=10.0, subs=np.arange(2, 10) * 0.1, numticks=10))

ax.grid(which='major', color='gray', linestyle='--', linewidth=0.5)
ax.grid(which='minor', color='lightgray', linestyle=':', linewidth=0.3)
ax.set_xlabel('Output (keV)')
ax.set_ylabel('Probability')
fig.legend(title='Input (keV)', loc='upper right', bbox_to_anchor=(0.99, 0.91))

# ファイル名をタイトルに設定し保存
ax.set_title('Resolve: Energy Dependence of RMF')
fig.savefig(out_filename, bbox_inches='tight')

plt.show()

Xtend

Resolve

ARFの可視化

ARFデータの特性

ARFには、エネルギーごとの有効面積（cm²）が格納されています。これは、観測装置の感度特性を示し、エネルギーごとにどれだけの光子を効率的に検出できるかを表します。有効面積は、観測装置の設計や構造、エネルギー帯域によって大きく異なります。

ARFの読み込みと可視化

astropy.io.fits を用いてARFを読み込み、有効面積データを可視化します。XtendとResolveのARFを比較することで、両装置の感度特性の違いが確認できます。

以下は、ARFから有効面積を抽出して可視化した結果の説明です。

ARFのデータ構造の確認

fstructでデータ構造を確認できます。

XtendのARF

fstruct xtd_source.arf
  No. Type     EXTNAME      BITPIX Dimensions(columns)      PCOUNT  GCOUNT
 
   0  PRIMARY                 16     0                           0    1
   1  BINTABLE SPECRESP        8     12(3) 6400                  0    1
 
      Column Name                Format     Dims       Units     TLMIN  TLMAX
      1 ENERG_LO                   E                   keV
      2 ENERG_HI                   E                   keV
      3 SPECRESP                   E                   cm^2

ResolveのARF

fstruct rsl_Hp.arf
  No. Type     EXTNAME      BITPIX Dimensions(columns)      PCOUNT  GCOUNT
 
   0  PRIMARY                 16     0                           0    1
   1  BINTABLE SPECRESP        8     12(3) 60000                 0    1
 
      Column Name                Format     Dims       Units     TLMIN  TLMAX
      1 ENERG_LO                   E                   keV
      2 ENERG_HI                   E                   keV
      3 SPECRESP                   E                   cm^2

RMFの読み込みと可視化

from astropy.io import fits
import matplotlib.pyplot as plt

# ARFのパスを指定
resolve_arf_file = "rsl_Hp.arf"
xtend_arf_file = "xtd_source.arf"

# ARFを読み込む関数
def load_arf(file_path):
    with fits.open(file_path) as hdul:
        arf_data = hdul[1].data
    energy_lo = arf_data['ENERG_LO']  # 入力の下限エネルギー
    energy_hi = arf_data['ENERG_HI']  # 出力の上限エネルギー
    sensitivity = arf_data['SPECRESP']  # 感度
    energy_mid = (energy_lo + energy_hi) / 2  # 中心エネルギー
    return energy_mid, sensitivity

# データをロード
resolve_energy, resolve_sensitivity = load_arf(resolve_arf_file)
xtend_energy, xtend_sensitivity = load_arf(xtend_arf_file)

# プロット
plt.figure(figsize=(8, 6))
plt.plot(xtend_energy, xtend_sensitivity, label="Xtend ARF", color="r")
plt.plot(resolve_energy, resolve_sensitivity, label="Resolve ARF", color="k")
plt.xlim(0, 14)
plt.xlabel("Energy (keV)")
plt.ylabel("Effective area (cm$^2$)")
plt.title("Comparison of Xtend and Resolve ARF")
plt.grid(True, which="both", linestyle="--", linewidth=0.5)
plt.legend()
plt.savefig("xtend_vs_resolve_arf.png")
plt.show()

出力結果

有効面積の比較
Resolveは全体的に有効面積が小さめですが、高エネルギー帯域でも一定の有効面積を維持しています。一方、Xtendは広視野で感度を最大化する設計のため、Resolveよりも有効面積が大きくなっています。
低エネルギー帯域での差異
Resolve では、低エネルギー（2 keV 以下）で有効面積が著しく減少していることが確認できます。これは、GV（Gate Valve⁴）が閉じている影響によるものです。
エッジやピークの構造
両装置の有効面積は、特定のエネルギー帯域で変化を示します。この変化は、観測装置のミラーやフィルターの特性を反映しており、例えば、望遠鏡に使用されている Au の M エッジや L エッジの影響を受けています。これにより、一部のエネルギー帯域では急激な減少（エッジの影響）や有効面積の最大値（ピーク）が生じます。

まとめ

本記事では、X線観測データ解析に必要な RMF（Response Matrix File）と ARF（Auxiliary Response File）について、その基本概念と可視化手法を解説しました。RMF は入射エネルギーと出力エネルギーの関係を記録し、観測装置のエネルギー応答を示します。一方、ARF はエネルギーごとの有効面積を表し、観測装置の感度特性を反映します。

また、Resolve と Xtend の RMF・ARF を可視化し、それぞれの特性や違いを比較しました。これにより、観測装置のレスポンスがどのようにデータに影響を与えるかを直感的に理解しやすくなります。

現在、XRISM 衛星の Early Release Data（初期公開データ）が一般向けに公開されており、どなたでも以下のデータにアクセスできます。

XRISM Early Release Data

ぜひ、実際のデータを活用し、X線観測データ解析に役立ててみてください。

補足

XRSIMガイド

XRISM ABC Guide： https://heasarc.gsfc.nasa.gov/docs/xrism/analysis/abc_guide/xrism_abc.html
XRISM Analysis：https://heasarc.gsfc.nasa.gov/docs/xrism/analysis/index.html
XRISM Proposal：https://heasarc.gsfc.nasa.gov/docs/xrism/proposals/index.html
XRISM Publication：https://heasarc.gsfc.nasa.gov/docs/xrism/results/xrism_publications.html
RMF・ARFの説明などPDF：https://heasarc.gsfc.nasa.gov/docs/xrism/analysis/workshops/doc_jan24/2_9_Mernier.pdf

検出器の説明

Resolveレスポンス特性：Megan E. Eckart et al. 2018
Xtendレスポンス特性： Inoue et al. 2016
Filter Wheel：de Vries, Cor P., et al.
Gate Valve：Midooka et al. 2020
Dewar の内側の filter：Caroline A. Kilbourne et al. 2018

Resolveのエネルギービンが60,000点ある理由

Resolveのエネルギービンは0.5 eV刻みで最大30 keVまであるため、60,000チャンネルになります。

これは XSPEC向けにeV単位のデフォルト値offset = 0.5, binwidth = 0.5で丸め込んでいるだけであり、実際のエネルギー情報は EPI2列にfloat精度で記録されています。

PI チャンネルは以下の式で整数化されます。

\text{PI} = \text{floor} \left( \frac{\text{EPI2} - 0.5}{0.5} + 1.0 \right)

したがって、キャリブレーション向けには EPI2 を使えば高精度なエネルギー情報が得られますが、通常の解析では PI を利用します。

参考: rslpha2piの公式ページ

Resolve RMFの種類

rslmkrmf (whichrmf=L) の他に、データ解析の目的に応じて様々なモードが用意されています。

現在、whichrmfは以下の4種類用意されています。

サイズ	オプション (`whichrmf`)	含まれる成分
Small	`S`	ガウシアンコアのみ
Medium	`M`	`S` + 指数関数的なテール + Si Kα線
Large	`L`	`M` + エスケープピーク
X-Large	`X`	`L` + 電子損失連続成分（ELC）

精密な解析が求められる場合には rslmkrmf (whichrmf=X) を使用できますが、データサイズが大きくなるため注意が必要です。
ここでは、whichrmf=X の注意点とLと比べた時の違いを可視化して確認します。

`whichrmf=X` の設定と注意点

whichrmf=X は、L のエスケープピークに加え、ELC（electron loss continumm） を考慮するため、より詳細な計算が可能です。しかし、データサイズが2GB を超えるため、処理が途中で停止することがあります。

実際に処理が停止した際に出るエラー

筆者が遭遇したエラーログ

INFO   ::rslrmf: Completed 44000 energy channels out of 60000
INFO   ::rslrmf: Completed 45000 energy channels out of 60000
INFO   ::rslrmf: Completed 46000 energy channels out of 60000
INFO   ::rslrmf: ERROR: During main: ahfits::Buffer::write() -> rsl_source_Hp_all_X.rmf[MATRIX] (row 46336): error writing column: N_CHAN [CFITSIO STATUS 412: datatype conversion overflow]
FATAL  ::Command failed : returned value 1

このログから、46,000 チャンネル目付近でエラーが発生していることが分かります。

この場合、生成ファイルはありますが、Python で解析しようとすると、途中までしか生成されていないことがわかります。

ValueError: could not broadcast input array from shape (0,) into shape (46331,)

この問題を回避するために、コア成分（comb） と ELC 成分（elc） を分割出力するオプションが用意されています。

筆者のモードXの参考コード

rslmkrmf infile=event.evt 
    outfileroot=rsl_Hp
    resolist="0"
    regmode=DET
    regionfile=NONE
    pixlist="0-11,13-26,28-35"
    eminin=0.0
    dein=0.5
    nchanin=60000
    useingrd=no
    eminout=0.0
    deout=0.5
    nchanout=60000
    clobber=yes
    splitrmf=yes  # 重要
    splitcomb=yes  # 重要
    elcbinfac=16  # 重要
    whichrmf=X

splitrmf=yes と splitcomb=yes を設定すると、Core 成分（comb）と ELC 成分（elc）を分割して出力できます。さらに、elcbinfac を指定することで ELC 成分の入力エネルギーを適切にビンまとめし、データサイズを抑えられます。以下で、elcbinfac について少し補足します。

`elcbinfac`の調整

ELC は 連続スペクトルを持つため、データ量が大きくなりやすい（非ゼロ要素が多い）特性があります。
そのため、elcbinfac（デフォルト: elcbinfac=32）を指定すると、入力チャンネルをビンまとめし、計算負荷を軽減できます。

例えば、上のコード（elcbinfac=16）を実行すると、以下のようなファイルが生成されます。

comb（Core 成分）
- 入力チャンネル数: 60,000
- 出力チャンネル数: 60,000
elc（ELC 成分）
- 入力チャンネル数: 3,750（elcbinfacに設定した値でビンまとめさることを確認）
- 出力チャンネル数: 60,000

各入力チャンネルに対して、comb と elc の合計が 1 になるように 設計されており、これらを組み合わせることで確率的な扱いが可能になります。
ただし、入力チャンネルのビンサイズが異なる点に注意が必要です。次節では、この違いを考慮し、可視化を行い各ファイルの特徴を確認します。

ELC 成分の可視化

ELC を加えた場合の影響を確認するために、可視化専用のコード（ここでは解析用途に使えるような厳密な目的ではないことに注意）を示します。

comb+elcの統合用の関数

def merge_x_rmf(matrix_comb, matrix_elc, elcbinfac=16):
    """
    rslmkrmf（whichrmf=X）用の
    core成分、ELC成分を統合した Matrix を返す関数。

    Parameters:
        matrix_comb (numpy.ndarray): 2D入力配列（Core 成分）。
        matrix_elc (numpy.ndarray): 2D入力配列（ELC 成分）。
        elcbinfac (int): ELC成分のbinsize（入力エネルギー方向のビニングサイズ）。

    Returns:
        numpy.ndarray: 統合された 2D 配列。
    """
    # matrix_comb の形状を取得
    rows, cols = matrix_comb.shape
    
    # bin_size の倍数でない場合、余った行を除外
    trimmed_rows = (rows // elcbinfac) * elcbinfac
    matrix_comb_trimmed = matrix_comb[:trimmed_rows, :]

    # y方向にbinning（平均を取る）
    matrix_comb_binned = matrix_comb_trimmed.reshape(-1, elcbinfac, cols).mean(axis=1)
    
    # matrix_comb_binned と matrix_elc を足し込み
    merged_matrix = matrix_comb_binned + matrix_elc

    return merged_matrix

combの入力の各ビンを平均することで、elc のビンに対応する平均的な特徴を計算し、そのまま足し込めるようにしています。

各成分の可視化

comb（Core 成分）	elc（ELC 成分）

core + elc（統合結果）

このように、図に示すことで ELC 成分の寄与を直感的に理解できます。

参考: rslmkrmfの公式ページ

Resolve: Line spread function ↩
Xtendレスポンス特性： Inoue et al. 2016 ↩
Resolveレスポンス特性（Cu Kα線を入射した際のレスポンスがFig. 14にあります。）：Megan E. Eckart et al. 2018 ↩
Gate Valve：Midooka et al. 2020 ↩