xspec のビンまとめを見た目通りにする方法

はじめに

xspec を用いて、見た目と同じように xspec でビンまとめする方法です。スペクトルの割り算や比較をしたいときに、ビン幅ばピタッと同じであると何かと便利なので、そういう用途が一番多いかと思います。

ビンまとめは、データをいろんな目で眺めるために大切な操作で、自分の目でデータの質感を感じ取り、最終的に論文などに乗せるビンニングを決めるので、いろんな角度からデータを見てみましょう。その上で、その見た目通りにデータをビンまとめを行うには、xspec の場合は、grppha というツールを使う必要がありますが、それが、xspec の setplot rebin と完璧に一致しない、という問題があるので、それを回避する方法を紹介します。

統計の説明の中で、「ビンまとめは不要」、や」ビンまとめしない方が良い」、を見た人で初学者の人は混乱するかもしれません。そういう主張は最小化の評価方法やモデルの計算方法、統計を語る文脈に依存するので注意してください。基本は適切にビンまとめしてフィットを行い、カイ２乗値を出す、なので、まずはそれをやって見ましょう。

下記に用いたファイル一式を置いています。

別の方法

xspec については、

などを参考にしてください。

最近では、ftgrouppha というツールもあります。

xspec ユーザーは、ftgrouppha を用いて、SPEX の rbin と同じアルゴリズムを、

grouptype=opt sets up a binning using the Kaastra & Bleeker (2016; A&A 587, A151) optimal binning scheme. Note that this option ignores any minchannel and maxchannel set.

grouptype=opt というオプションで実現できます。

その他、ftrbnrmf というツールで、レスポンスのビンまとめができます。

SPEX ユーザーは、

で、optimal binning できます。

ステップ１ : xspec で ch が詰まった qdp ファイルの作成

ビンまとめを全くしていない pha ファイルを使います。

(py2021) [syamada] $ xspec             [~/work/ana/rikkyo/2023/horio/qiita_autogrppha]

		XSPEC version: 12.11.0
	Build Date/Time: Tue Aug  4 22:49:16 2020

XSPEC12>data 1:1 ni0100320106_0mpu7_time_1.pha 

1 spectrum  in use
 
Spectral Data File: ni0100320106_0mpu7_time_1.pha  Spectrum 1
Net count rate (cts/s) for Spectrum:1  1.887e+03 +/- 6.486e-01
 Assigned to Data Group 1 and Plot Group 1
  Noticed Channels:  1-1501
  Telescope: NICER Instrument: XTI  Channel Type: PI
  Exposure Time: 4487 sec
 Using fit statistic: chi
 No response loaded.

***Warning!  One or more spectra are missing responses,
               and are not suitable for fit.
XSPEC12>resp 1 ni0100320106mpu7.rmf 
arResponse successfully loaded.
XSPEC12>arf 1 ni0100320106mpu7.arf 
Arf successfully loaded.
XSPEC12>cpd /xs 
XSPEC12>setp e 
XSPEC12>ignore **-4.0 10.0-** 
   400 channels (1-400) ignored in spectrum #     1
   502 channels (1000-1501) ignored in spectrum #     1

XSPEC12>model power 

Input parameter value, delta, min, bot, top, and max values for ...
              1       0.01(      0.01)         -3         -2          9         10
1:powerlaw:PhoIndex>
              1       0.01(      0.01)          0          0      1e+20      1e+24
2:powerlaw:norm>

========================================================================
Model powerlaw<1> Source No.: 1   Active/On
Model Model Component  Parameter  Unit     Value
 par  comp
   1    1   powerlaw   PhoIndex            1.00000      +/-  0.0          
   2    1   powerlaw   norm                1.00000      +/-  0.0          
________________________________________________________________________


Fit statistic  : Chi-Squared                534488.8     using 599 bins.

Test statistic : Chi-Squared                534488.8     using 599 bins.
 Null hypothesis probability of 0.0e+00 with 597 degrees of freedom
 Current data and model not fit yet.
XSPEC12>
XSPEC12>
XSPEC12>fit
                                   Parameters
Chi-Squared  |beta|/N    Lvl    1:PhoIndex        2:norm
12788.8      13860        -2       1.28113      0.832525
4527.51      69080        -2       1.42598       1.11569
2752.21      33302.6      -3       1.57944       1.45523
1573.85      27584.4      -4       1.58864       1.53277
1573.45      527.257      -5       1.58826       1.53283
1573.45      0.347935     -2       1.58827       1.53284
==============================
 Variances and Principal Axes
                 1        2  
 3.5526E-07| -0.9323   0.3617  
 1.8734E-04|  0.3617   0.9323  
------------------------------

========================
  Covariance Matrix
        1           2   
   2.481e-05   6.305e-05
   6.305e-05   1.629e-04
