Arf_BG_corrector.py

お察し下さい (pyfits練習) 第2弾。前よりは多少ファイルI/Oの書き方が賢くなっている...かも。あくまで個人的に作ったスクリプト。

SCRIPT_SHELL="python" 
SCRIPT_NAME="Arf_BG_corrector.py"
SCRIPT_HELP_ARG="1) Input Rate Bg PI file made by mathpha 2) Src Arf 3) Bg Arf 4) Output Bg PI file"
SCRIPT_HELP_SENTENCE="Modiefy Bg PI Channel with Src / Bg  arf ratio."
SCRIPT_NUM_ARG=4

###IMPORT###
import sys
import pyfits
import os
###HELP###

if len(sys.argv[1:])!=SCRIPT_NUM_ARG:
    print 'Name: '+SCRIPT_NAME
    print 'Arguments: '+SCRIPT_HELP_ARG
    print 'Explanation: '+SCRIPT_HELP_SENTENCE
    sys.exit()


###MAIN###
ev_per_ch=3.65

bg_pif=str(sys.argv[1])
src_arf=str(sys.argv[2])
bg_arf=str(sys.argv[3])
out_pif=str(sys.argv[4])

def main() :
    
    
    if os.path.isfile(out_pif) :
        os.remove(out_pif)
    
    data_src_arf = (pyfits.open(src_arf))['SPECRESP'].data
    data_bg_arf = (pyfits.open(bg_arf))['SPECRESP'].data
   
    ratio=[] 
    for n in range(len(data_src_arf)) :
        ratio.append([data_src_arf['ENERG_LO'][n],data_src_arf['ENERG_HI'][n],data_src_arf['SPECRESP'][n] / data_bg_arf['SPECRESP'][n]])
    
    hdulist=pyfits.open(bg_pif)
    
    sp=hdulist['SPECTRUM'].data
    ch_num=len(hdulist['SPECTRUM'].data)
    
    hdulist['SPECTRUM'].header['COMMENT']='Modefied with Arf_BG_corrector.py'
    hdulist['SPECTRUM'].header['COMMENT']='Src arf : '+src_arf
    hdulist['SPECTRUM'].header['COMMENT']='Bg arf : '+bg_arf
    hdulist['SPECTRUM'].header['COMMENT']='Bg pi : '+bg_pif

        
    for n in xrange(ch_num) :
        enegy=ev_per_ch*n/1e3
        for erange in ratio :
            if erange[0] <= enegy and enegy < erange[1] :
                sp['RATE'][n]=sp['RATE'][n]*erange[2]
                sp['STAT_ERR'][n]=sp['STAT_ERR'][n]*erange[2]
                break

    hdulist.writeto(out_pif)

main()

PI Shifter for Suzaku XIS

お察し下さい。個人的に作ったスクリプト。pyfitsの練習というか。

SCRIPT_SHELL="python" 
SCRIPT_NAME="PiShift.py"
SCRIPT_HELP_ARG="1) input Count PI file made by mathpha 2)output PI file 3)c0 4)c1 5)c2"
SCRIPT_HELP_SENTENCE="Modify PI Channel with f(x)=c2*x**2+c1*x+c0. Only for Suzaku XIS."
SCRIPT_NUM_ARG=5

###IMPORT###
import sys
import pyfits
import os
###HELP###

if len(sys.argv[1:])!=SCRIPT_NUM_ARG:
    print 'Name: '+SCRIPT_NAME
    print 'Arguments: '+SCRIPT_HELP_ARG
    print 'Explanation: '+SCRIPT_HELP_SENTENCE
    sys.exit()


###MAIN###
inf=str(sys.argv[1])
outf=str(sys.argv[2])
c0=float(sys.argv[3])
c1=float(sys.argv[4])
c2=float(sys.argv[5])


ch_num=4096

hdulist = pyfits.open(inf)

if os.path.isfile(outf) :
    os.remove(outf)

sp=hdulist[1].data

