Chandraのバンド分割をもっと楽に：Pythonで作る「統計が揃う」エネルギーバンド

はじめに

Chandra の画像解析では、エネルギー依存性を見るためにバンド分割を行う場面が多くあります。
しかし、CIAO の一般的な soft / medium / hard のような固定バンドでは、

• 天体によってカウントが大きく偏る
• PSF 補正やデコンボリューションで統計誤差が揃わない
• バンド間で構造を比較しづらい

といった問題が生じることがあります。

特にイメージデコンボリューションでは、各バンドの統計量（count）を揃えることでノイズの影響を揃えられるため、
バンド間でカウントが大きく偏っているかどうかは実務上とても重要になります。

そこで本記事では、CIAO の標準機能では柔軟に扱いにくい
「統計がほぼ均等になるエネルギーバンド」を Python で自動生成する方法 を紹介します。

なお、この記事のコードは SS 433（ObsID 659）のイベントファイルで
すぐに動作確認できるようにサンプルデータも用意しています。

分割の考え方

今回の方法では、エネルギー全域を最小ビンから厳密に探索するわけではありません。

実務では、エネルギー境界が細かすぎても扱いにくいため、まず

• 適度な刻み幅 dE (例：0.1 keV) を決める
  → バンド境界をこの刻み上で探す

という前提を置きます。

そのうえで、

• 累積カウントが 1/n, 2/n, … に近い位置
  → その dE 刻みの中で 最も近い場所 を採用する

という 軽量なヒューリスティックにしています。

ここで強調したいのは、

• 全ての境界の組み合わせを調べて最適な分割を求めるような
  組合せ最適化は行っていない
• 標準偏差を最小化する“完全な均等割り” を求めているわけでもない

という点です。

あくまで「実務で使いやすく」「十分に均等な」エネルギーバンドを、
シンプルな計算で得られる方法として実装しています。

サンプルデータ（SS 433 / ObsID 659）

この記事の動作確認用に、軽量化した SS 433 のイベントファイル一式を
Google Drive に置いてあります。コードの確認にご活用ください。

ファイル
https://drive.google.com/file/d/1j5OqOt6hiBPPhU5naYsstzdbbMJGNzwh/view?usp=sharing

内容物

• acisf00659_repro_evt2.fits（ObsID 659 の evt2）
• region.reg（CIAO 用の領域ファイル）
• ds9_region.reg（DS9 表示用／スクリプトでは使用しません）

領域のイメージ（DS9表示例）：

この緑の領域内に含まれるイベントを取り出し、統計がほぼ等しくなるようにエネルギーバンドへ分割するのが本記事の目的です。

実装（pycrates）

実装には CIAO の pycrates （CIAO v 4.14で導入された）を使います。
region・energy・time を CIAO のフィルタ構文そのままで Python コード中に指定できるため、イベントファイルの抽出・加工を 柔軟かつ一貫した方法で扱えます。

今回の関数の仕様は次のとおりです：

• energy の範囲 [minE, maxE] を指定
• nbands 個に分割して各バンドの統計を揃える
• 境界の探索幅 dE を指定
• count と mean energy を返す
• mode="manual" で境界を手動指定できる

コード全体は以下の通りです。

import numpy as np
from pycrates import read_file
import matplotlib.pyplot as plt

# ============================================================
# ★ パラメータ（ユーザーが調整すべき部分）
# ============================================================

# 解析するエネルギー範囲 [keV]
minE = 0.5
maxE = 7.0

# 分割数（equal-count の区間をいくつ作るか）
nbands = 3

# エネルギー境界を丸める幅 [keV]
dE = 0.1

# 対象となるイベントファイル・領域ファイル
evtfile = "acisf00659_repro_evt2.fits"   # ObsID を変更する場合はここを編集
regfile = "region.reg"                   # 使用する領域ファイル

# ============================================================
# ★ イベントファイル読み込み
# ============================================================

evt = read_file(
    f"{evtfile}[sky=region({regfile}),energy={minE*1e3}:{maxE*1e3}]"
)

# energy 列（eV → keV に変換）
energy = evt.get_column("energy").values / 1000.0


# ============================================================
# ★ メイン：エネルギーバンド分割関数
# ============================================================
def split_bands_with_stats(
    energy, minE, maxE,
    nbands=3, dE=0.1,
    mode="auto",
    custom_edges=None,
):
    """
    Energy band splitter.

    mode = "auto"   → 統計均等割りで自動決定
    mode = "manual" → custom_edges に従う
    """

    # --------------------------------------------------------
    # 0. 指定されたエネルギー範囲にフィルタリング
    # --------------------------------------------------------
    mask = (energy >= minE) & (energy <= maxE)
    energy = energy[mask]

    if len(energy) == 0:
        raise ValueError("指定したエネルギー範囲内にイベントがありません。")

    # ========================================================
    # ★ mode="manual" : ユーザー指定の edges を使う
    # ========================================================
    if mode == "manual":
        if custom_edges is None:
            raise ValueError("mode='manual' の場合 custom_edges を指定してください。")

        # minE < edges < maxE を保証
        edges = [minE] + [e for e in custom_edges if minE < e < maxE] + [maxE]
        edges = sorted(edges)

        # 統計計算
        band_means = []
        band_counts = []

        for i in range(len(edges) - 1):
            lo, hi = edges[i], edges[i+1]
            mask = (energy >= lo) & (energy < hi)
            c = np.sum(mask)
            band_counts.append(c)
            band_means.append(np.mean(energy[mask]) if c > 0 else np.nan)

        bands = [(edges[i], edges[i+1]) for i in range(len(edges)-1)]
        return bands, band_means, band_counts, edges

    # ========================================================
    # ★ mode="auto" : 統計から自動で決める
    # ========================================================
    if mode == "auto":

        # dE ヒストグラム
        bins = int((maxE - minE) / dE)
        bins = max(bins, 1)

        counts, edges_hist = np.histogram(energy, range=(minE, maxE), bins=bins)
        cum = np.cumsum(counts)
        total = cum[-1]

        targets = np.linspace(0, total, nbands + 1)[1:-1]
        idx = [np.searchsorted(cum, t) for t in targets]

        edges = [minE] + [edges_hist[i] for i in idx] + [maxE]

        # 幅ゼロバンドの補正
        for i in range(1, len(edges)):
            if edges[i] <= edges[i - 1]:
                edges[i] = edges[i - 1] + dE

        if edges[-1] < maxE:
            edges[-1] = maxE

        # 統計計算
        band_means = []
        band_counts = []

        for i in range(len(edges) - 1):
            lo, hi = edges[i], edges[i+1]
            mask = (energy >= lo) & (energy < hi)
            c = np.sum(mask)
            band_counts.append(c)
            band_means.append(np.mean(energy[mask]) if c > 0 else np.nan)

        bands = [(edges[i], edges[i+1]) for i in range(len(edges)-1)]
        return bands, band_means, band_counts, edges

    else:
        raise ValueError("mode は 'auto' か 'manual' を指定してください。")



