はじめに
1999年に打ち上げられたChandraは、25年以上にわたり高分解能X線分光を支えてきました。その中心的な役割を担ってきた装置の一つが、High Energy Transmission Grating（HETG）です。
HETGは、回折格子によってX線を波長ごとに分散させる分散型のX線分光装置です。通常は±1次光が広く用いられますが、2次光や3次光も存在し、一般に回折次数が高いほど高いエネルギー分解能が得られます。
そして2023年にはXRISMが打ち上げられ、新しい高分解能X線分光の時代が始まりました。XRISMに搭載されたResolveは、X線光子一つひとつのエネルギーを温度上昇として測定するX線マイクロカロリメータです。回折格子を用いない非分散型の装置でありながら、広いエネルギー帯域で5 eVを下回る高いエネルギー分解能を実現しています。
このような精密X線分光によって複雑な輝線構造を分離することで、高温プラズマの温度や元素組成だけでなく、ガスの速度、乱流、電離状態などを詳しく調べることができます。またResolveは、高いエネルギー分解能に加えて比較的大きな有効面積を持ち、高い統計精度を得やすいことも大きな科学的強みです。
HETGとResolveは、いずれも高分解能X線分光を実現する装置ですが、その測定原理やエネルギー依存性は大きく異なります。両者の性能はFWHMなどの数値として比較されることがありますが、実際のResponse Matrix File（RMF）に含まれるLine Spread Function（LSF）を、回折次数や入力エネルギーを変えながら同じ図上で比較する機会は、意外と多くありません。
では、HETGの1次光・高次光とResolveでは、実際のLSFはどのように異なるのでしょうか。
この記事では、実際の観測データに対して作成されたRMFを読み込み、その中身を可視化して比較します。比較するのは次のレスポンスです。
- Chandra/HETG HEGの1次・2次・3次光
- Chandra/HETG MEGの1次・2次・3次光
- XRISM/Resolve
具体的には、次の3点を確認します。
- 6.395 keVにおけるLSF
- 0.5 keV刻みで見たLSFの変化
- 入力エネルギーとFWHMの関係
ここで比較するのは、特定の観測データに対して作成されたRMFです。ミッション全体を代表する絶対的な性能値ではなく、実際の解析で使用するRMFを用いた一例としてご覧ください。
エネルギー分解能を決めるもの
XRISM/ResolveとChandra/HETGでは、X線光子のエネルギーを測定する仕組みが異なります。まずは、その点について簡単に整理します。
XRISM/Resolve
Resolveは、入射したX線光子のエネルギーを熱に変換し、その温度上昇からエネルギーを測定するX線マイクロカロリメータです。最も単純には、
\Delta T \simeq \frac{E}{C}
と表されます。ここで、$E$は入射した光子のエネルギー、$C$は検出器の熱容量、$\Delta T$は光子吸収によって生じる温度上昇です。
熱容量が小さいほど、同じエネルギーの光子に対して大きな温度上昇が得られます。ただし、実際のエネルギー分解能には、温度計の感度、熱揺らぎ、電気ノイズ、パルス処理なども関わります。Resolveは、これらの微小な温度変化を測定するため、約50 mKという極低温で動作しています。
Resolveのエネルギー分解能は、一般に「約5 eV」と表現されます。ただし、RMFに含まれるFWHMは完全な定数ではなく、高エネルギー側で緩やかに増加します。
ResolveのGaussian coreのFWHMは、概念的には、
\Delta E(E)
=
\sqrt{
\Delta E_{\mathrm{B}}^2
+
\left(C_{\mathrm{x}}E\right)^2
}
と表されます。ここで、$\Delta E_{\mathrm{B}}$は、熱揺らぎや電気雑音などに由来する、エネルギーにほとんど依存しないbaseline resolutionです。一方、$C_{\mathrm{x}}E$は、吸収体内での熱化やパルス形状など、イベントごとの応答の微小な割合のばらつきによって生じると考えられるexcess broadeningです（詳しくは、山田さんのこちらの記事やLeutenegger et al. 2025 を参照）。
同じ割合のばらつきでも、eV単位の誤差は入射エネルギーに比例して大きくなります。そのため、ResolveのFWHMは高エネルギー側でわずかに増加します。
Chandra/HETG
HETGでは、回折格子によってX線を波長ごとに分散します。回折条件は
m\lambda=p\sin\beta
で表されます。
ここで、
- $m$：回折次数
- $\lambda$：X線の波長
- $p$：回折格子の周期（ピッチ）
- $\beta$：回折角
です。
HETGでは、検出器上の位置から波長を求めています。したがって、位置をどれだけ正確に決められるかが波長分解能を決めます。
この関係は概念的に
\Delta\lambda_m\propto\frac{p}{|m|}
と表せます。
したがって、回折次数$|m|$が高いほど、また格子周期$p$が小さいほど、波長幅$\Delta\lambda$は小さくなります。HETGを構成するHEGはMEGよりも格子周期が小さいため、同じ回折次数ではHEGの方が高い波長分解能を持ちます。
一方、エネルギーと波長の関係は
E=\frac{hc}{\lambda}
です。ここで、
- $E$：X線のエネルギー
- $h$：プランク定数
- $c$：光速
です。
この式を微分すると、エネルギー幅と波長幅の関係は
\Delta E\simeq\frac{E^2}{hc}\Delta\lambda
となります。
したがって、波長幅$\Delta\lambda$がほぼ一定であれば、エネルギー単位のFWHMは高エネルギーほど大きくなります。さらに、格子周期と回折次数への依存性を含めると、点源に対する概念的な関係は
\Delta E_m\propto\frac{pE^2}{|m|}
となります。
使用するRMF
今回は、いずれもSS 433の観測から作成したRMFを使用します。
XRISM/Resolve
- 天体：SS 433
- ObsID：300041010
- 観測期間：2024年4月10日〜15日
- 使用ピクセル：較正用ピクセルを除く全ピクセル
- イベントグレード：Hp
- RMF：
whichrmf=L
Chandra/HETG
- 天体：SS 433
- ObsID：106
- 観測時期：1999年9月23日
- 使用レスポンス：HEG・MEGの|m|=1,2,3（正負回折次数を結合）
ObsID 106は、Chandra/HETGSによるSS 433初期観測の一つです。
本記事で比較するのは、これらの観測用RMFを用いた一例です。実際の解析では、対象データに対応したRMFを使用してください。
今回のデータは以下から取得できます。
エネルギー分解能の評価方法
RMFは、真のエネルギーを持つ光子が検出器上でどのエネルギーに記録されるかを表す行列です。
ある入力エネルギーに対応する1行を取り出すことで、その単色X線に対するLSFを得ることができます。
今回は、RMFの値を出力エネルギービン幅で割り、
P_E=\frac{P_i}{\Delta E_i}
として確率密度（keV$^{-1}$）で描画します。こうすることで、出力ビン幅が異なるRMF同士でもLSFの形状を直接比較できます。
入力エネルギーとRMFの入力ビン中心は必ずしも一致しないため、本記事では指定したエネルギーを含む入力ビンを選択し、その入力ビン中心エネルギーも併せて表示します。
また、LSFをFWHM（Full Width at Half Maximum）で評価します。FWHMは、主ピークの最大値の半分（Half Maximum）となる左右のエネルギーを線形補間によって求め、その差として計算しています。副次的なピークやescape peakなどの構造は含めず、主ピークのみを対象としています。
6.395 keVのLSFを比較する
まず入力エネルギーを6.395 keVとし、各RMFのLSFを比較します。6.4 keV付近には低電離鉄のFe Kα線があり、X線分光で特に重要なエネルギー帯です。
コードはこちら
import os
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
# ============================================================
# Matplotlib settings
# ============================================================
plt.rcParams["xtick.direction"] = "in"
plt.rcParams["ytick.direction"] = "in"
plt.rcParams["font.size"] = 12
# ============================================================
# RMF matrix construction
# This function is intentionally unchanged.
# ============================================================
def construct_rmf_matrix(filename, ext=0):
"""
RMFからレスポンスマトリックスを構築する関数（Chandra用）。
Parameters:
filename (str): FITSファイルのパス。
ext: indexの始まり
Returns:
dict: 構築されたレスポンスマトリックスと関連データを含む辞書。
"""
with fits.open(filename) as f:
# 必要なデータを抽出
f_chan = f[1].data['F_CHAN']-ext # チャンネルの開始位置（Chandraは1始まり、XRISMは0始まりなので注意）
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
# ============================================================
# Input-energy bin selection
# ============================================================
def find_input_energy_bin(in_ene_lo, in_ene_hi, target_energy):
"""Return the RMF input-bin index containing target_energy."""
indices = np.where(
(in_ene_lo <= target_energy)
& (target_energy < in_ene_hi)
)[0]
# Include the final upper boundary when target_energy is exactly equal to it.
if len(indices) == 0 and np.isclose(target_energy, in_ene_hi[-1]):
return len(in_ene_hi) - 1
if len(indices) == 0:
return None
return int(indices[0])
# ============================================================
# FWHM calculation
# ============================================================
def calc_fwhm(energy, response_density):
"""
Calculate the FWHM of the main response peak using linear interpolation.
Only the contiguous half-maximum region containing the global maximum is
used. This prevents a distant escape peak or secondary response component
from artificially broadening the measured FWHM.
Parameters
----------
energy : numpy.ndarray
Output-energy bin centers [keV].
response_density : numpy.ndarray
Response probability density [keV^-1].
Returns
-------
float
FWHM [keV]. Returns np.nan when it cannot be determined.
"""
energy = np.asarray(energy, dtype=float)
response_density = np.asarray(response_density, dtype=float)
valid = np.isfinite(energy) & np.isfinite(response_density)
energy = energy[valid]
response_density = response_density[valid]
if len(energy) < 3:
return np.nan
peak_index = int(np.argmax(response_density))
peak_value = response_density[peak_index]
if not np.isfinite(peak_value) or peak_value <= 0:
return np.nan
half_maximum = 0.5 * peak_value
left_inside = peak_index
while (
left_inside > 0
and response_density[left_inside - 1] >= half_maximum
):
left_inside -= 1
right_inside = peak_index
while (
right_inside < len(response_density) - 1
and response_density[right_inside + 1] >= half_maximum
):
right_inside += 1
if left_inside == 0 or right_inside == len(response_density) - 1:
return np.nan
def interpolate_crossing(i1, i2):
x1, x2 = energy[i1], energy[i2]
y1, y2 = response_density[i1], response_density[i2]
if np.isclose(y1, y2):
return 0.5 * (x1 + x2)
return x1 + (
(half_maximum - y1)
* (x2 - x1)
/ (y2 - y1)
)
left_crossing = interpolate_crossing(left_inside - 1, left_inside)
right_crossing = interpolate_crossing(right_inside, right_inside + 1)
fwhm = right_crossing - left_crossing
return fwhm if fwhm > 0 else np.nan
# ============================================================
# Shared plotting definitions
# ============================================================
def build_series_definitions(chandra_rmf_directory, resolve_rmf_filename):
"""Build consistent labels and plotting styles for all three scripts."""
colors = {1: "C0", 2: "C1", 3: "C2"}
definitions = []
for abs_order in (1, 2, 3):
definitions.append(
{
"key": f"heg{abs_order}",
"label": f"HEG |m|={abs_order}",
"short_label": f"HEG{abs_order}",
"csv_label": f"HEG_m{abs_order}_FWHM_eV",
"filename": os.path.join(
chandra_rmf_directory,
f"combo_heg_abs{abs_order}.rmf",
),
"ext": 1,
"instrument": "HEG",
"order": abs_order,
"color": colors[abs_order],
"linestyle": "-",
"marker": "o",
"linewidth": 1.7,
"alpha": 1.0,
}
)
definitions.append(
{
"key": f"meg{abs_order}",
"label": f"MEG |m|={abs_order}",
"short_label": f"MEG{abs_order}",
"csv_label": f"MEG_m{abs_order}_FWHM_eV",
"filename": os.path.join(
chandra_rmf_directory,
f"combo_meg_abs{abs_order}.rmf",
),
"ext": 1,
"instrument": "MEG",
"order": abs_order,
"color": colors[abs_order],
"linestyle": ":",
"marker": "s",
"linewidth": 1.7,
"alpha": 0.8,
}
)
definitions.append(
{
"key": "resolve",
"label": "XRISM Resolve",
"short_label": "Resolve",
"csv_label": "Resolve_FWHM_eV",
"filename": resolve_rmf_filename,
"ext": 0,
"instrument": "Resolve",
"order": None,
"color": "black",
"linestyle": "-",
"marker": "D",
"linewidth": 2.3,
"alpha": 1.0,
}
)
return definitions
# ============================================================
# Extract one RMF response
# ============================================================
def load_response_density(rmf_filename, target_energy, ext, label):
"""Load one RMF and extract the response at target_energy."""
if not os.path.isfile(rmf_filename):
raise FileNotFoundError(f"{label}: RMF file not found:\n{rmf_filename}")
matrix, in_lo, in_hi, out_lo, out_hi = construct_rmf_matrix(
rmf_filename,
ext=ext,
)
input_index = find_input_energy_bin(in_lo, in_hi, target_energy)
if input_index is None:
raise ValueError(
f"{label}: target energy {target_energy:.6f} keV is outside "
"the RMF input-energy range."
)
response = np.asarray(matrix[input_index], dtype=float)
edges = np.append(out_lo, out_hi[-1])
centers = 0.5 * (out_lo + out_hi)
widths = out_hi - out_lo
if np.any(widths <= 0):
raise ValueError(f"{label}: non-positive output-energy bin width found.")
density = response / widths
fwhm = calc_fwhm(centers, density)
input_center = 0.5 * (in_lo[input_index] + in_hi[input_index])
input_offset_ev = (input_center - target_energy) * 1.0e3
probability_sum = np.sum(response)
print(f"[{label}]")
print(f" target energy : {target_energy:.6f} keV")
print(f" input-bin range : [{in_lo[input_index]:.6f}, {in_hi[input_index]:.6f}] keV")
print(f" input-bin center : {input_center:.6f} keV")
print(f" center offset : {input_offset_ev:+.2f} eV")
print(f" response integral: {probability_sum:.6f}")
if np.isfinite(fwhm):
print(f" main-peak FWHM : {fwhm * 1.0e3:.2f} eV")
else:
print(" main-peak FWHM : N/A")
print()
return edges, density, fwhm
# ============================================================
# Main
# ============================================================
def main():
target_energy = 6.395 # keV
x_half_width = 0.05 # keV
output_filename = "rmf_response_6p395kev_heg_meg_abs123_resolve.png"
chandra_rmf_directory = (
"ss433/chandra/"
"repro/merge_all_tg"
)
resolve_rmf_filename = (
"ss433/resolve/"
"xa300041010rsl_source_Hp_all.rmf"
)
definitions = build_series_definitions(
chandra_rmf_directory,
resolve_rmf_filename,
)
fig, ax = plt.subplots(figsize=(7.4, 5.2))
for definition in definitions:
edges, density, fwhm = load_response_density(
rmf_filename=definition["filename"],
target_energy=target_energy,
ext=definition["ext"],
label=definition["label"],
)
fwhm_label = (
f"{fwhm * 1.0e3:.1f} eV"
if np.isfinite(fwhm)
else "N/A"
)
ax.step(
edges[:-1],
density,
where="post",
color=definition["color"],
linestyle=definition["linestyle"],
linewidth=definition["linewidth"],
alpha=definition["alpha"],
label=f'{definition["label"]} (FWHM={fwhm_label})',
)
ax.axvline(
target_energy,
color="0.5",
linestyle="--",
linewidth=0.8,
alpha=0.7,
zorder=0,
)
ax.set_xlim(target_energy - x_half_width, target_energy + x_half_width)
ax.set_ylim(bottom=0)
ax.set_xlabel("Output Energy (keV)")
ax.set_ylabel(r"Probability Density (keV$^{-1}$)")
ax.set_title(
f"RMF response at {target_energy:.3f} keV: "
"Chandra HETG vs XRISM Resolve",
pad=5,
)
ax.grid(visible=True, which="major", alpha=0.4)
ax.tick_params(which="both", top=True, right=True, pad=3)
ax.minorticks_on()
ax.legend(
fontsize="x-small",
ncol=1,
frameon=True,
handlelength=2.8,
)
fig.tight_layout(pad=0.7)
fig.savefig(output_filename, dpi=200, bbox_inches="tight", pad_inches=0.04)
print(f"Saved figure: {output_filename}")
plt.show()
if __name__ == "__main__":
main()
図では、実線がHEG、点線がMEG、色が回折次数、黒線がResolveを表しています。
このエネルギー帯域では、Resolveが最も高いエネルギー分解能を示しています。 HETGでは、回折次数が高いほどLSFの主ピークが細くなり、さらにHEGはMEGよりも分解能が高いことが分かります。
他のエネルギーでも確認する
続いて、入力エネルギーを0.5〜12 keVまで0.5 keV刻みで変化させ、各LSFをタイル状に表示します。
コードはこちら
import os
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
# ============================================================
# Matplotlib settings
# ============================================================
plt.rcParams["xtick.direction"] = "in"
plt.rcParams["ytick.direction"] = "in"
plt.rcParams["font.size"] = 12
# ============================================================
# RMF matrix construction
# ============================================================
def construct_rmf_matrix(filename, ext=0):
"""
RMFからレスポンスマトリックスを構築する関数（Chandra用）。
Parameters:
filename (str): FITSファイルのパス。
ext: indexの始まり
Returns:
dict: 構築されたレスポンスマトリックスと関連データを含む辞書。
"""
with fits.open(filename) as f:
# 必要なデータを抽出
f_chan = f[1].data['F_CHAN']-ext # チャンネルの開始位置（Chandraは1始まり、XIRMは0始まりなので注意）
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
# ============================================================
# Optional output-energy binning
# ============================================================
def binning_out_ene(matrix, bin_size, ene_lo, ene_hi):
"""
2D配列の列方向をビニングし、出力エネルギー境界も更新する。
Parameters
----------
matrix : numpy.ndarray
RMFの2次元レスポンスマトリックス。
bin_size : int
まとめる出力エネルギービンの数。
ene_lo : numpy.ndarray
出力エネルギービンの下限。
ene_hi : numpy.ndarray
出力エネルギービンの上限。
Returns
-------
binned_matrix : numpy.ndarray
列方向にビニングしたレスポンスマトリックス。
binned_out_ene_edges : numpy.ndarray
ビニング後の出力エネルギー境界。
"""
if bin_size < 1:
raise ValueError("bin_size must be a positive integer.")
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
# ============================================================
# Input-energy bin selection
# ============================================================
def check_selected_bin(in_lo, in_hi, target_energy, label=""):
"""
指定した入力エネルギーを含むRMF入力ビンを返す。
Parameters
----------
in_lo, in_hi : numpy.ndarray
RMF入力エネルギービンの下限と上限。
target_energy : float
選択する入力エネルギー [keV]。
label : str
標準出力に表示するラベル。
Returns
-------
idx : int
指定エネルギーを含む入力ビンのインデックス。
"""
selected = np.where(
(in_lo <= target_energy) & (in_hi > target_energy)
)[0]
if len(selected) == 0:
raise ValueError(
f"{label}: target energy {target_energy:.6f} keV "
"is outside the RMF input-energy range."
)
idx = selected[0]
lo = in_lo[idx]
hi = in_hi[idx]
center = 0.5 * (lo + hi)
offset = center - target_energy
print(f"[{label}]")
print(f" target E : {target_energy:.6f} keV")
print(f" bin index : {idx}")
print(f" bin range : [{lo:.6f}, {hi:.6f}] keV")
print(f" bin center : {center:.6f} keV")
print(f" offset : {offset * 1e3:+.2f} eV")
print()
return idx
# ============================================================
# FWHM calculation
# ============================================================
def calc_fwhm(energy, response_density):
"""
レスポンス曲線のFWHMを線形補間によって計算する。
Parameters
----------
energy : numpy.ndarray
出力エネルギービンの中心 [keV]。
response_density : numpy.ndarray
レスポンスの確率密度 [keV^-1]。
Returns
-------
float
FWHM [keV]。計算できない場合はnp.nan。
"""
energy = np.asarray(energy, dtype=float)
response_density = np.asarray(response_density, dtype=float)
valid = np.isfinite(energy) & np.isfinite(response_density)
energy = energy[valid]
response_density = response_density[valid]
if len(energy) < 3:
return np.nan
maximum = np.max(response_density)
if maximum <= 0:
return np.nan
half_maximum = 0.5 * maximum
above_half = np.where(response_density >= half_maximum)[0]
if len(above_half) == 0:
return np.nan
left_inside = above_half[0]
right_inside = above_half[-1]
# 左右の半値交点を補間するため、曲線の内外に点が必要
if left_inside == 0 or right_inside >= len(energy) - 1:
return np.nan
def interpolate_crossing(i1, i2):
x1 = energy[i1]
x2 = energy[i2]
y1 = response_density[i1]
y2 = response_density[i2]
if np.isclose(y1, y2):
return 0.5 * (x1 + x2)
return x1 + (
(half_maximum - y1)
* (x2 - x1)
/ (y2 - y1)
)
left_crossing = interpolate_crossing(
left_inside - 1,
left_inside,
)
right_crossing = interpolate_crossing(
right_inside,
right_inside + 1,
)
return right_crossing - left_crossing
# ============================================================
# Extract one RMF response
# ============================================================
def load_response_density(
rmf_filename,
target_energy,
ext,
label,
):
"""
RMFから指定入力エネルギーに対する出力レスポンスを取り出す。
Parameters
----------
rmf_filename : str
RMFファイル。
target_energy : float
入力エネルギー [keV]。
ext : int
F_CHANに対して引く値。
Chandra RMFでは通常1、Resolve RMFでは通常0。
label : str
標準出力用ラベル。
Returns
-------
edges : numpy.ndarray
出力エネルギービン境界 [keV]。
centers : numpy.ndarray
出力エネルギービン中心 [keV]。
density : numpy.ndarray
確率密度 [keV^-1]。
fwhm : float
FWHM [keV]。
"""
if not os.path.isfile(rmf_filename):
raise FileNotFoundError(
f"{label}: RMF file not found:\n{rmf_filename}"
)
matrix, in_lo, in_hi, out_lo, out_hi = construct_rmf_matrix(
rmf_filename,
ext=ext,
)
idx = check_selected_bin(
in_lo,
in_hi,
target_energy,
label=label,
)
response = np.asarray(matrix[idx], dtype=float)
edges = np.append(out_lo, out_hi[-1])
centers = 0.5 * (out_lo + out_hi)
widths = out_hi - out_lo
if np.any(widths <= 0):
raise ValueError(
f"{label}: non-positive output-energy bin width was found."
)
density = response / widths
fwhm = calc_fwhm(centers, density)
probability_sum = np.sum(response)
density_integral = np.sum(density * widths)
print(f" response sum : {probability_sum:.6f}")
print(f" density integral: {density_integral:.6f}")
if np.isfinite(fwhm):
print(f" FWHM : {fwhm * 1e3:.2f} eV")
else:
print(" FWHM : could not be calculated")
print()
return edges, centers, density, fwhm
import os
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
# ============================================================
# 複数の入力エネルギーについてレスポンスを取得
# ============================================================
def collect_rmf_responses(
rmf_filename,
target_energies,
ext,
label,
):
"""
1つのRMFを一度だけ読み込み、複数の入力エネルギーに対する
出力レスポンスとFWHMを取得する。
Parameters
----------
rmf_filename : str
RMFファイルのパス。
target_energies : array-like
評価する入力エネルギー [keV]。
ext : int
F_CHANから差し引く値。
Chandraでは通常1、Resolveでは通常0。
label : str
標準出力に表示する名前。
Returns
-------
dict
edges:
出力エネルギービン境界 [keV]。
responses:
各target energyに対応するレスポンス情報のリスト。
RMF範囲外の場合はNone。
"""
if not os.path.isfile(rmf_filename):
raise FileNotFoundError(
f"{label}: RMF file not found:\n{rmf_filename}"
)
print("=" * 72)
print(f"Loading: {label}")
print(rmf_filename)
matrix, in_lo, in_hi, out_lo, out_hi = construct_rmf_matrix(
rmf_filename,
ext=ext,
)
out_edges = np.append(out_lo, out_hi[-1])
out_centers = 0.5 * (out_lo + out_hi)
out_widths = out_hi - out_lo
if np.any(out_widths <= 0):
raise ValueError(
f"{label}: non-positive output-energy bin width was found."
)
responses = []
for target_energy in target_energies:
indices = np.where(
(in_lo <= target_energy)
& (in_hi > target_energy)
)[0]
# RMF入力エネルギー範囲外
if len(indices) == 0:
responses.append(None)
print(
f" {target_energy:4.1f} keV : "
"outside RMF input-energy range"
)
continue
input_index = int(indices[0])
response = np.asarray(
matrix[input_index],
dtype=float,
)
# 確率密度 [keV^-1]
response_density = response / out_widths
fwhm = calc_fwhm(
out_centers,
response_density,
)
input_bin_center = 0.5 * (
in_lo[input_index]
+ in_hi[input_index]
)
responses.append(
{
"density": response_density,
"fwhm": fwhm,
"input_bin_center": input_bin_center,
}
)
if np.isfinite(fwhm):
print(
f" {target_energy:4.1f} keV : "
f"FWHM = {fwhm * 1e3:8.2f} eV "
f"(input-bin center = "
f"{input_bin_center:.6f} keV)"
)
else:
print(
f" {target_energy:4.1f} keV : "
"FWHM could not be calculated"
)
print()
return {
"edges": out_edges,
"responses": responses,
}
# ============================================================
# Grid layout
# ============================================================
def choose_panel_grid(
n_panels,
ncols=None,
max_ncols=4,
):
"""
パネル数から適切な行数と列数を決定する。
Parameters
----------
n_panels : int
作成するパネルの総数。
ncols : int or None
列数を手動指定する場合の値。
Noneの場合は自動決定。
max_ncols : int
自動決定時の最大列数。
Returns
-------
nrows : int
行数。
ncols : int
列数。
"""
if n_panels < 1:
raise ValueError("n_panels must be at least 1.")
if ncols is None:
if n_panels <= 3:
ncols = n_panels
else:
ncols = int(np.ceil(np.sqrt(n_panels)))
ncols = min(ncols, max_ncols)
else:
ncols = int(ncols)
if ncols < 1:
raise ValueError("ncols must be at least 1.")
ncols = min(ncols, n_panels)
nrows = int(np.ceil(n_panels / ncols))
return nrows, ncols
# ============================================================
# Main
# ============================================================
def main():
# --------------------------------------------------------
# Energy settings
# --------------------------------------------------------
energy_min = 0.5
energy_max = 12.0
energy_step = 0.5
target_energies = np.arange(
energy_min,
energy_max + 0.5 * energy_step,
energy_step,
)
# 各パネルで表示するtarget energy周辺の範囲
x_half_width = 0.05
# 各パネル内にFWHMを表示
show_fwhm_text = True
output_filename = (
"rmf_response_tiles_"
"heg_meg_abs123_resolve_energy_depend.png"
)
# --------------------------------------------------------
# Figure layout settings
# --------------------------------------------------------
# Noneの場合はパネル数から列数を自動決定
# 例：3を指定すると必ず3列になる
ncols_override = None
# 自動決定時の最大列数
max_ncols = 4
# 1パネル当たりのおおよそのサイズ [inch]
panel_width = 3.35
panel_height = 2.45
# --------------------------------------------------------
# Input RMF directories
# --------------------------------------------------------
chandra_rmf_directory = (
"ss433/chandra/"
"repro/merge_all_tg"
)
resolve_rmf_filename = (
"ss433/resolve/"
"xa300041010rsl_source_Hp_all.rmf"
)
abs_orders = [1, 2, 3]
colors = {
1: "C0",
2: "C1",
3: "C2",
}
# --------------------------------------------------------
# RMF definitions
# --------------------------------------------------------
series_definitions = []
for abs_order in abs_orders:
series_definitions.append(
{
"key": f"heg{abs_order}",
"label": f"HEG |m|={abs_order}",
"short_label": f"HEG{abs_order}",
"filename": os.path.join(
chandra_rmf_directory,
f"combo_heg_abs{abs_order}.rmf",
),
"ext": 1,
"color": colors[abs_order],
"linestyle": "-",
"linewidth": 1.5,
"alpha": 1.0,
}
)
series_definitions.append(
{
"key": f"meg{abs_order}",
"label": f"MEG |m|={abs_order}",
"short_label": f"MEG{abs_order}",
"filename": os.path.join(
chandra_rmf_directory,
f"combo_meg_abs{abs_order}.rmf",
),
"ext": 1,
"color": colors[abs_order],
"linestyle": ":",
"linewidth": 1.5,
"alpha": 0.75,
}
)
series_definitions.append(
{
"key": "resolve",
"label": "Resolve",
"short_label": "Resolve",
"filename": resolve_rmf_filename,
"ext": 0,
"color": "black",
"linestyle": "-",
"linewidth": 2.0,
"alpha": 1.0,
}
)
# --------------------------------------------------------
# Load RMFs
# --------------------------------------------------------
all_response_data = {}
for definition in series_definitions:
key = definition["key"]
all_response_data[key] = collect_rmf_responses(
rmf_filename=definition["filename"],
target_energies=target_energies,
ext=definition["ext"],
label=definition["label"],
)
# --------------------------------------------------------
# Automatically determine figure layout
# --------------------------------------------------------
n_panels = len(target_energies)
nrows, ncols = choose_panel_grid(
n_panels=n_panels,
ncols=ncols_override,
max_ncols=max_ncols,
)
# 凡例の列数も図の横幅に合わせて決定
legend_ncols = min(
len(series_definitions),
max(4, 2 * ncols - 1),
)
legend_nrows = int(
np.ceil(len(series_definitions) / legend_ncols)
)
# タイトル・凡例・軸ラベル分を追加
header_height = 0.85 + 0.25 * legend_nrows
footer_height = 0.35
figure_width = panel_width * ncols
figure_height = (
panel_height * nrows
+ header_height
+ footer_height
)
print()
print("=" * 72)
print("Figure layout")
print(f" Number of panels : {n_panels}")
print(f" Grid : {nrows} rows x {ncols} columns")
print(
f" Figure size : "
f"{figure_width:.2f} x {figure_height:.2f} inch"
)
print("=" * 72)
print()
# --------------------------------------------------------
# Figure
# --------------------------------------------------------
fig, axes = plt.subplots(
nrows=nrows,
ncols=ncols,
figsize=(figure_width, figure_height),
sharex=False,
sharey=False,
squeeze=False,
)
axes = axes.ravel()
# --------------------------------------------------------
# Plot each target energy
# --------------------------------------------------------
for panel_index, target_energy in enumerate(target_energies):
ax = axes[panel_index]
fwhm_text_lines = []
for definition in series_definitions:
key = definition["key"]
rmf_data = all_response_data[key]
response_info = rmf_data["responses"][panel_index]
# RMF範囲外
if response_info is None:
fwhm_text_lines.append(
f"{definition['short_label']:<7s}: N/A"
)
continue
density = response_info["density"]
fwhm = response_info["fwhm"]
edges = rmf_data["edges"]
ax.step(
edges[:-1],
density,
where="post",
color=definition["color"],
linestyle=definition["linestyle"],
linewidth=definition["linewidth"],
alpha=definition["alpha"],
)
if np.isfinite(fwhm):
fwhm_text_lines.append(
f"{definition['short_label']:<7s}: "
f"{fwhm * 1e3:5.1f} eV"
)
else:
fwhm_text_lines.append(
f"{definition['short_label']:<7s}: N/A"
)
# 入力エネルギー位置
ax.axvline(
target_energy,
color="0.5",
linestyle="--",
linewidth=0.8,
alpha=0.7,
zorder=0,
)
ax.set_xlim(
target_energy - x_half_width,
target_energy + x_half_width,
)
ax.set_ylim(bottom=0)
# タイトルとグラフの距離を縮める
ax.set_title(
f"{target_energy:.1f} keV",
fontsize=10.5,
pad=2,
)
ax.grid(
visible=True,
which="major",
alpha=0.35,
)
ax.tick_params(
which="both",
direction="in",
top=True,
right=True,
labelsize=8,
pad=2,
)
ax.minorticks_on()
# ----------------------------------------------------
# FWHM values
# ----------------------------------------------------
if show_fwhm_text:
ax.text(
0.025,
0.975,
"\n".join(fwhm_text_lines),
transform=ax.transAxes,
ha="left",
va="top",
fontsize=5.7,
family="monospace",
linespacing=1.05,
bbox={
"facecolor": "white",
"edgecolor": "0.7",
"alpha": 0.78,
"boxstyle": "round,pad=0.20",
},
)
# --------------------------------------------------------
# Hide unused panels
# --------------------------------------------------------
for panel_index in range(
len(target_energies),
len(axes),
):
axes[panel_index].set_visible(False)
# --------------------------------------------------------
# Global legend
# --------------------------------------------------------
legend_handles = []
for definition in series_definitions:
legend_handles.append(
Line2D(
[0],
[0],
color=definition["color"],
linestyle=definition["linestyle"],
linewidth=definition["linewidth"],
alpha=definition["alpha"],
label=definition["label"],
)
)
# --------------------------------------------------------
# Global title
# --------------------------------------------------------
fig.suptitle(
"RMF Responses of Chandra HETG and XRISM Resolve\n"
f"Input energies: {energy_min:.1f}–{energy_max:.1f} keV "
f"in {energy_step:.1f} keV steps",
fontsize=14,
y=0.995,
linespacing=1.0,
)
# 凡例をタイトルの下、パネルの上に配置
fig.legend(
handles=legend_handles,
loc="upper center",
bbox_to_anchor=(0.5, 0.955),
ncol=legend_ncols,
fontsize=8.5,
frameon=True,
handlelength=2.5,
columnspacing=1.2,
handletextpad=0.5,
borderpad=0.4,
)
# --------------------------------------------------------
# Global axis labels
# --------------------------------------------------------
fig.supxlabel(
"Output Energy (keV)",
fontsize=12.5,
y=0.012,
)
fig.supylabel(
r"Probability Density (keV$^{-1}$)",
fontsize=12.5,
x=0.008,
)
# --------------------------------------------------------
# Compact spacing
# --------------------------------------------------------
axes_top = 0.915 - 0.025 * (legend_nrows - 1)
fig.subplots_adjust(
left=0.065,
right=0.992,
bottom=0.050,
top=axes_top,
wspace=0.17,
hspace=0.23,
)
fig.savefig(
output_filename,
dpi=200,
bbox_inches="tight",
pad_inches=0.04,
)
print(f"Saved figure: {output_filename}")
plt.show()
if __name__ == "__main__":
main()
タイル表示にすると、装置ごとのエネルギー依存性が一目で分かります。
Resolveでは広いエネルギー帯域で数eV程度の幅が維持されます。一方HETGでは、高エネルギーほどFWHMが増加し、高次光ほど主ピークが細くなることが確認できます。
入力エネルギーとFWHMの関係
最後に、各入力エネルギーについて求めたFWHMをプロットしてみます。
コードはこちら
import csv
import os
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import MultipleLocator
from astropy.io import fits
# ============================================================
# Matplotlib settings
# ============================================================
plt.rcParams["xtick.direction"] = "in"
plt.rcParams["ytick.direction"] = "in"
plt.rcParams["font.size"] = 12
# ============================================================
# RMF matrix construction
# This function is intentionally unchanged.
# ============================================================
def construct_rmf_matrix(filename, ext=0):
"""
RMFからレスポンスマトリックスを構築する関数（Chandra用）。
Parameters:
filename (str): FITSファイルのパス。
ext: indexの始まり
Returns:
dict: 構築されたレスポンスマトリックスと関連データを含む辞書。
"""
with fits.open(filename) as f:
# 必要なデータを抽出
f_chan = f[1].data['F_CHAN']-ext # チャンネルの開始位置（Chandraは1始まり、XIRMは0始まりなので注意）
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
# ============================================================
# Input-energy bin selection
# ============================================================
def find_input_energy_bin(in_ene_lo, in_ene_hi, target_energy):
"""Return the RMF input-bin index containing target_energy."""
indices = np.where(
(in_ene_lo <= target_energy)
& (target_energy < in_ene_hi)
)[0]
# Include the final upper boundary when target_energy is exactly equal to it.
if len(indices) == 0 and np.isclose(target_energy, in_ene_hi[-1]):
return len(in_ene_hi) - 1
if len(indices) == 0:
return None
return int(indices[0])
# ============================================================
# FWHM calculation
# ============================================================
def calc_fwhm(energy, response_density):
"""
Calculate the FWHM of the main response peak using linear interpolation.
Only the contiguous half-maximum region containing the global maximum is
used. This prevents a distant escape peak or secondary response component
from artificially broadening the measured FWHM.
Parameters
----------
energy : numpy.ndarray
Output-energy bin centers [keV].
response_density : numpy.ndarray
Response probability density [keV^-1].
Returns
-------
float
FWHM [keV]. Returns np.nan when it cannot be determined.
"""
energy = np.asarray(energy, dtype=float)
response_density = np.asarray(response_density, dtype=float)
valid = np.isfinite(energy) & np.isfinite(response_density)
energy = energy[valid]
response_density = response_density[valid]
if len(energy) < 3:
return np.nan
peak_index = int(np.argmax(response_density))
peak_value = response_density[peak_index]
if not np.isfinite(peak_value) or peak_value <= 0:
return np.nan
half_maximum = 0.5 * peak_value
left_inside = peak_index
while (
left_inside > 0
and response_density[left_inside - 1] >= half_maximum
):
left_inside -= 1
right_inside = peak_index
while (
right_inside < len(response_density) - 1
and response_density[right_inside + 1] >= half_maximum
):
right_inside += 1
if left_inside == 0 or right_inside == len(response_density) - 1:
return np.nan
def interpolate_crossing(i1, i2):
x1, x2 = energy[i1], energy[i2]
y1, y2 = response_density[i1], response_density[i2]
if np.isclose(y1, y2):
return 0.5 * (x1 + x2)
return x1 + (
(half_maximum - y1)
* (x2 - x1)
/ (y2 - y1)
)
left_crossing = interpolate_crossing(left_inside - 1, left_inside)
right_crossing = interpolate_crossing(right_inside, right_inside + 1)
fwhm = right_crossing - left_crossing
return fwhm if fwhm > 0 else np.nan
# ============================================================
# Shared plotting definitions
# ============================================================
def build_series_definitions(chandra_rmf_directory, resolve_rmf_filename):
"""Build consistent labels and plotting styles for all three scripts."""
colors = {1: "C0", 2: "C1", 3: "C2"}
definitions = []
for abs_order in (1, 2, 3):
definitions.append(
{
"key": f"heg{abs_order}",
"label": f"HEG |m|={abs_order}",
"short_label": f"HEG{abs_order}",
"csv_label": f"HEG_m{abs_order}_FWHM_eV",
"filename": os.path.join(
chandra_rmf_directory,
f"combo_heg_abs{abs_order}.rmf",
),
"ext": 1,
"instrument": "HEG",
"order": abs_order,
"color": colors[abs_order],
"linestyle": "-",
"marker": "o",
"linewidth": 1.7,
"alpha": 1.0,
}
)
definitions.append(
{
"key": f"meg{abs_order}",
"label": f"MEG |m|={abs_order}",
"short_label": f"MEG{abs_order}",
"csv_label": f"MEG_m{abs_order}_FWHM_eV",
"filename": os.path.join(
chandra_rmf_directory,
f"combo_meg_abs{abs_order}.rmf",
),
"ext": 1,
"instrument": "MEG",
"order": abs_order,
"color": colors[abs_order],
"linestyle": ":",
"marker": "s",
"linewidth": 1.7,
"alpha": 0.8,
}
)
definitions.append(
{
"key": "resolve",
"label": "XRISM Resolve",
"short_label": "Resolve",
"csv_label": "Resolve_FWHM_eV",
"filename": resolve_rmf_filename,
"ext": 0,
"instrument": "Resolve",
"order": None,
"color": "black",
"linestyle": "-",
"marker": "D",
"linewidth": 2.3,
"alpha": 1.0,
}
)
return definitions
# ============================================================
# Calculate an FWHM curve for one RMF
# ============================================================
def calculate_fwhm_curve(rmf_filename, target_energies, ext, label):
"""Calculate the main-peak FWHM at each requested input energy."""
if not os.path.isfile(rmf_filename):
raise FileNotFoundError(f"{label}: RMF file not found:\n{rmf_filename}")
print("=" * 72)
print(label)
print(rmf_filename)
matrix, in_lo, in_hi, out_lo, out_hi = construct_rmf_matrix(
rmf_filename,
ext=ext,
)
out_centers = 0.5 * (out_lo + out_hi)
out_widths = out_hi - out_lo
if np.any(out_widths <= 0):
raise ValueError(f"{label}: non-positive output-energy bin width found.")
fwhm_values_ev = np.full(len(target_energies), np.nan, dtype=float)
input_centers = np.full(len(target_energies), np.nan, dtype=float)
for i, target_energy in enumerate(target_energies):
input_index = find_input_energy_bin(in_lo, in_hi, target_energy)
if input_index is None:
print(f" {target_energy:5.2f} keV : outside RMF input range")
continue
response = np.asarray(matrix[input_index], dtype=float)
response_density = response / out_widths
fwhm_kev = calc_fwhm(out_centers, response_density)
input_centers[i] = 0.5 * (in_lo[input_index] + in_hi[input_index])
if np.isfinite(fwhm_kev):
fwhm_values_ev[i] = fwhm_kev * 1.0e3
print(
f" {target_energy:5.2f} keV : "
f"FWHM={fwhm_values_ev[i]:8.2f} eV, "
f"input-bin center={input_centers[i]:.6f} keV"
)
else:
print(f" {target_energy:5.2f} keV : FWHM=N/A")
print()
return fwhm_values_ev, input_centers
# ============================================================
# CSV output
# ============================================================
def save_results_csv(filename, target_energies, definitions, results, input_centers):
"""Save target energies, selected RMF-bin centers, and FWHM values."""
with open(filename, "w", newline="") as file_obj:
writer = csv.writer(file_obj)
header = ["Target_Energy_keV"]
for definition in definitions:
header.extend(
[
f'{definition["key"]}_InputBinCenter_keV',
definition["csv_label"],
]
)
writer.writerow(header)
for i, target_energy in enumerate(target_energies):
row = [f"{target_energy:.3f}"]
for definition in definitions:
key = definition["key"]
center = input_centers[key][i]
fwhm_value = results[key][i]
row.append(f"{center:.6f}" if np.isfinite(center) else "")
row.append(f"{fwhm_value:.6f}" if np.isfinite(fwhm_value) else "")
writer.writerow(row)
print(f"Saved CSV: {filename}")
# ============================================================
# Main
# ============================================================
def main():
energy_min = 0.5
energy_max = 12.0
energy_step = 0.5
target_energies = np.arange(
energy_min,
energy_max + 0.5 * energy_step,
energy_step,
)
# Use "linear" or "log". The scale is applied before saving the figure.
y_scale = "linear"
output_figure = f"rmf_fwhm_vs_energy_heg_meg_abs123_resolve_{y_scale}.png"
output_csv = "rmf_fwhm_vs_energy_heg_meg_abs123_resolve.csv"
chandra_rmf_directory = (
"ss433/chandra/"
"repro/merge_all_tg"
)
resolve_rmf_filename = (
"ss433/resolve/"
"xa300041010rsl_source_Hp_all.rmf"
)
definitions = build_series_definitions(
chandra_rmf_directory,
resolve_rmf_filename,
)
results = {}
input_centers = {}
for definition in definitions:
values, centers = calculate_fwhm_curve(
rmf_filename=definition["filename"],
target_energies=target_energies,
ext=definition["ext"],
label=definition["label"],
)
results[definition["key"]] = values
input_centers[definition["key"]] = centers
save_results_csv(
filename=output_csv,
target_energies=target_energies,
definitions=definitions,
results=results,
input_centers=input_centers,
)
fig, ax = plt.subplots(figsize=(8.2, 5.6))
for definition in definitions:
ax.plot(
target_energies,
results[definition["key"]],
color=definition["color"],
linestyle=definition["linestyle"],
linewidth=definition["linewidth"],
alpha=definition["alpha"],
marker=definition["marker"],
markersize=4,
label=definition["label"],
zorder=10 if definition["instrument"] == "Resolve" else None,
)
ax.set_xlim(energy_min - 0.1 * energy_step, energy_max + 0.1 * energy_step)
ax.xaxis.set_major_locator(MultipleLocator(1))
ax.set_yscale(y_scale)
ax.set_xlabel("Input Energy (keV)")
ax.set_ylabel("Main-peak FWHM (eV)")
ax.set_title("Energy resolution: Chandra HETG vs XRISM Resolve", pad=5)
ax.tick_params(which="both", top=True, right=True, pad=3)
ax.minorticks_on()
ax.grid(visible=True, which="major", alpha=0.4)
ax.grid(visible=True, which="minor", alpha=0.15)
ax.legend(
fontsize="small",
ncol=1,
frameon=True,
handlelength=2.6,
)
fig.tight_layout(pad=0.7)
fig.savefig(output_figure, dpi=200, bbox_inches="tight", pad_inches=0.04)
print(f"Saved figure: {output_figure}")
plt.show()
if __name__ == "__main__":
main()
-
縦軸リニア
- ResolveのFWHMあたりの拡大図
- ResolveのFWHMあたりの拡大図
-
縦軸ログ
この図から、ResolveとHETGではエネルギー分解能のエネルギー依存性が大きく異なることが分かります。
| 傾向 | 主な理由 |
|---|---|
| ResolveのFWHMはほぼ一定だが、高エネルギー側でわずかに増加する | ほぼ一定のbaseline resolutionに、エネルギーに比例するexcess broadeningが加わるため |
| HEGはMEGよりFWHMが小さい | 格子周期が短く、波長分散が大きいため |
| 高次光ほどFWHMが小さい | 回折次数が高いほど波長分散が大きいため |
| HETGは高エネルギーほどFWHMが増加する | 波長幅をエネルギー幅へ変換すると、概ねエネルギーの二乗に比例するため |
補足：Resolving power で比較する
エネルギー分解能は、FWHMだけでなく
R=\frac{E}{\Delta E}
というresolving powerで表すこともできます。
波長で書けば
R=\frac{\lambda}{\Delta\lambda}
であり、これは光学分光などでもよく使われる指標です。そのため、多波長でline profileを比較するときには、装置ごとの分解能を同じ無次元量として整理できます。
ここでは、RMFから求めた主ピークのFWHMを$\Delta E$として、
R(E)=\frac{E}{\mathrm{FWHM}(E)}
を計算してみます。
コードはこちら
import csv
import os
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import MultipleLocator
from astropy.io import fits
# ============================================================
# Matplotlib settings
# ============================================================
plt.rcParams["xtick.direction"] = "in"
plt.rcParams["ytick.direction"] = "in"
plt.rcParams["font.size"] = 12
# ============================================================
# RMF matrix construction
# ============================================================
def construct_rmf_matrix(filename, ext=0):
"""
RMFからレスポンスマトリックスを構築する関数。
Parameters
----------
filename : str
FITS RMF filename.
ext : int
F_CHAN offset.
Chandra: ext=1
XRISM: ext=0
Returns
-------
matrix
in_ene_lo
in_ene_hi
out_ene_lo
out_ene_hi
"""
with fits.open(filename) as f:
f_chan = f[1].data["F_CHAN"] - ext
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=(
num_in_ene_bins,
num_out_ene_bins,
),
dtype="float32",
)
for j in range(
f_chan.shape[0]
):
tot_n = 0
for f_start, n in zip(
f_chan[j],
n_chan[j],
):
matrix[
j,
f_start:f_start + 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,
)
# ============================================================
# Input-energy bin selection
# ============================================================
def find_input_energy_bin(
in_ene_lo,
in_ene_hi,
target_energy,
):
"""
Return the RMF input-bin index containing target_energy.
"""
indices = np.where(
(in_ene_lo <= target_energy)
& (target_energy < in_ene_hi)
)[0]
# Include final upper boundary.
if (
len(indices) == 0
and np.isclose(
target_energy,
in_ene_hi[-1],
)
):
return len(
in_ene_hi
) - 1
if len(
indices
) == 0:
return None
return int(
indices[0]
)
# ============================================================
# FWHM calculation
# ============================================================
def calc_fwhm(
energy,
response_density,
):
"""
Calculate the FWHM of the main response peak
using linear interpolation.
Only the contiguous half-maximum region containing
the global maximum is used.
Parameters
----------
energy : numpy.ndarray
Output-energy bin centers [keV].
response_density : numpy.ndarray
Response probability density [keV^-1].
Returns
-------
float
FWHM [keV].
"""
energy = np.asarray(
energy,
dtype=float,
)
response_density = np.asarray(
response_density,
dtype=float,
)
valid = (
np.isfinite(
energy
)
& np.isfinite(
response_density
)
)
energy = energy[
valid
]
response_density = (
response_density[
valid
]
)
if len(
energy
) < 3:
return np.nan
peak_index = int(
np.argmax(
response_density
)
)
peak_value = (
response_density[
peak_index
]
)
if (
not np.isfinite(
peak_value
)
or peak_value <= 0
):
return np.nan
half_maximum = (
0.5
* peak_value
)
# ========================================================
# Search left
# ========================================================
left_inside = (
peak_index
)
while (
left_inside > 0
and response_density[
left_inside - 1
] >= half_maximum
):
left_inside -= 1
# ========================================================
# Search right
# ========================================================
right_inside = (
peak_index
)
while (
right_inside
< len(
response_density
) - 1
and response_density[
right_inside + 1
] >= half_maximum
):
right_inside += 1
if (
left_inside == 0
or right_inside
== len(
response_density
) - 1
):
return np.nan
# ========================================================
# Linear interpolation at half maximum
# ========================================================
def interpolate_crossing(
i1,
i2,
):
x1 = energy[
i1
]
x2 = energy[
i2
]
y1 = response_density[
i1
]
y2 = response_density[
i2
]
if np.isclose(
y1,
y2,
):
return (
0.5
* (
x1
+ x2
)
)
return (
x1
+ (
(
half_maximum
- y1
)
* (
x2
- x1
)
/ (
y2
- y1
)
)
)
left_crossing = (
interpolate_crossing(
left_inside - 1,
left_inside,
)
)
right_crossing = (
interpolate_crossing(
right_inside,
right_inside + 1,
)
)
fwhm = (
right_crossing
- left_crossing
)
if fwhm <= 0:
return np.nan
return fwhm
# ============================================================
# Shared plotting definitions
# ============================================================
def build_series_definitions(
chandra_rmf_directory,
resolve_rmf_filename,
):
"""
Build plotting definitions for
HEG/MEG orders 1--3 and XRISM Resolve.
"""
colors = {
1: "C0",
2: "C1",
3: "C2",
}
definitions = []
for abs_order in (
1,
2,
3,
):
# ====================================================
# HEG
# ====================================================
definitions.append(
{
"key":
f"heg{abs_order}",
"label":
f"HEG |m|={abs_order}",
"short_label":
f"HEG{abs_order}",
"filename":
os.path.join(
chandra_rmf_directory,
f"combo_heg_abs{abs_order}.rmf",
),
"ext":
1,
"instrument":
"HEG",
"order":
abs_order,
"color":
colors[
abs_order
],
"linestyle":
"-",
"marker":
"o",
"linewidth":
1.7,
"alpha":
1.0,
}
)
# ====================================================
# MEG
# ====================================================
definitions.append(
{
"key":
f"meg{abs_order}",
"label":
f"MEG |m|={abs_order}",
"short_label":
f"MEG{abs_order}",
"filename":
os.path.join(
chandra_rmf_directory,
f"combo_meg_abs{abs_order}.rmf",
),
"ext":
1,
"instrument":
"MEG",
"order":
abs_order,
"color":
colors[
abs_order
],
"linestyle":
":",
"marker":
"s",
"linewidth":
1.7,
"alpha":
0.8,
}
)
# ========================================================
# XRISM Resolve
# ========================================================
definitions.append(
{
"key":
"resolve",
"label":
"XRISM Resolve",
"short_label":
"Resolve",
"filename":
resolve_rmf_filename,
"ext":
0,
"instrument":
"Resolve",
"order":
None,
"color":
"black",
"linestyle":
"-",
"marker":
"D",
"linewidth":
2.3,
"alpha":
1.0,
}
)
return definitions
# ============================================================
# Calculate FWHM and resolving-power curves
# ============================================================
def calculate_resolution_curve(
rmf_filename,
target_energies,
ext,
label,
):
"""
Calculate:
FWHM [eV]
and
R = E / DeltaE_FWHM
at each requested input energy.
E is taken from the actual selected RMF input-bin center.
"""
if not os.path.isfile(
rmf_filename
):
raise FileNotFoundError(
f"{label}: "
f"RMF file not found:\n"
f"{rmf_filename}"
)
print(
"=" * 72
)
print(
label
)
print(
rmf_filename
)
(
matrix,
in_lo,
in_hi,
out_lo,
out_hi,
) = construct_rmf_matrix(
rmf_filename,
ext=ext,
)
# ========================================================
# Output-energy grid
# ========================================================
out_centers = (
0.5
* (
out_lo
+ out_hi
)
)
out_widths = (
out_hi
- out_lo
)
if np.any(
out_widths <= 0
):
raise ValueError(
f"{label}: "
"non-positive output-energy bin width found."
)
# ========================================================
# Output arrays
# ========================================================
fwhm_values_ev = np.full(
len(
target_energies
),
np.nan,
dtype=float,
)
resolving_power = np.full(
len(
target_energies
),
np.nan,
dtype=float,
)
input_centers = np.full(
len(
target_energies
),
np.nan,
dtype=float,
)
# ========================================================
# Energy loop
# ========================================================
for i, target_energy in enumerate(
target_energies
):
input_index = (
find_input_energy_bin(
in_lo,
in_hi,
target_energy,
)
)
if input_index is None:
print(
f" {target_energy:5.2f} keV : "
"outside RMF input range"
)
continue
# ====================================================
# RMF response
# ====================================================
response = np.asarray(
matrix[
input_index
],
dtype=float,
)
# Convert probability / bin
# to probability density / keV.
response_density = (
response
/ out_widths
)
# ====================================================
# FWHM
# ====================================================
fwhm_kev = calc_fwhm(
out_centers,
response_density,
)
# ====================================================
# Actual RMF input-bin center
# ====================================================
input_center_kev = (
0.5
* (
in_lo[
input_index
]
+ in_hi[
input_index
]
)
)
input_centers[
i
] = input_center_kev
# ====================================================
# Resolving power
#
# R = E / DeltaE
#
# both quantities must have same unit.
# ====================================================
if np.isfinite(
fwhm_kev
):
fwhm_ev = (
fwhm_kev
* 1.0e3
)
fwhm_values_ev[
i
] = fwhm_ev
resolving_power[
i
] = (
input_center_kev
/ fwhm_kev
)
print(
f" {target_energy:5.2f} keV : "
f"FWHM = "
f"{fwhm_ev:8.3f} eV, "
f"R = "
f"{resolving_power[i]:8.1f}, "
f"input-bin center = "
f"{input_center_kev:.6f} keV"
)
else:
print(
f" {target_energy:5.2f} keV : "
"FWHM=N/A, R=N/A"
)
print()
return (
fwhm_values_ev,
resolving_power,
input_centers,
)
# ============================================================
# CSV output
# ============================================================
def save_results_csv(
filename,
target_energies,
definitions,
fwhm_results,
resolving_power_results,
input_centers,
):
"""
Save target energies, actual RMF-bin centers,
FWHM, and resolving power.
"""
with open(
filename,
"w",
newline="",
) as file_obj:
writer = csv.writer(
file_obj
)
# ====================================================
# Header
# ====================================================
header = [
"Target_Energy_keV"
]
for definition in definitions:
key = definition[
"key"
]
header.extend(
[
f"{key}_InputBinCenter_keV",
f"{key}_FWHM_eV",
f"{key}_ResolvingPower_R",
]
)
writer.writerow(
header
)
# ====================================================
# Rows
# ====================================================
for i, target_energy in enumerate(
target_energies
):
row = [
f"{target_energy:.3f}"
]
for definition in definitions:
key = definition[
"key"
]
center = (
input_centers[
key
][
i
]
)
fwhm_value = (
fwhm_results[
key
][
i
]
)
resolving_power = (
resolving_power_results[
key
][
i
]
)
row.append(
f"{center:.6f}"
if np.isfinite(
center
)
else ""
)
row.append(
f"{fwhm_value:.6f}"
if np.isfinite(
fwhm_value
)
else ""
)
row.append(
f"{resolving_power:.6f}"
if np.isfinite(
resolving_power
)
else ""
)
writer.writerow(
row
)
print(
f"Saved CSV: {filename}"
)
# ============================================================
# Main
# ============================================================
def main():
# ========================================================
# Energy grid
# ========================================================
energy_min = 0.5
energy_max = 12.0
energy_step = 0.5
target_energies = np.arange(
energy_min,
energy_max
+ 0.5
* energy_step,
energy_step,
)
# ========================================================
# Y-axis scale
#
# "linear"
# "log"
# ========================================================
y_scale = "linear"
# ========================================================
# Output
# ========================================================
output_figure = (
"rmf_resolving_power_vs_energy_"
"heg_meg_abs123_resolve_"
f"{y_scale}.png"
)
output_csv = (
"rmf_resolving_power_vs_energy_"
"heg_meg_abs123_resolve.csv"
)
# ========================================================
# RMF paths
# ========================================================
chandra_rmf_directory = (
"ss433/chandra/"
"repro/merge_all_tg"
)
resolve_rmf_filename = (
"ss433/resolve/"
"xa300041010rsl_source_Hp_all.rmf"
)
# ========================================================
# Definitions
# ========================================================
definitions = (
build_series_definitions(
chandra_rmf_directory,
resolve_rmf_filename,
)
)
# ========================================================
# Calculate curves
# ========================================================
fwhm_results = {}
resolving_power_results = {}
input_centers = {}
for definition in definitions:
(
fwhm_values,
resolving_power,
centers,
) = calculate_resolution_curve(
rmf_filename=definition[
"filename"
],
target_energies=target_energies,
ext=definition[
"ext"
],
label=definition[
"label"
],
)
key = definition[
"key"
]
fwhm_results[
key
] = fwhm_values
resolving_power_results[
key
] = resolving_power
input_centers[
key
] = centers
# ========================================================
# Save CSV
# ========================================================
save_results_csv(
filename=output_csv,
target_energies=target_energies,
definitions=definitions,
fwhm_results=fwhm_results,
resolving_power_results=(
resolving_power_results
),
input_centers=input_centers,
)
# ========================================================
# Plot
# ========================================================
fig, ax = plt.subplots(
figsize=(
8.2,
5.6,
)
)
for definition in definitions:
key = definition[
"key"
]
ax.plot(
target_energies,
resolving_power_results[
key
],
color=definition[
"color"
],
linestyle=definition[
"linestyle"
],
linewidth=definition[
"linewidth"
],
alpha=definition[
"alpha"
],
marker=definition[
"marker"
],
markersize=4,
label=definition[
"label"
],
zorder=(
10
if definition[
"instrument"
] == "Resolve"
else None
),
)
# ========================================================
# Axis
# ========================================================
ax.set_xlim(
energy_min
- 0.1
* energy_step,
energy_max
+ 0.1
* energy_step,
)
ax.xaxis.set_major_locator(
MultipleLocator(
1
)
)
ax.set_yscale(
y_scale
)
ax.set_xlabel(
"Input Energy (keV)"
)
ax.set_ylabel(
r"Resolving power $R=E/\Delta E_{\rm FWHM}$"
)
ax.set_title(
"Spectral resolving power: "
"Chandra HETG vs XRISM Resolve",
pad=5,
)
ax.tick_params(
which="both",
top=True,
right=True,
pad=3,
)
ax.minorticks_on()
# ========================================================
# Grid
# ========================================================
ax.grid(
visible=True,
which="major",
alpha=0.4,
)
ax.grid(
visible=True,
which="minor",
alpha=0.15,
)
# ========================================================
# Legend
# ========================================================
ax.legend(
fontsize="small",
ncol=1,
frameon=True,
handlelength=2.6,
)
# ========================================================
# Save
# ========================================================
fig.tight_layout(
pad=0.7
)
fig.savefig(
output_figure,
dpi=200,
bbox_inches="tight",
pad_inches=0.04,
)
print(
f"Saved figure: "
f"{output_figure}"
)
plt.show()
# ============================================================
# Execute
# ============================================================
if __name__ == "__main__":
main()
ResolveはFWHMが数eV程度でほぼ一定なので、高エネルギーほど$R$が大きくなります。一方、HETGではエネルギー依存性が異なり、高次光ほど高い$R$を持ちます。
また、$R$はドップラー速度のスケールに換算すると、
\Delta v_{\rm FWHM}\simeq\frac{c}{R}
と表せます。
そのため、多波長で輝線profileを比較する際には、各装置がどの程度の速度幅や複数の輝線成分を分離できるかを比較する目安としても利用できます。
ただし、これは装置のFWHMに対応する分解能の指標です。十分な光子数が得られている場合には、一本の輝線の中心位置（centroid）はFWHMよりも高い精度で決定できるため、$c/R$はドップラーシフトの測定精度そのものを表すわけではない点に注意が必要です。
補足：RMFだけでは観測性能を語れない
今回比較したのは、RMFから分かるline responseとそのFWHMです。
ただし、エネルギー分解能だけを比較しても、それぞれの装置の観測性能を十分に評価できない場合があります。高いエネルギー分解能を持っていても、得られる光子数が少なく統計が十分でなければ、期待するサイエンスにつなげることが難しくなります。
どれだけ効率よくX線光子を集めて検出できるかは、X線望遠鏡の結像光学系、回折格子、検出器などの装置設計や、観測時の条件によって決まります。
スペクトル解析では、これらの集光効率や検出効率をエネルギーごとの仮想的な面積に換算し、有効面積として表します。その情報を格納したレスポンスファイルがAuxiliary Response File（ARF）です。
したがって、装置の分光性能を比較する際には、RMFから分かるエネルギー分解能だけでなく、装置全体の光子収集能力を反映したARFも併せて確認することが重要です。
以下に、今回比較した観測データ（ObsID）に対して作成されたARFを示します。
有効面積とARFについて
ARFに記録された有効面積は、装置そのものの幾何学的な面積ではありません。望遠鏡の集光効率や、回折格子、検出器などの効率を、エネルギーごとの仮想的な面積として集約したものです。また、ARFは観測位置や抽出領域などにも依存するため、ここで示すのは今回使用した観測データに対する一例です。
コードはこちら
import os
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import MultipleLocator
from astropy.io import fits
# ============================================================
# Matplotlib settings
# ============================================================
plt.rcParams["xtick.direction"] = "in"
plt.rcParams["ytick.direction"] = "in"
plt.rcParams["font.size"] = 12
# ============================================================
# ARF reader
# ============================================================
def load_arf(filename):
"""
ARFからエネルギービンと有効面積を読み込む。
Parameters
----------
filename : str
ARFファイルのパス。
Returns
-------
energy_lo : numpy.ndarray
エネルギービン下端 [keV]。
energy_hi : numpy.ndarray
エネルギービン上端 [keV]。
effective_area : numpy.ndarray
有効面積 [cm^2]。
"""
if not os.path.isfile(filename):
raise FileNotFoundError(
f"ARF file not found:\n{filename}"
)
with fits.open(filename) as hdul:
arf_hdu = None
for hdu in hdul:
if hdu.data is None:
continue
if not hasattr(hdu, "columns"):
continue
column_names = {
name.upper()
for name in hdu.columns.names
}
required_columns = {
"ENERG_LO",
"ENERG_HI",
"SPECRESP",
}
if required_columns.issubset(column_names):
arf_hdu = hdu
break
if arf_hdu is None:
raise ValueError(
"Could not find an ARF extension containing "
f"ENERG_LO, ENERG_HI, and SPECRESP:\n{filename}"
)
energy_lo = np.asarray(
arf_hdu.data["ENERG_LO"],
dtype=float,
)
energy_hi = np.asarray(
arf_hdu.data["ENERG_HI"],
dtype=float,
)
effective_area = np.asarray(
arf_hdu.data["SPECRESP"],
dtype=float,
)
valid = (
np.isfinite(energy_lo)
& np.isfinite(energy_hi)
& np.isfinite(effective_area)
& (energy_hi > energy_lo)
& (effective_area >= 0)
)
energy_lo = energy_lo[valid]
energy_hi = energy_hi[valid]
effective_area = effective_area[valid]
order = np.argsort(energy_lo)
return (
energy_lo[order],
energy_hi[order],
effective_area[order],
)
# ============================================================
# Shared plotting definitions
# ============================================================
def build_series_definitions(
chandra_arf_directory,
resolve_arf_filename,
):
"""
RMF図と同じ色・線種でARF系列を定義する。
次数:
|m|=1 -> C0
|m|=2 -> C1
|m|=3 -> C2
装置:
HEG -> 実線
MEG -> 点線
Resolve -> 黒の太線
"""
colors = {
1: "C0",
2: "C1",
3: "C2",
}
definitions = []
for abs_order in (1, 2, 3):
definitions.append(
{
"key": f"heg{abs_order}",
"label": f"HEG |m|={abs_order}",
"filename": os.path.join(
chandra_arf_directory,
f"combo_heg_abs{abs_order}.arf",
),
"instrument": "HEG",
"order": abs_order,
"color": colors[abs_order],
"linestyle": "-",
"linewidth": 1.7,
"alpha": 1.0,
"zorder": 4,
}
)
definitions.append(
{
"key": f"meg{abs_order}",
"label": f"MEG |m|={abs_order}",
"filename": os.path.join(
chandra_arf_directory,
f"combo_meg_abs{abs_order}.arf",
),
"instrument": "MEG",
"order": abs_order,
"color": colors[abs_order],
"linestyle": ":",
"linewidth": 1.7,
"alpha": 0.8,
"zorder": 3,
}
)
definitions.append(
{
"key": "resolve",
"label": "XRISM Resolve",
"filename": resolve_arf_filename,
"instrument": "Resolve",
"order": None,
"color": "black",
"linestyle": "-",
"linewidth": 2.3,
"alpha": 1.0,
"zorder": 5,
}
)
return definitions
# ============================================================
# Main
# ============================================================
def main():
energy_min = 0.3
energy_max = 12.0
# Use "linear" or "log". The scale is applied before saving the figure.
y_scale = "linear"
output_filename = (
f"arf_heg_meg_abs123_resolve_{y_scale}.png"
)
chandra_arf_directory = (
"ss433/chandra/"
"repro/merge_all_tg"
)
resolve_arf_filename = (
"ss433/resolve/"
"xa300041010rsl_source_Hp_all.arf"
)
definitions = build_series_definitions(
chandra_arf_directory,
resolve_arf_filename,
)
fig, ax = plt.subplots(figsize=(7.4, 5.2))
for definition in definitions:
energy_lo, energy_hi, effective_area = load_arf(
definition["filename"]
)
mask = (
(energy_hi >= energy_min)
& (energy_lo <= energy_max)
)
energy_lo = energy_lo[mask]
energy_hi = energy_hi[mask]
effective_area = effective_area[mask]
if len(effective_area) == 0:
print(
f'WARNING: No ARF bins in selected range: '
f'{definition["label"]}'
)
continue
# ARFはエネルギービンごとの値なのでstepで描画
ax.step(
energy_lo,
effective_area,
where="post",
color=definition["color"],
linestyle=definition["linestyle"],
linewidth=definition["linewidth"],
alpha=definition["alpha"],
zorder=definition["zorder"],
label=definition["label"],
)
print(f'[{definition["label"]}]')
print(
f" energy range : "
f"{energy_lo.min():.3f}–{energy_hi.max():.3f} keV"
)
print(
f" maximum area : "
f"{np.nanmax(effective_area):.3f} cm^2"
)
print()
ax.set_xlim(energy_min, energy_max)
ax.xaxis.set_major_locator(MultipleLocator(1))
if y_scale == "linear":
y_min = 0
elif y_scale == "log":
y_min = 0.1
ax.set_ylim(bottom=y_min)
ax.set_yscale(y_scale)
ax.set_xlabel("Energy (keV)")
ax.set_ylabel(r"Effective Area (cm$^2$)")
ax.set_title(
"ARF: Chandra HETG vs XRISM Resolve",
pad=5,
)
ax.grid(
visible=True,
which="major",
alpha=0.4,
)
ax.tick_params(
which="both",
top=True,
right=True,
pad=3,
)
ax.minorticks_on()
ax.legend(
fontsize="x-small",
ncol=1,
frameon=True,
handlelength=2.8,
)
fig.tight_layout(pad=0.7)
fig.savefig(
output_filename,
dpi=200,
bbox_inches="tight",
pad_inches=0.04,
)
print(f"Saved figure: {output_filename}")
plt.show()
if __name__ == "__main__":
main()
-
縦軸リニア
-
縦軸ログ
結果を見ると、Resolveは広いエネルギー帯域で高いエネルギー分解能を持つだけでなく、HETGの各回折次数と比べて大きな有効面積を持っています。そのため、精密な輝線構造を分離しながら、比較的多くの光子を集め、高い統計精度を得やすいことがResolveの大きな科学的強みだと思います。
まとめ
この記事では、実際の観測データから作成したRMFを用いて、XRISM/ResolveとChandra/HETGのLSFを比較しました。
RMFを直接可視化することで、Resolveでは広いエネルギー帯域で高いエネルギー分解能が維持される一方、HETGではHEG・MEGや回折次数、入力エネルギーによってエネルギー分解能が大きく変化することを確認しました。特に高次光では、より高いエネルギー分解能が得られることを示しました。
一方、実際の観測性能はRMFだけでは決まりません。ARFを比較すると、Resolveは高いエネルギー分解能に加え、HETGよりもはるかに大きな有効面積を持っていることが分かります。このように、分光性能を評価する際には、RMFによるエネルギー分解能だけでなく、ARFによる有効面積も併せて考えることが重要です。
普段はスペクトルフィッティングの中で意識することの少ないRMFやARFですが、その中身を可視化することで、各装置の特性やスペクトル解析で何が行われているのかをより直感的に理解できるでしょう。本記事が、X線精密分光で用いられるレスポンスファイルへの理解を深めるきっかけになれば幸いです。
参考記事
M. A. Leutenegger et al., “Core line spread function calibration of the X-ray Imaging and Spectroscopy Mission Resolve X-ray calorimeter spectrometer,” Journal of Astronomical Telescopes, Instruments, and Systems 11, 042024 (2025).
https://doi.org/10.1117/1.JATIS.11.4.042024
関連記事
RMFとARFの基本的な意味や可視化方法については、以下の記事で紹介していますので、よろしければご覧ください。