スペクトル有効面積補正ツール

広がった天体のスペクトル解析をする際には、source と background での有効面積のちがいが重要になる場合があります。

たとえば右図のように視野中心付近からsource領域を、視野の端からbackgroundをえらぶ場合、off-axis anglesに大きく依存するvignetting effectにより２領域間の有効面積が数十%も変わります。この効果はenergy dependentなため、high energy側ほど「落ち」が著しくなるという厄介な性質をもちます。

右図の場合、background領域でhigh energy側の有効面積をunder estimateしてしまうので、なにも考えず差し引くと、ありもしないhard emissionが浮かび上がってきてしまうことになり、非常に危険です。実際にそういう場面を目撃したことが何度もあります。

これを避けるには、各領域で作成したarf fileをもちいて、有効面積 (specresp) の比を計算し、backgroundの pi fileを補正してやる必要があります。

残念ながら、実際にこの補正をするための公式なftoolは存在しません。

というわけで書きました。

土台となっているのは学生の頃に某先輩が書いてくれたスクリプトですが (Hyodo et al. 2008)、細かいbugを修正し、すざくXIS以外でも (原理的には) つかえるようにし、pythonに移植しました。使用した結果の責任は一切とりませんが、ご自由にお使い下さい。Bugがあれば連絡ください。

使用上の注意：
pythonのmodule, numpyとpyfitsが必要。
結果のまとめをqdp形式で出力するので、そのplotにはqdpが必要 (出力例を下に示します)。
望遠鏡を経ない成分 (=NXB) はあらかじめ引いておくこと。
補正後のbackground spectrumをsource spectrumから引く際、source regionからもNXBを引くことをお忘れなく。

#!/usr/bin/env python

u'''

Written by Makoto Sawada
Version 0.1 (2012/06/12)

mularfratio scales the background spectrum with 
the energy-dependent arf ratio between 
the source and background regions.

ARGUMENTS
bgdarf          ... arf for background
srcarf          ... arf for source
bgdrmf          ... rmf for background
bgdpi           ... original background pi
corrected_bgdpi ... corrected background pi

NOTE
An nxb-subtracted background spectrum is needed.
The pi file must have a RATE column in the 1st extension.

The script requires the pyfits module for FITS I/O 
and the numpy module for calculations.
You need qdp to plot the results.

This is primarily made for Suzaku/XIS spectra, 
although it may work for others.

'''

import pyfits
import numpy as np

eb_acc = 3 #in unit of -log10(E/eV). 1e-{eb_acc} eV accuracy.

def load_arf(arffile):

    hdulist = pyfits.open(arffile)
    arfdata = hdulist[1].data
    en_lo = np.round(arfdata.field('ENERG_LO') * 1e3, eb_acc)
    en_hi = np.round(arfdata.field('ENERG_HI') * 1e3, eb_acc)
    specresp = arfdata.field('SPECRESP')
    hdulist.close()

    return([en_lo,en_hi,specresp])    

def load_rmf(rmffile):

    hdulist = pyfits.open(rmffile)
    extnames = [hdulist[i].name for i in range(len(hdulist))]
    extnum = extnames.index('EBOUNDS')
    ebounds = hdulist[extnum].data
    channel = ebounds.field('CHANNEL')
    e_min = np.round(ebounds.field('E_MIN') * 1e3, eb_acc)
    e_max = np.round(ebounds.field('E_MAX') * 1e3, eb_acc)
    hdulist.close()

    return([channel,e_min,e_max])

def load_pi(pifile):

    pi_hdulist = pyfits.open(pifile)
    pidata = pi_hdulist[1].data
    channel = pidata.field('CHANNEL')
    rate = pidata.field('RATE')

    return([channel,rate,pi_hdulist])

def save_pi(newpifile,newrate,hdulist):

    pidata = hdulist[1].data
    pidata.field('RATE')[:] = newrate
    hdulist.writeto(newpifile,clobber=True)
    hdulist.close()

    return(0)