# ============================================================
# ★ 関数実行
# ============================================================
bands, means, counts, edges = split_bands_with_stats(
    energy, minE, maxE, nbands=nbands, dE=dE
)

total_count = np.sum(counts)

# ============================================================
# ★ 結果を出力
# ============================================================
print("Bands:")
for (lo, hi), m, c in zip(bands, means, counts):
    pct = 100 * c / total_count
    print(f"{lo:.1f}–{hi:.1f} keV | mean={m:.2f} keV | count={c} ({pct:.1f}%)")


# ============================================================
# ★ ヒストグラム + バンド境界の可視化
# ============================================================
plt.figure(figsize=(8, 4))
plt.hist(energy, bins=200, histtype='step', color='black')
plt.xlabel("Energy [keV]")
plt.ylabel("Counts")
plt.title("Energy Distribution and Band Boundaries")

for i, e in enumerate(edges):
    # 最初（i=0）と最後（i=len(edges)-1）は実線
    if i == 0 or i == len(edges) - 1:
        plt.axvline(e, color='tab:red', linestyle='-', alpha=0.9, linewidth=1.5)
    else:
        plt.axvline(e, color='tab:red', linestyle='--', alpha=0.7)

info = "Bands:\n"
for (lo, hi), m, c in zip(bands, means, counts):
    pct = 100 * c / total_count
    info += f"{lo:.1f}–{hi:.1f} keV | {pct:.1f}%\n"

plt.text(
    0.98, 0.95, info,
    ha='right', va='top',
    transform=plt.gca().transAxes,
    fontsize=10,
    bbox=dict(facecolor='white', alpha=0.7, edgecolor='gray')
)

plt.tight_layout()
plt.savefig(f"energy_band_split.png", dpi=140, bbox_inches="tight")
plt.show()

ターミナル出力例

Bands:
0.5–1.7 keV | mean=1.32 keV | count=972 (32.4%)
1.7–2.9 keV | mean=2.17 keV | count=995 (33.1%)
2.9–7.0 keV | mean=4.34 keV | count=1036 (34.5%)

このように minE=0.5, maxE=7.0, nbands=3, dE=0.1 と指定しているため、0.5–7 keV をほぼ均等に 3 分割したエネルギーバンドが得られます。

補足

fluximage で実際に分割する方法

実際にイメージを作成する際には fluximage を利用するのが一般的です。そのとき必要となる representative energy（代表エネルギー） には、本スクリプトがログとして出力する 各バンドの平均エネルギー をそのまま使用できます。

たとえば今回の例では 0.5–1.7 keV | mean=1.32 keV に対応して、fluximage では band=0.5:1.7:1.32 と指定すれば OK です。
（本記事に付属のデータセットには fluximage に必要な他のファイルが含まれていないので、実行は chandra_repro の出力フォルダで行ってください。）

マニュアルバンド（自由に境界を設定したいとき）

均等割りをしたあとに「少しだけ境界を動かしたい」場面があります。
例えば、2.9 keV では中途半端なので 3.0 keV の方が都合が良い、といったケースです。

使用例

custom = [1.7, 3.0]  # 好きな境界

bands, means, counts, edges = split_bands_with_stats(
    energy, minE, maxE,
    mode="manual",
    custom_edges=custom
)

ターミナル出力例

Bands:
0.5–1.7 keV | mean=1.32 keV | count=972 (32.4%)
1.7–3.0 keV | mean=2.22 keV | count=1064 (35.4%)
3.0–7.0 keV | mean=4.44 keV | count=967 (32.2%)

3.0 keV を手動指定しても統計のバランスは大きく崩れず、実務ではこのようなケースで mode="manual" が活用できると思います。

まとめ

本記事では、Chandra イメージング解析で重要となる
「統計の揃ったエネルギーバンド」を自動生成する Python 関数を紹介しました。

累積カウントを用いたシンプルな手法により、
重い最適化を行わずとも 実務上十分にバランスの良いバンド を得られる点が特徴です。

さらに、代表エネルギーの取得や境界の手動調整（manual mode）など、
画像作成やデコンボリューションへ自然に接続できる柔軟な設計になっています。

Chandra 画像解析におけるエネルギーバンド選定が、少しでも軽く・効率的になる助けとなれば幸いです。

あわせて、時間方向で統計を揃える方法については以下の記事も書いていますので、興味があればご覧ください：