hdulist[1].header['COMMENT']='Modefied with PiSift.py'
hdulist[1].header['COMMENT']='c0='+str(c0)
hdulist[1].header['COMMENT']='c1='+str(c1)
hdulist[1].header['COMMENT']='c2='+str(c2)

def f(x) :
    x=float(x)
    return c2*x**2+c1*x+c0


def main() :
    new_counts=[]
    for n in range(ch_num) :
        new_counts.append(0.0)  
        
    old_counts=sp['COUNTS']
    
    for n in range(ch_num) :
        new_ch0=f(n)
        new_ch1=f(n+1)
        if new_ch1 < new_ch0 :
            sys.exit('Correction function is not monotone.')
        elif new_ch0 > 0 and new_ch1 < ch_num :
            for m in range(int(new_ch0+1),int(new_ch1)) :
                new_counts[m]+=old_counts[n]/(new_ch1-new_ch0)

            new_counts[int(new_ch0)]+=old_counts[n]/(new_ch1-new_ch0)*(int(new_ch0+1)-new_ch0)
            new_counts[int(new_ch1)]+=old_counts[n]/(new_ch1-new_ch0)*(new_ch1-int(new_ch1))
         
  

    for n in range(ch_num) :
        sp['COUNTS'][n]=int(round(new_counts[n]))
        sp['STAT_ERR'][n]=float(new_counts[n]**0.5)

    hdulist.writeto(outf)

main()

各XIS (0,1,3) とHXDのNormが食い違う理由

XISとHXDのNormが食い違う理由についての考察?のようなもの。
後輩から質問を受けたのでまとめた。
間違えているかもしれないので注意。指摘していただけると助かります。


-------

以下、簡単に各XIS (0,1,3) とHXDのNormが食い違う理由になりうる較正上の
不定性について簡単に書いておきます。

1. XISのnon-X-ray backgroundの不定性
宇宙線由来のバックグラウンドは観測（衛星軌道）毎・検出器毎、更に検出器上の
場所毎に異なります。XIS 0,1,3間のnormalizationが観測毎に異なるのは主にこのNXBの
揺らぎによる筈です。暗い天体のハードな帯域 (5-10 keVあたり) のfluxはこの不
定性の影響を大きく受けます。 BGを視野内の暗い所から引くという解析では、
この不定性の影響が大きくなります。 xisnxbgenというftoolを使い、「(1) SRCとBG
の領域毎にNXBのスペクトルを作り、(2) SRCとBGの各々NXBを差し引いた上で更に
SRCからBGを引く」といったことを行えば、NXBの不定性は小さくなりますが少し
統計誤差が大きくなります。

2. XISのcontaminationの較正の不定性
XISには可視光遮断フィルターがついていますが、この表面にX線を吸収してしまう汚染
物質（ゴム等からアウトガス）が付着しつつあることが知られています。この汚染物質は
観測時期が進むにつれ増加しており、検出器毎、更に検出器上の場所毎に異なります。
この汚染物質量の較正の不定性はXIS 0,1,3間のnormalizationが異なる原因になります。
ソフトな帯域 (0.5-2 keVあたり) のfluxはこの不定性の影響を大きく受けます。


3. XRTの望遠鏡の較正の不定性
XRTの有効面積レスポンス (arf) を計算するモンテカルロシミュレーション
 (xissimarfgen) 内では設計図どおりの理想的な望遠鏡の形状が仮定されてい
ます。しかし、実際には打ち上げ (+経年劣化?) によりXRTの形状がいくらか変化し
ていることが分かっており、XIS 0,1,3間、およびPINとのnormalizationの
食い違いの原因になります。この不定性はXIS 0,1,3で異なり、更に光軸（視野）
中心から離れるにつれて大きくなります。

4. 衛星の姿勢の揺れの不定性
観測中に「すざく」衛星の姿勢が揺れ、それが補正できないために天体が
光軸中心からずれてしまい、有効面積の見積もりを間違えてしまうことが
わかっています。XIS 0,1,3間、およびPINとのnormalizationの食い違いの
原因になります。この不定性は光軸（視野）中心から離れるにつれて大き
くなります。これまで割と良く補正されていたのですが、最近になってまた
姿勢の揺れが問題になりつつあるようです。aeattcor2というftoolを使うと少し
改善されるかもしれません。