def result_qdp(src_srs,bgd_srs,bgd_e1s,bgd_e2s,arf_ratios,rmf_e1s,rmf_e2s,corrected_bgdpi):

    emin = min(np.min(bgd_e1s),np.min(rmf_e1s))
    emax = max(np.max(bgd_e2s),np.max(rmf_e2s))

    qdpfile = corrected_bgdpi.replace('.pi','.qdp')
    f = open(qdpfile,'w')

    f.write(
        'skip sing\n'
        'fo r\n cs 1.5\n lw 4\n ti off\n la rot\n'
        'win 1\n'
        'yplot 1 2\n'
        'loc 0.05 0.05 1 0.55\n'
        'la x Energy (eV)\n'
        'la oy SPECRESP (cm\u2\d)\n'
        'r x '+str(emin)+' '+str(emax)+'\n'
        'la f Black: source, Red: background\n'
        'win 2\n'
        'yplot 3 4\n'
        'loc 0.05 0.5 1 1\n'
        'la oy Scaling factor\n'
        'r x '+str(emin)+' '+str(emax)+'\n'
        'la f Green: in arf bin, Blue: in rmf bin\n'
        'lab nx off\n'
        'win all\n'
        )

    for src_sr,bgd_e1,bgd_e2 in zip(src_srs,bgd_e1s,bgd_e2s):
        f.write('%f %f\n' % (bgd_e1,src_sr))
        f.write('%f %f\n' % (bgd_e2,src_sr))
    f.write('no\n')
    for bgd_sr,bgd_e1,bgd_e2 in zip(bgd_srs,bgd_e1s,bgd_e2s):
        f.write('%f %f\n' % (bgd_e1,bgd_sr))
        f.write('%f %f\n' % (bgd_e2,bgd_sr))
    f.write('no\n')
    for src_sr,bgd_sr,bgd_e1,bgd_e2 in zip(src_srs,bgd_srs,bgd_e1s,bgd_e2s):
        f.write('%f %f\n' % (bgd_e1,src_sr/bgd_sr))
        f.write('%f %f\n' % (bgd_e2,src_sr/bgd_sr))
    f.write('no\n')
    for arf_ratio,rmf_e1,rmf_e2 in zip(arf_ratios,rmf_e1s,rmf_e2s):
        f.write('%f %f\n' % (rmf_e1,arf_ratio))
        f.write('%f %f\n' % (rmf_e2,arf_ratio))

    f.close()

    return(0)

def main(bgdarf,srcarf,bgdrmf,bgdpi,corrected_bgdpi):

    ### Loading input files.
    src_e1s,src_e2s,src_srs = load_arf(srcarf)
    bgd_e1s,bgd_e2s,bgd_srs = load_arf(bgdarf)
    rmf_chs,rmf_e1s,rmf_e2s = load_rmf(bgdrmf)
    channels,rates,pi_hdulist = load_pi(bgdpi)

    ### Checking channel consistency of arf files.
    if not np.all(src_e1s == bgd_e1s) or not np.all(src_e2s == bgd_e2s):
        print('Inconsistent channels between the source and background arf files.')
        exit()
        
    ### Checking channel consistency between rmf and pi files.
    if not np.all(channels == rmf_chs):
        print('Inconsistent channels between the rmf and pi files.')
        exit()

    arf_ratios = np.ones(len(channels))

    ### Calculating arf_ratio for each pi channel.
    for channel in channels:
        
        rmf_e1 = rmf_e1s[channel]
        rmf_e2 = rmf_e2s[channel]

        inside_index = np.greater(bgd_e1s, rmf_e1) * np.less(bgd_e1s, rmf_e2)
        arf_inside_es = bgd_e1s[inside_index]
        arf_inside_chs = np.arange(len(inside_index))[inside_index]
        split_num = len(arf_inside_es) + 1
        #print('ch=%f, rmf_ebounds=[%f,%f], split_num=%s' % (channel,rmf_e1,rmf_e2,split_num))

        if split_num == 1:
            ### No arf bin defined.
            arf_ratios[channel] = 1.0

        else:

            if arf_inside_chs[0] == 0:
                arf_ratios[channel] = 1.0
                continue

            arf_prev_ch = arf_inside_chs[0] - 1
            arf_next_ch = arf_inside_chs[-1] + 1
            
            arf_about_chs = np.concatenate([np.array([arf_prev_ch]),arf_inside_chs,np.array([arf_next_ch])])
            arf_about_es = bgd_e1s[arf_about_chs]
            split_ebounds = np.concatenate([np.array([rmf_e1]),arf_inside_es,np.array([rmf_e2])])

            weights = np.ones(split_num)
            src_arfs = np.ones(split_num)
            bgd_arfs = np.ones(split_num)

            for i in range(split_num):

                arf_ch = arf_about_chs[i]
                weights[i] = (split_ebounds[i+1] - split_ebounds[i]) / (rmf_e2 - rmf_e1)
                src_arfs[i] = src_srs[arf_ch]
                bgd_arfs[i] = bgd_srs[arf_ch]

            arf_ratios[channel] = np.sum(weights * src_arfs) / np.sum(weights * bgd_arfs)

        #print('Channel=%d, Ratio=%1.5f' % (channel,arf_ratios[channel]))

    corrected_rates = rates * arf_ratios

    save_pi(corrected_bgdpi,corrected_rates,pi_hdulist)

    result_qdp(src_srs,bgd_srs,bgd_e1s,bgd_e2s,arf_ratios,rmf_e1s,rmf_e2s,corrected_bgdpi)

    return(0)

