Pythonで解析するXRISM衛星のRMF・ARFの可視化

## **はじめに**
宇宙実験や天文学においては、観測装置の特性を理解し、正確なデータ解析を行うことが重要です。この記事では、特に、X線スペクトルの解析を例として、**観測データと物理モデルを比較するための応答関数（XPECという汎用X線解析ツールでRMF・ARFと呼ばれるもの）** について、XRISM衛星の例をつかって紹介したいと思います。  

[XRISM X線天文衛星](https://www.xrism.jaxa.jp)には、２つの観測装置（X線精密分光装置 Resolve とX線撮像装置 Xtend）を搭載しています。この記事では、Python を用いた RMF・ARF の可視化手法を通じて、各装置の特性や違いが自然に理解できるように紹介したいと思います。深い装置の理解は、新しいサイエンスの結果を見逃さないためにも役立つと思いますので、直感的に理解できる可視化技術も心かげています。  

XRISM衛星のARF・RMFの詳細について深く知りたい場合は、日本語記事を探すのではなく、下記の公式ページから情報を集めてください。

https://heasarc.gsfc.nasa.gov/docs/xrism/analysis/abc_guide/Resolve_Data_Analysis.html

https://heasarc.gsfc.nasa.gov/docs/xrism/analysis/abc_guide/Xtend_Data_Analysis.html


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

### **RMF（Response Matrix File）**
RMFは、Response Matrix File、の略語で、[xspec](https://heasarc.gsfc.nasa.gov/xanadu/xspec/)というX線天文のNASAが開発管理しているX線解析プラットフォームで定義されるものです。一般的な、観測応答や装置レスポンスと、xspec の terminology での RMF が対応しない場合もあるので注意が必要です。

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

```math
R_{ij} = P(E_i | E_j)
```

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

```math
\sum_{i} R_{ij} = 1
```

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

RMFファイルが一つの行列でよいとは限らない点にも注意が必要です。例えば、Resolveでは、36個の異なる個性をもつ画素があり、36種類の個性をCALDBというデータベースファイルに保存し、 

https://heasarc.gsfc.nasa.gov/lheasoft/help/rslmkrmf.html

というソフトウェアでユーザーが使う画素に応じて重みづけしたレスポンスファイルを生成する仕組みになっています。Resolve のカロリメータの場合は、吸収体のHgTeの厚みや接着の強さにも個性があり(全自動ではなくハンドマニピュレーターつける作業)、温度計としてのSi半導体カロリメータも分解能をよくするためにリンとドープしており(ref. https://ui.adsabs.harvard.edu/abs/2016SPIE.9905E..0VK/abstract) 36画素あり、場所による個性も僅かに存在し、全画素が同じレスポンス関数には原理的になりません。
一方で、XtendのCCDの方は、約164万画素(640x640x4chip)のあるのに、同じことをしてるかというとそうではなく、カロリメータよりも均質性がよいので画素ごとに作ることはしてません。

---

### **ARF（Auxiliary Response File）**
ARFは、Auxiliary Response Fileという、Xspec用語で、「アーフ」とジャーゴンで呼ばれますが、Auxiliary という英語は、"補助の"という意味しかないので、「アーフ」はX線天文以外の方には絶対に通じません。検出器応答には、検出された量を比較可能な物理量に直す必要があるので、入射エネルギーと出力エネルギーを繋ぐ確率遷移関数のような役割をするファイルとRMFとして定義しただけでは無次元なので物理量には直らないので、補助ファイル(ARF)が物理次元を保持して、観測量を物理量に繋げる役割をするものとして棲み分けましょう、という気持ちだと理解するとよいでしょう。

ARFは、X線望遠鏡がある場合は望遠鏡を含んだ検出器システム全体の有効面積（effective area）をエネルギー方向の一次元ベクトルで表したものです。各エネルギーにおける感度を示し、望遠鏡の集光効率や検出器の感度、X線入射窓、フィルター、保護膜などの効果を含んでいます。つまり、ARFは「入射エネルギー $E_{j}$ ごとの有効面積」という1次元の情報（$A_{j}$）を持ち、入射エネルギーに対する関数として表されます。

逆に言えば、ARFだけを補正するのであれば、単に$A_{j}$の逆数を観測データにかければよいだけですが、宇宙観測ではこのような逆演算には注意が必要です。例えば、有効面積がほぼ 0 に近いエネルギー帯域であれば、$A_{j}$の値が $10^{-5}$とか$10^{-6}$とか、どこまでも小さい値になり得ますが、実際の観測データには装置のノイズや non X-ray background があるため、0 ではなく、1 cts, 2cts など有限のカウント数がありえます。有限な観測量(主にノイズ) を小さい有効面積で割ると大きな観測値になりますが、これは単に有効面積の計算では考慮してない background成分が存在することを意味しています。

---

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

```math
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）](https://heasarc.gsfc.nasa.gov/docs/software/heasoft/)のパッケージの一部として提供されているツールで、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
```

<details><summary>ResolveのRMF</summary>

```: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
```

</details>

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

```python:トイデータでスパース行列の可視化
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()
```
**出力結果**

![sparse_to_full_rmf.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/244df0da-6130-f531-a23a-1decda33a894.png)

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

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

**コードの説明**
1. **`j`（ループのインデックス）**  
   - **入力チャンネル**を指します（例: 0番目のチャンネル、1番目のチャンネルなど）。
   - `f_chan[j]` と `n_chan[j]` は、このビンに対応する出力チャンネル情報を持っています。

2. **`f_chan[j]`**  
   - この入力チャンネルにおける**出力チャンネルの開始位置**を示します。  
   - 例: `f_chan[j] = [2, 6]` なら、このビンの非ゼロ応答データはチャンネル2から開始し、次にチャンネル6から開始する応答が続いています。

3. **`n_chan[j]`**  
   - **`f_chan[j]` の各開始位置から、何出力チャンネル分の応答を含んでいるか**を示しています。  
   - 例:  `f_chan[j] = [2, 6]`で、`n_chan[j] = [2, 3]` なら、出力チャンネル2から始まり2チャンネル分、チャンネル6から始まり、3チャンネル分の応答を含むことを意味します。

4. **`compressed_matrix[j]`**  
   - スパース形式で格納された応答データ（非ゼロ部分のみ）。  
   - `tot_n` を利用して、現在どの部分のスパースデータを展開すべきかを管理しています。

:::note
**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  

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

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

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

```python: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 fits

# 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 fits.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`）
![xtend_rmf_matrix_bin_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/a0ae29ea-f267-e971-0a23-fffd01400709.png)

**Resolve**（RMF: `rslmkrmf (whichrmf=L)`, 出力: `bin_size=200`）
（`binsize`が小さいと描画に時間がかかるので注意）
![resolve_rmf_matrix_bin_200_L.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/155764c2-b99c-f02f-8573-5a97caf2df44.png)


**RMFマトリックスの説明**
- **対角線の強度**: 入射エネルギーと出力エネルギーが一致する対角線上には強いピークが見られます。これらのピークは、検出器が入力エネルギーを正確に検出した際のスペクトルのエネルギー分解能を反映しており、LSF（Line spread function[^1]）を示しています。特に、Xtend に比べて Resolve はエネルギー分解能が桁で高く、よりシャープなピークを持つことがわかります。LSFは、は、検出器の分解能や材質の特性に依存します。特に、エネルギーが高くなるにつれて検出効率が低下し、分布が広がる傾向が確認できます。
- **エスケープピーク**: Xtend・Resolve の両方で、対角成分よりも少し低いエネルギー帯に斜め方向の分布が見られます。これは、Si 検出器特有のエスケープピークによるものです[^2]。さらに、Resolve では、検出器の  HgのM線やTeのL線などに起因するエスケープピークも確認できます[^3]。

- **エネルギー分解能**: 分解能が低い場合、対角線から外れる成分が広がる傾向があります。Resolveのような高分解能検出器ではピークがシャープに集中しますが、Xtendでは分解能が低いため、広がったピークが現れます。この違いにより、両者のエネルギー応答特性を可視化で比較することが可能です。

- **非ゼロ要素の分布**: 非ゼロ要素の分布は、検出器が観測可能なエネルギー範囲を示します。Siを用いた検出器のため、1.7 keV付近のSiの特性X線も顕著に見られます。

:::note
**元素成分のデータベース**
- X線スペクトルに含まれる輝線成分は、NISTの[X-ray Transition Energies Database](https://physics.nist.gov/PhysRefData/XrayTrans/Html/search.html)などのデータベースで調べることができます。
- 例えば、Siの特性X線には以下のような輝線があります：
  - Kα2 (KL2): 1.7397 keV  
  - Kα1 (KL3): 1.7412 keV  
  - Kβ3 (KM2): 1.8418 keV  
- 一般に、Kα : Kβ の強度比はおおよそ 10:1 であり、主に Kα1・Kα2 が観測されます。
- カロリメータの温度計部分や基板にはSiが多く使用されているため、入射X線や宇宙線によって励起されたSiの蛍光X線が観測されます。
:::

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

![xtend_resolve_rmf_specific_ene_output.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/cfedea92-11c7-37e3-39f5-60edaf8c54bc.png)
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}$以下）はリニア、大きい成分は対数**で表示 
- **カラースケール** → 入力エネルギーごとに色分けし、分布の違いを視覚化  

<details><summary>RMFエネルギー依存性の可視化コード</summary>


```python: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 fits

# 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()
```

```python: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 fits

# 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()
```

---

</details>

**Xtend**
![xtend_rmf_specific_ene_output.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/fe0a374b-b57e-0806-5134-3f001a27fb87.png)

**Resolve**
![resolve_rmf_specific_ene_output.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/624edd60-748d-f7d2-4230-4d0e565782e3.png)

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

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

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

<details><summary>ARFのデータ構造の確認</summary>

`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
```

---

</details>



```python: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()
```
**出力結果**
![xtend_vs_resolve_arf.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/71c731e7-7197-35dc-7215-b203b2744fc5.png)

- **有効面積の比較**  
  Resolveは全体的に有効面積が小さめですが、高エネルギー帯域でも一定の有効面積を維持しています。一方、Xtendは広視野で感度を最大化する設計のため、Resolveよりも有効面積が大きくなっています。

- **低エネルギー帯域での差異**  
  Resolve では、低エネルギー（2 keV 以下）で有効面積が著しく減少していることが確認できます。これは、GV（Gate Valve[^4]）が閉じている影響によるものです。

- **エッジやピークの構造**  
  両装置の有効面積は、特定のエネルギー帯域で変化を示します。この変化は、観測装置のミラーやフィルターの特性を反映しており、例えば、望遠鏡に使用されている 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](https://heasarc.gsfc.nasa.gov/docs/xrism/results/erdata/index.html)  

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

実際にXRISMの実データのスペクトル解析の記事をご用意しましたのでよろしければご覧ください。

https://qiita.com/yusuke_s_yusuke/items/15f32f780c8a83b970a5

---

## 補足

### 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](https://doi.org/10.1117/1.JATIS.4.2.021406)
- Xtendレスポンス特性： [Inoue et al. 2016](https://doi.org/10.1016/j.nima.2016.03.071)
- Filter Wheel：[de Vries, Cor P., et al.](https://doi.org/10.1117/1.JATIS.4.1.011204)
- Gate Valve：[Midooka et al. 2020](https://doi.org/10.1117/12.2559451)
- Dewar の内側の filter：[Caroline A. Kilbourne et al. 2018](https://doi.org/10.1117/1.JATIS.4.1.011215)

---

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

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

PI チャンネルは以下の式で整数化されます。
```math
\text{PI} = \text{floor} \left( \frac{\text{EPI2} - 0.5}{0.5} + 1.0 \right)
```
したがって、キャリブレーション向けには EPI2 を使えば高精度なエネルギー情報が得られますが、通常の解析では PI を利用します。


参考: [`rslpha2pi`の公式ページ](https://heasarc.gsfc.nasa.gov/docs/software/lheasoft/help/rslpha2pi.html)

---

### 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 を超えるため、処理が途中で停止することがあります。

<details><summary>実際に処理が停止した際に出るエラー</summary>

```:筆者が遭遇したエラーログ
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,)
```

----


</details>



この問題を回避するために、**コア成分（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 を加えた場合の影響を確認するために、**可視化専用**のコード（ここでは解析用途に使えるような厳密な目的ではないことに注意）を示します。

```python: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 のビンに対応する平均的な特徴を計算し、そのまま足し込めるようにしています。

**各成分の可視化**


<table border="1">
    <thead>
        <tr>
            <th style="text-align:center;">comb（Core 成分）</th>
            <th style="text-align:center;">elc（ELC 成分）</th>
        </tr>
    </thead>
    <tbody>
        <tr>
            <td style="text-align:center;">
                <img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/f89bb434-812e-6e5f-8cd2-b3fb86d3c5bc.png" alt="comb">
            </td>
            <td style="text-align:center;">
                <img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/36af6b6f-c648-49e5-4b7c-49df418d85ac.png" alt="elc">
            </td>
        </tr>
    </tbody>
    <thead>
        <tr>
            <th colspan="2" style="text-align:center;">core + elc（統合結果）</th>
        </tr>
    </thead>
    <tbody>
        <tr>
            <td colspan="2" style="text-align:center;">
                <img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/52865b11-5b13-f9b2-45b3-28644c9352d7.png" alt="core+elc">
            </td>
        </tr>
    </tbody>
</table>


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

参考: [`rslmkrmf`の公式ページ](https://heasarc.gsfc.nasa.gov/lheasoft/help/rslmkrmf.html)

[^1]: [Resolve: Line spread function](https://heasarc.gsfc.nasa.gov/docs/xrism/proposals/POG/Resolve.html#sec:resolve_LSF)
[^2]: Xtendレスポンス特性： [Inoue et al. 2016](https://doi.org/10.1016/j.nima.2016.03.071)
[^3]: Resolveレスポンス特性（Cu Kα線を入射した際のレスポンスがFig. 14にあります。）：[Megan E. Eckart et al. 2018](https://doi.org/10.1117/1.JATIS.4.2.021406)
[^4]: Gate Valve：[Midooka et al. 2020](https://doi.org/10.1117/12.2559451)