はじめに
XRISM衛星のXtendのQL解析プログラムについて紹介します。主目的は pileup の cal なので、明るい点原を想定したツール群で構成されています。pileupの判定ツールだけではなく、そもそもデータにどの程度の影響があるかまでを end-to-end でチェックできることを目的としたソフトウェアになります。
「pileup」効果とは
「pileup」効果とは、複数の光子が検出器でほぼ同時に到達することで、これらが1つのイベントとして記録される現象を指します。このような記録の結果をもとに元の光子数やエネルギー分布を推定することは、逆問題にあたります。この逆問題では、元の光子数分布が複数の異なる組み合わせから生成される可能性があり、解が一意に定まりにくいという性質を持っています。特に、X線CCDの場合は、一画素に信号が留まらず、Grade分岐のパターンが様々存在するため、pileupの複雑性が増すセンスです。望遠鏡の癖や、dust scattering などの効果も絡むと、厳密な推定はさらに難しいです。そのため、実データとシミュレーションの比較検証が大切になります。ソフトウェアの使い方は、基本的に xtend_autorun_docal.py というスクリプトに OBSID を与えるだけで動くようにしてあります。コードは、
にあります。(まだ開発中なので、なんかおかしかったらコードを読んでね、という部分が多いです。)
Resolve の解析については、
を参照ください。
概要説明
xtend_autorun_docal.py の概要について説明します。主には下記のステップと、最終的に結果のHTML出力から構成されます。
-
Step 1: Check pileup
-
xtend_pileup_check_quick.sh
- パイルアップの判定、pileup している場合は region ファイルの生成
-
xtend_pileup_check_quick.sh
-
Step 2: Create region
-
xtend_util_genregion.py
- パイルアップによらず、半径60秒角の円、annulusの region ファイルの生成
-
xtend_util_genregion.py
-
Step 3: Create PHA, RMF, ARF
-
xtend_auto_gen_phaarfrmf.py
- Step1, 2 で生成した領域ファイルに従って、pha,rmf, arf を生成
- pha ファイルは ftgrouppha で optimal binning したものも生成
-
xtend_auto_gen_phaarfrmf.py
-
Step 4: Check PHA, RMF, ARF
-
xrism_util_plot_arf.py
- Step3 で arf ファイルの比較の図を生成
-
xrism_spec_qlfit_many.py
- Step3 生成した pha ファイルを全部 xspec に読み込み、powerlaw 比を生成
-
xrism_util_plot_arf.py
-
HTML生成
-
xrism_autorun_png2html.py
- プログラムの最後に、png ファイル、xspecのlogを全部集めて、html 化する。
-
xrism_autorun_png2html.py
まで自動的に走ります。
- 明るい天体を想定しているので、flickering pixel の除去はしてません。
- 点原を想定して作成してますが、広がったソースの場合でも pileup image は参考になるかもしれません。
- XISの時は、
cellsizeを 3x3 = 9 でしたが、今は 5 にしています。これは今後 cal をして決めていく数字なりますし、pileup 判定で大事な量になります。 - calフェーズなので、全部unzipしたデータを使ってコード生成してます。gzip のままのデータがあるときの動作検証はしてません。
実行結果はHTMLで表示されます。参考例は下記から見えます。(暫定版なのでパスワードで限定中です。)
使い方
PATH を通す
git の rksysoft 以下の resolve, xtend, xrism 以下のファイルをダウンロードして、RESOLVETOOLS=/Users/syamada/work/software/mytool/gitsoft/rksysoft のように、rksysoft にPATHを通して、
# add for resolve tool
export RESOLVETOOLS=/Users/syamada/work/software/mytool/gitsoft/rksysoft
# RESOLVETOOLS/resolve以下のすべてのディレクトリをPATHに追加
for dir in $(find "$RESOLVETOOLS/resolve" -type d); do
PATH="$dir:$PATH"
done
export PATH
# RESOLVETOOLS/xtend以下のすべてのディレクトリをPATHに追加
for dir in $(find "$RESOLVETOOLS/xtend" -type d); do
PATH="$dir:$PATH"
done
export PATH
# RESOLVETOOLS/xrism以下のすべてのディレクトリをPATHに追加
for dir in $(find "$RESOLVETOOLS/xrism" -type d); do
PATH="$dir:$PATH"
done
export PATH
のようにごそっと PATH を通してください。(pileup toolだけに絞った minimal なコードセットは最終的には作った方が良いですが、今はまだ開発中なので。。)
OSBIDを指定して実行する
000125000 というデータがある場合、このOBSIDが見えるディレクトリから、下記を実行する。
000125000 <== この場所で、xtend_autorun_docal.py 000125000 を実行する
----------- auxil/ log/ resolve/ xtend/ <== 下にこれらが存在しているはず。
たとえば、
%ls
auxil log resolve xtend
の４つのディレクトリが ls で見れる場所にいる想定で、
xtend_autorun_docal.py 000125000
と実行する。(xtend_autorun_docal.py -h で使い方がでます。) arfgen にちょっと時間かかりますが、終わると、
A_XRISM_QL_html_000125000/A_QL_html_000125000_checkpileup_v0.html
という html ファイルが生成されます。
たとえば、NGC4151について走らせた場合は、
のような結果になります。
途中生成物は、全て、
000125000/xtend/event_cl/checkpileup_std/
以下に保存されています。
パイルアップの判定だけしたい場合
cl.evt を用いたパイルアップ判定方法
xtend_pileup_check_quick.sh xa300041010xtd_p0300000a0_cl.evt
のように、cleaned event を指定すると、img fits の生成と、pileup 判定ツールが走ります。
img ファイルを用いたパイルアップ判定方法
すでにimageを生成している場合は、
に、image の fits ファイルを与えると、piupup 判定されます。
結果の生成
html の生成
プログラムの最後に、png ファイルを全部集めて、html 化する。
というプログラムで、png を収集して、html の生成をしています。
高速化の tips
OBSID と numphoton を指定して実行する場合
動作確認をしたい場合は、arfgen の numphoton 数を小さくすることで、早く計算が終わるので、-n というオプションで、numphoton を制限することができます。
デフォルトは 100万光子 default=1000000 なので、たとえばこれを1万光子にすると、
xtend_autorun_docal.py 900001010 -n 10000 -s 3 4 -c yes
のように、-n 10000 で光子数を10000個にして、-s 3 4 とすると、arfgen とスペクトルフィットだけ再実行できます。-c yes をつけないと、xaarfgenで指定している xrtevtfile=xtd_raytrace_ptsrc.fits の xtd_raytrace_ptsrc.fitsが存在する限りは、再計算はしないので、-c yes をつけると xtd_raytrace_ptsrc.fits を削除して再計算が走ります。
光子数が少なすぎると、
INFO ::xaxmaarfgen: Insufficient number of raytrace photons
in detector region: check coordinates and region files for errors
before increasing the number of photons- aborting
というメッセージが出て、abort して arf ファイルが生成されないので、注意が必要です。
パイルアップのついての詳細説明
pileup fraction の定義
pileup fraction の式は、
の式(1) に従って、
f_{\mathrm{pl}}(x)=\frac{\sum_{k=2}^{\infty} P(k, x)}{\sum_{k=1}^{\infty} P(k, x)}
を用いています。この式は無限和が分母分子に出てきますが、次のように簡単な式なります。
分子の整理
分子の項 $\sum_{k=2}^{\infty} P(k, x)$ を考えます。ポアソン分布の表式を代入すると：
\sum_{k=2}^{\infty} P(k, x) = \sum_{k=2}^{\infty} \frac{x^k e^{-x}}{k!}.
ポアソン分布の確率の総和の性質より、すべての $k $ についての総和は 1 です：
\sum_{k=0}^{\infty} P(k, x) = 1.
したがって、分子は次のように $k=0$ と $k=1$ の項を除いた形で書けます：
\sum_{k=2}^{\infty} P(k, x) = \sum_{k=0}^{\infty} P(k, x) - P(0, x) - P(1, x).
\sum_{k=2}^{\infty} P(k, x) = 1 - P(0, x) - P(1, x).
ポアソン分布の具体的な表式を用いて $P(0, x)$ と $P(1, x)$ を計算します：
P(0, x) = \frac{x^0 e^{-x}}{0!} = e^{-x},
P(1, x) = \frac{x^1 e^{-x}}{1!} = x e^{-x}.
したがって、分子は次のようになります：
\sum_{k=2}^{\infty} P(k, x) = 1 - e^{-x} - x e^{-x}.
分母の整理
分母の項 $\sum_{k=1}^{\infty} P(k, x)$ を考えます。同様にポアソン分布の表式を代入すると：
\sum_{k=1}^{\infty} P(k, x) = \sum_{k=0}^{\infty} P(k, x) - P(0, x).
ここで $\sum_{k=0}^{\infty} P(k, x) = 1$ を用いると：
\sum_{k=1}^{\infty} P(k, x) = 1 - P(0, x).
さらに、$P(0, x) = e^{-x}$ を代入すると：
\sum_{k=1}^{\infty} P(k, x) = 1 - e^{-x}.
元の式の整理
これまでの結果を用いて、元の式を簡単化します：
f_{\mathrm{pl}}(x) = \frac{\sum_{k=2}^{\infty} P(k, x)}{\sum_{k=1}^{\infty} P(k, x)}.
分子と分母を代入すると：
f_{\mathrm{pl}}(x) = \frac{1 - e^{-x} - x e^{-x}}{1 - e^{-x}}.
丁寧に式変形とポアソン分布の性質を利用することで、簡単な形に変形できます。このように整理することで、分母・分子が有限項の計算で表現され、数値計算にも適した形式となります。下記のコードでは、この式を用いて pileup fraction を計算しています。
def calc_pileupfraction(r):
# Convert input 'r' to a numpy array if it is not already
r = np.asarray(r) # Ensure 'r' is a numpy array for element-wise operations
# Initialize an array of zeros with the same shape as 'r'
result = np.zeros_like(r) # To store the final pile-up fraction values
# Create a mask for elements where 'r' > 0
mask = (r > 0) # Boolean mask to identify valid (positive) values in 'r'
# Calculate the pile-up fraction for elements where 'r' > 0
result[mask] = (1. - (1. + r[mask]) * np.exp(-r[mask])) / (1. - np.exp(-r[mask]))
# Formula explanation:
# (1 - (1 + r) * exp(-r)) / (1 - exp(-r)) calculates the pile-up fraction
# Only applied to elements where 'r' > 0, using the mask
# Return the array of pile-up fractions
return result # 'result' contains the computed values for all elements in 'r'
この関数は、numpy 配列と互換性のある方法で定義された関数 なので、引数の r はスカラーでも、一次元のベクトルでも、２次元配列でも動きます。python の numpy の使い方の大事な部分なので、もし理解が曖昧な方は、
を一読ください。
この性質を使うことで、radial profile でも、image でも同じ関数を使えます。
- radial profile (一次元ベクトル)の pileup fraction の計算
counts_per_cellsize_timedel = radial_profile * cellsize * timedel / (areas * exposure)
radial_pileupfraction = calc_pileupfraction(counts_per_cellsize_timedel)
- image (2次元行列)の pileup fraction の計算
pileupfraction_img = calc_pileupfraction(cellsize_per_exposure)
この pileup fraction の定義において、最後に曖昧さが残ります。平均カウントをどう定義すべきか？という問題です。 CCDの場合は、一つのX線イベントがほぼ1画素だけに信号を落とす装置ではなく、3x3程度の画素数に広がるため、pileup の計算に使うべき平均カウントは何画素分を使うべきか、という点に曖昧さが生まれてしまいます。
現状は、cellsize を 5 にして計算していますが、これは、XISの時の cellsize = 9 を用いると、pileup を overestimate 気味な気がするので、現象論的に 5 にしていますが、cal を行い決めていく量になります。
コードについての補足情報
中心座標の計算
ある領域内の重心を計算しています。
中心座標の初期値は、(x,y)=(1215,1215) としています(オプションで変更できます)。
parser.add_argument('-x', '--x_center', type=int, help='x coordinate of center', default=1215)
parser.add_argument('-y', '--y_center', type=int, help='y coordinate of center', default=1215)
具体的なコードについては、
で紹介しています。
ライトカーブの生成
これで生成した evt ファイルを使って、ライトカーブも簡単に生成したい場合は、次のスクリプトが使える。
xtend_genlc.py ds9_region_0_60.evt -i cygx1
デフォルトで 128秒ビンのライトカーブを生成する。
まとめ
パイルアップは不可逆な現象であり、その推定は本質的に難しいです。サイエンスの目標を見定め、影響を適切に取り除き、光子を最大限活用するためには知恵と経験の積み重ねが大切です。今回は、効率的に多くのデータを解析するため、本ツールを作成しました。現在開発中のため、デバッグ情報などのフィードバックは大歓迎です。