if __name__ == '__main__':

    import sys

    try:
        bgdarf,srcarf,bgdrmf,bgdpi,corrected_bgdpi = sys.argv[1:]
    except:
        print('bgdarf srcarf bgdrmf bgdpi corrected_bgdpi')
        exit()

    main(bgdarf,srcarf,bgdrmf,bgdpi,corrected_bgdpi)

多段sshを介したMacの画面共有

gateway-Aの下にあるlocal-Aの画面を、gateway-Bの下にあるlocal-B (=目の前の端末)で画面共有したいとき。多段sshをつかいましょう。

参考: SSH多段トンネリング接続でMacを画面共有する

gateway-Bへのsshを介す必要がなければ、上記の参考ページの方法そのままでいけます。 僕の環境では目の前の端末から外に直接出られないので、もう１段sshが必要ですが、基本的にはおなじ方法で大丈夫です。

具体的には以下のようになります。

ssh -t -L 10000:127.0.0.1:5900 \
username@gateway-B "ssh -t -L 5900:127.0.0.1:5900 \
username@gateway-A "ssh  -L 5900:127.0.0.1:5900 username@local-A""

あとは画面共有.appを起動してlocalhost:10000を叩けばokです。 もちろんserver (local-A)での画面共有の許可、sshごとのpassword/passphraseの入力が必要ですのでお気をつけを。

D論お役立ち？スクリプト

○Labelのチェックをおこなう

Thesis全体で使用されている\labelコマンドを列挙し、chapter間の重複をしらべる。
\label{labtype:labname}という形式を仮定。labels.txtにchapter, labtype, labname, 重複の有無を列記。重複がなければ最後のcolumnは'ok'、あれば該当するchapter名がcsvで書かれる。重複がある場合も、初出時は'ok'になる仕様。

#!/usr/bin/env python

skipnames = [] #List tex filenames which should be excluded label search.
outfile = 'labels.txt'
maintex = 'ms.tex' #The main tex filename in which your articles are included.

def add_elem(elem,key,dic):

    if dic.has_key(key):

        dic[key].append(elem)

    else:
        
        dic[key] = [elem]

    return(0)

texnames = []
labels = {}

for line in open(maintex,'r'):

    if line.startswith('\include{') or line.startswith('%\include{'):
        texname = line.replace('\include{','').replace('}\n','').replace('%','')
        if not texname in skipnames:
            texnames.append(texname)

f = open(outfile,'w')

for texname in texnames:

    texfile = texname+'.tex'

    for line in open(texfile,'r'):

        if line.find('\label{') != -1:

            templine = line[line.find('\label{')+len('\label{'):]
            label = templine[:templine.find('}')]
            labtype, labname = label.split(':')
            add_elem(label,texname,labels)
            dupflag = 'ok'
            for texkey in labels.keys():
                if texkey == texname:
                    continue
                lablist = labels[texkey]
                if label in lablist:
                    if dupflag == 'ok':
                        dupflag = texkey
                    else:
                        dupflag += ','+texkey
                else:
                    pass
            f.write('%10s %10s %15s %20s\n' % (texname, labtype, labname, dupflag))
f.close()

○Referenceをsortする

reference.texに記述したreference listをアルファベット順にsortする。
コメントアウトしたものは消されます。Archiveは残すので間違えて消したりしても安心のはず。

#!/usr/bin/env python

import os,datetime

archive = 'ref_archives'
if not os.path.exists(archive):
    os.mkdir(archive)
else:
    pass

reffile = 'reference.tex'
date  = datetime.date.today().isoformat().replace('-','')
biblist = []

i = 0
row = 0
for line in open(reffile,'r'):

    if line.startswith('%') or line.startswith('\\begin') or line.startswith('\\end') or line == '\n':
        pass
    elif line.startswith('\\bibitem'):
        try:
            biblist.append(line.replace('\n',''))
        except:
            print(line)
            exit()
        i += 1
    else:
        try:
            biblist[i-1] = biblist[i-1]+line.replace('\n','')
        except:
            print(line)
            exit()
        pass
    row += 1