------------------------

========================================================================
Model powerlaw<1> Source No.: 1   Active/On
Model Model Component  Parameter  Unit     Value
 par  comp
   1    1   powerlaw   PhoIndex            1.58827      +/-  4.98144E-03  
   2    1   powerlaw   norm                1.53284      +/-  1.27623E-02  
________________________________________________________________________


Fit statistic  : Chi-Squared                 1573.45     using 599 bins.

Test statistic : Chi-Squared                 1573.45     using 599 bins.
 Null hypothesis probability of 5.58e-89 with 597 degrees of freedom
XSPEC12>pl
XSPEC12>pl ld delc
XSPEC12>pl ld ra
XSPEC12>setp rebin 45 40 

で rebin しましょう。ここでは、setp rebin 45 sigma 40 counts でビンまとめになります。

を参考にしてください。

setp com col 2 on 2

モデルを赤線にして、モデルの形状を確認する。

setp com r y2 0.8 1.2

y2軸の範囲を 0.8 から 1.2 にして、ビンまとめの良し悪しを判断する。
必要なら、setp rebin のパラメータを変更してやり直す。

XSPEC12>iplot
PLT> hard spec_rebin_45_40.ps/cps
で、図を保存しておく。

このような図ができればOK。

XSPEC12>setp list
   1: col 2 on 2 
   2: r y2 0.8 1.2 
   3: hard spec_rebin_45_40.ps/cps 
XSPEC12>setp delete 1-3 

で、リストを空にします。

XSPEC12>model none
XSPEC12>notice **-**
  1501 channels (1,1501) noticed in spectrum #     1

モデルを空にして、ignore したエネルギーも元に戻します。

これで、横軸が ch で、カットが掛かってないスペクトルになるのを確認します。

このようになればOKです。

XSPEC12>iplot
PLT> we ch_rebin_45_40

として、ch_rebin_45_40.qdp を作成します。これで、ビンまとめされた ch のヒストグラムのテキストファイルができました。

ステップ2 : grppha の sh script を生成

python gen_grpphacode_fromqdp_v1.py ch_rebin_45_40.qdp

を実行すると、 sh script (grppha_tmp.sh) が生成される。

gen_grpphacode_fromqdp_v1.py の中身はこちらです。

gen_grpphacode_fromqdp_v1.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import argparse

# Constants
HEADER_LINES = 3

def create_grppha_script(qdp_filename, output_script_name):
    """
    Create a grppha script based on the input QDP file.
    """
    cr_raw = []

    with open(qdp_filename, 'r') as fin:
        for i, oneline in enumerate(fin):
            if i < HEADER_LINES:
                continue
            data = list(map(float, oneline.split()))
            cr_raw.append(data)

    with open(output_script_name, "w") as fout:
        fout.write("#!/bin/sh\n")
        fout.write("infile=$1\n")
        fout.write("outfile=$2\n")
        fout.write("grppha << EOF\n")
        fout.write("$infile\n")
        fout.write("$outfile\n")

        for i in range(len(cr_raw) - 1):
            bounds = (cr_raw[i][0] - cr_raw[i][1], cr_raw[i][0] + cr_raw[i][1] - 1, 2 * cr_raw[i][1])
            fout.write("group %12.0d %12.0d %12.0d\n" % bounds)

        last_bound = (cr_raw[-1][0] + cr_raw[-1][1], cr_raw[-1][0] + 5 * cr_raw[-1][1] - 1, 2 * cr_raw[-1][1])
        fout.write("group %12.0d %12.0d %12.0d\n" % last_bound)
        fout.write("exit\n")
        fout.write("quit\n")
        fout.write("EOF\n")

    print(f"{output_script_name} is created")


def main():
    parser = argparse.ArgumentParser(description='Generate grppha script from QDP file.')
    parser.add_argument('qdp', type=str, help='Re-binned QDP file name (need Channel for vertical axis)')
    args = parser.parse_args()
    print(f"QDP filename = {args.qdp}")

    create_grppha_script(args.qdp, "grppha_tmp.sh")


if __name__ == "__main__":
    main()

ステップ3 : grppha でビンまとめする

ビンまとめを全くしていない pha ファイル (ni0100320106_0mpu7_time_1.pha)をここでも使います。

./grppha_tmp.sh ni0100320106_0mpu7_time_1.pha

実行すると、

 ......  Upper channel reset to last channel : 1500
 ......  Channel 1500 is less than channel 1501 !
GRPPHA[group         1501         1518            9] exit
GRPPHA[exit] quit
 ** grppha 3.1.0 completed successfully

のように、completed successfully を確認する。

ステップ4 : ステップ１の見た目と比較して確認する

最後に、xspec で生成された pha ファイルと、ステップ１で setplot rebin を比較しましょう。