XISRMFPARAM

サワDに対するメールより

ae_xi?_rmfparam_*.fits
の中を見るとParam_s*_p* という32行のパラメータ列
が並んでいると思います。これらがレスポンスの
詳細情報（メインピークの幅とメインピークに対
する他コンポーネントの幅・強度）を示しています。
各行の意味については
src/function/xisRespUtil/xisRespUtil.c
を読んで下さい。（これ以外ドキュメントがなさ
そうです・・・）

現在はメインピークの幅のみを変更しています。
Param_s*_p* の
2行目=s12
3行目=s13
4行目=s14
と定義してE(keV)に対して

FWHM(eV)=(S12+S13*E+S14*E**2)**0.5

です。(Param_s*_p* の1行目は1/(3.65*2.35)でPIとSigmaに変換するための係数
です）Eの単位はkeV、計算結果はeVなので注意して下さい。

なぜ、Param_s*_p*が複数あるかというと
s*は各セグメント毎にp*は同一セグメント内の場所(4カ所)毎に
異なったレスポンスを設定できるからです。
現在(ae_xi?_rmfparam_20080901.fits)は同一センサー
内では16列に全て同じ値を詰めています。

MITによるSpaced-row Charge Injection (SCI) の本論文

みんな Uchiyama et al. (2007) or (2009) を引いてくれるのを嬉しいのだが、あれはsaw-tooth correctionの論文であって、「すざく」の機上SCI論文としては下記を引くべき。

• CCD Charge Injection Structure at Very Small Signal Levels

Prigozhin, Gregory; Burke, Barry; Bautz, Marshall; Kissel, Steve; Lamarr, Beverly
IEEE Transactions on Electron Devices, vol. 55, issue 8, pp. 2111-2120 (2008)

\bibitem[Prigozhin et al.(2008)]{2008ITED...55.2111P} Prigozhin, G., Burke, B., Bautz, M., Kissel, S.,  Lamarr, B.\ 2008, IEEE Transactions on Electron Devices, 55, 2111

MITの方が主著者で、査読付きのIEEE論文になっている、「筋としては」SCIの機上運用の本論文として引くべき一番の論文。…CIそのものについてが主でSCIの記述は少ないけど。

• Mitigating CCD radiation damage with charge injection: first flight results from Suzaku

Bautz, M. W.; LaMarr, B. J.; Miller, E. D.; Kissel, S. E.; Prigozhin, G. Y.; Burke, B. E.; Gregory, J. A.; Uchiyama, H.; Hyodo, Y.; Yamaguchi, H.; Mori, H.; Tsuru, T.; Matsumoto, H.; Koyama, K.; Torii, K.; Katsuda, S.; Hasuike, K.; Nakajima, H.; Hayashida, K.; Tsunemi, H.; Murakami, H.; Ozaki, M.; Dotani, T.

UV, X-Ray, and Gamma-Ray Space Instrumentation for Astronomy XV. Edited by Siegmund, Oswald H. Proceedings of the SPIE, Volume 6686, pp. 66860Q-66860Q-11 (2007)

\bibitem[Bautz et al.(2007)]{2007SPIE.6686E..23B} Bautz, M.~W., LaMarr,
B.~J., Miller, E.~D., et al.\ 2007, \procspie, 6686

「すざく」でのSCI機上運用についての情報、見せたい（見せるべき）図などがまとまったBautz先生のSPIE論文。SCIの機上運用の論文としては内容的に最適。学生さんにも読んでもらいたいところ。…SPIE論文な上、Astro-phにも上がっていないので入手しにくいのが玉に瑕…。


でもUchiyama et al.も引いてくれると嬉しい。
\bibitem[Uchiyama et al.(2009)]{2009PASJ...61S...9U} Uchiyama, H., et al. 2009, \pasj, 61, 9