os.system('mv '+reffile+' '+archive+'/'+reffile.replace('.tex','_'+date+'.tex'))

f = open(reffile,'w')

bibdoc = {}

for bib in biblist:

    title = bib.replace(']','[').split('[')[1]
    citekey = bib.replace(']','[').split('[')[-1].replace('}','{').split('{')[1]
    content = bib.replace('\\bibitem[','').replace(title,'').replace(']{','').replace(citekey+'} ','')
    author = title.split('(')[0]
    year = title.replace(')','(').split('(')[1]

    bibdoc[title] = bib

keys = bibdoc.keys()

keys.sort()

f.write('\\begin{thebibliography}{}\n')
for key in keys:
    f.write(bibdoc[key]+'\n')
f.write('\end{thebibliography}\n')
f.close()

print('Done.')

○Abbreviationの各章での引用回数をcheckする

description環境で書かれたabbreviation listに各章での略語および元の用語の引用回数をコメントで追記します。
元の用語に関しては複数行にまたがるものは数え落とします。先頭の語のcapitalizeの有無はどっちでも拾います。

#!/usr/bin/env python

import os,datetime

skipnames = [
    'abbreviations',
    'reference',
    'acknowledgement',
    ] #List tex filenames which should be excluded label search.

maintex = 'ms.tex' #The main tex filename in which your articles are included.

texnames = []

for line in open(maintex,'r'):

    if line.startswith('\include{') or line.startswith('%\include{'):
        texname = line.replace('\include{','').replace('}\n','').replace('%','')
        if not texname in skipnames:
            texnames.append(texname)

archive = 'abb_archives'
if not os.path.exists(archive):
    os.mkdir(archive)
else:
    pass

def find_keyinchaps(key):

    chaplab = ''
    ch = 1
    for texname in texnames:
        if texname in skipnames:
            continue
        count = 0
        texfile = texname+'.tex'
        for line in open(texfile,'r'):
            if line.find(key+' ') != -1 or line.find(key+'s ') != -1 or line.find(key+')') != -1:
                count += 1
        if count != 0:
            chaplab += str(ch)+' x'+str(count)+', '
        ch += 1
    return(chaplab[:-2])

def find_continchaps(content):

    content = content.replace('-',' ').lower()

    chaplab = ''
    ch = 1
    for texname in texnames:
        if texname in skipnames:
            continue
        count = 0
        texfile = texname+'.tex'
        for line in open(texfile,'r'):
            line = line.replace('-',' ').lower()
            if line.find(content) != -1:
                count += 1
        if count != 0:
            chaplab += str(ch)+' x'+str(count)+', '
        ch += 1
    return(chaplab[:-2])

abbfile = 'abbreviations.tex'
date  = datetime.date.today().isoformat().replace('-','')
abblist = []

i = 0
row = 0
for line in open(abbfile,'r'):

    if line.startswith('\\chapter') or line.startswith('\\begin') or line.startswith('\\end') or line == '\n':
        pass
    elif line.startswith('\\item'):
        try:
            abblist.append(line.replace('\n',''))
        except:
            print(line)
            exit()
        i += 1
    else:
        try:
            abblist[i-1] = abblist[i-1]+line.replace('\n','')
        except:
            print(line)
            exit()
        pass
    row += 1

os.system('mv '+abbfile+' '+archive+'/'+abbfile.replace('.tex','_'+date+'.tex'))

f = open(abbfile,'w')

abbdoc = {}

for abb in abblist:

    abbkey = abb[abb.find('[')+1:abb.find(']')]
    maxnum = abb.find('%')
    if maxnum == -1:
        maxnum = len(abb)
    content = abb[abb.find('$\\ldots$')+len('$\\ldots$'):maxnum].strip(' ')

    abbdoc[abbkey] = content

keys = abbdoc.keys()

keys.sort()

f.write('\\chapter{Abbreviation}\n')
f.write('\\begin{description}\n')
for key in keys:
    usedchaps = find_keyinchaps(key)
    usedchaps2 = find_continchaps(abbdoc[key])

    f.write('\\item [%s] $\\ldots$ %s %s\n' % (key, abbdoc[key], '% '+usedchaps+' ; '+usedchaps2))

f.write('\\end{description}\n')
f.close()

print('Done.')