たとえば右図のように視野中心付近から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)
