aboutsummaryrefslogtreecommitdiffstats
path: root/astro/spectrum
diff options
context:
space:
mode:
authorAaron LI <aly@aaronly.me>2019-10-16 13:21:08 +0800
committerAaron LI <aly@aaronly.me>2019-10-16 13:21:08 +0800
commitf55d3b51de559ea68f5d63bbf1a3caf92e1eb931 (patch)
treeb26465fdedf3f4852c231983ef661bb0fdb632e6 /astro/spectrum
parent1b2c8eb3b653ecc05e332be16a79e740cd679a9f (diff)
downloadatoolbox-f55d3b51de559ea68f5d63bbf1a3caf92e1eb931.tar.bz2
astro: Add sfr_xmm.py by Junhua GU
This script is used to do XMM automatic spectrum extraction and fitting for a list of regions.
Diffstat (limited to 'astro/spectrum')
-rwxr-xr-xastro/spectrum/sfr_xmm.py197
1 files changed, 197 insertions, 0 deletions
diff --git a/astro/spectrum/sfr_xmm.py b/astro/spectrum/sfr_xmm.py
new file mode 100755
index 0000000..a9ec4da
--- /dev/null
+++ b/astro/spectrum/sfr_xmm.py
@@ -0,0 +1,197 @@
+#!/usr/bin/env python
+
+import sys
+import os
+import exceptions
+import re
+
+grppha_comm="aa"
+
+
+def extract_mos(regdes,table_name,temp_pi,bg_table_name,out_bg_name,rmf_name,arf_name,spec_bin_cmd,det='mos1'):
+ expr="'((X,Y) IN %s)'"%regdes
+ out_pi_name = "_a.pi"
+ if not re.search(".pi$",bg_table_name):
+ cmd="evselect"
+ cmd+=" table=%s withspectrumset=yes spectrumset=%s energycolumn=PI spectralbinsize=5 withspecranges=yes specchannelmin=0 specchannelmax=11999 expression=%s "%(bg_table_name,out_bg_name,expr)
+ print(cmd)
+ cmd+=">/dev/null"
+ print("extracting background...")
+ os.system(cmd)
+ else:
+ out_bg_name=bg_table_name
+
+ cmd="evselect"
+ cmd+=" table=%s withspectrumset=yes spectrumset=%s energycolumn=PI spectralbinsize=5 withspecranges=yes specchannelmin=0 specchannelmax=11999 expression=%s"%(table_name,out_pi_name,expr)
+ print(cmd)
+ cmd+=">/dev/null"
+ print("extracting spectrum...")
+ os.system(cmd)
+
+
+ try:
+ os.unlink(temp_pi)
+ except exceptions.OSError:
+ pass
+
+ cmd="grppha infile='%s' outfile='%s' comm='%s & chkey BACKFILE %s & chkey RESPFILE %s & chkey ANCRFILE %s & exit &'"%(out_pi_name,temp_pi,spec_bin_cmd,out_bg_name,rmf_name,arf_name)
+# print cmd
+ cmd+=">/dev/null"
+ os.system(cmd)
+
+
+
+def extract_pn(regdes,table_name,temp_pi,bg_table_name,out_bg_name,rmf_name,arf_name,spec_bin_cmd,det='mos1'):
+ expr="'((X,Y) IN %s)'"%regdes
+ out_pi_name = "_a.pi"
+ if not re.search(".pi$",bg_table_name):
+ cmd="evselect"
+ cmd+=" table=%s withspectrumset=yes spectrumset=%s energycolumn=PI spectralbinsize=5 withspecranges=yes specchannelmin=0 specchannelmax=20479 expression=%s "%(bg_table_name,out_bg_name,expr)
+ print(cmd)
+ cmd+=">/dev/null"
+ print("extracting background...")
+ os.system(cmd)
+ else:
+ out_bg_name=bg_table_name
+
+ cmd="evselect"
+ cmd+=" table=%s withspectrumset=yes spectrumset=%s energycolumn=PI spectralbinsize=5 withspecranges=yes specchannelmin=0 specchannelmax=20479 expression=%s"%(table_name,out_pi_name,expr)
+ print(cmd)
+ cmd+=">/dev/null"
+ print("extracting spectrum...")
+ os.system(cmd)
+
+
+ try:
+ os.unlink(temp_pi)
+ except exceptions.OSError:
+ pass
+
+ cmd="grppha infile='%s' outfile='%s' comm='%s & chkey BACKFILE %s & chkey RESPFILE %s & chkey ANCRFILE %s & exit &'"%(out_pi_name,temp_pi,spec_bin_cmd,out_bg_name,rmf_name,arf_name)
+# print cmd
+ cmd+=">/dev/null"
+ os.system(cmd)
+
+
+
+def make_fit_script(fit_script_template,pi_list,fit_script):
+ fit_script.write("set ff [open %s w]\n"%temp_result_name)
+# print "fit_scriptfit_script"
+ for i in fit_script_template:
+ if re.search("^#PARA",i):
+ para=re.split("\ +",i)
+ fit_script.write("tclout para %s\n"%para[2])
+ fit_script.write('puts $ff "$xspec_tclout %s"\n'%para[1])
+ elif re.search("^#STAT",i):
+ fit_script.write("tclout stat\n")
+ fit_script.write('puts $ff "$xspec_tclout statistic"\n')
+ fit_script.write("tclout dof\n")
+ fit_script.write('puts $ff "$xspec_tclout dof"\n')
+ elif re.search("^#GRPPHA",i):
+ grppha_comm=str.strip(re.sub("#GRPPHA","",i))
+ return grppha_comm
+ elif re.search("^#MODEL",i):
+ fit_script.write(str.strip(re.sub("#MODEL","",i))+"\n")
+ elif re.search("^#DATA",i):
+ for j in range(0,len(pi_list)):
+ fit_script.write("data %d:%d %s\n"%(j+1,j+1,pi_list[j]))
+ elif re.search("^#IGNORE",i):
+ for j in range(0,len(pi_list)):
+ fit_script.write("ignore %d:%s\n"%(j+1,str.strip(re.sub("#IGNORE","",i))))
+ else:
+ fit_script.write(i)
+
+
+ fit_script.write("tclexit\n")
+# print "fdsf"
+
+
+def parse_result(result_file):
+ result=""
+ for i in result_file:
+ str.strip(i)
+ sl=re.split("\ +",i)
+
+ if re.search('statistic',i):
+ result+=sl[0].strip()
+ result+=" "
+ elif re.search('dof',i):
+ result+=sl[0].strip()
+ result+=" "
+ else:
+ result+=sl[0].strip()
+ result+=" "
+ result+=sl[6].strip()
+ result+=" "
+ return result
+
+
+
+
+temp_mos1_pi="_mos1.pi"
+temp_mos2_pi="_mos2.pi"
+temp_pn_pi="_pn.pi"
+temp_mos1_bg="_mos1.bg"
+temp_mos2_bg="_mos2.bg"
+temp_pn_bg="_pn.bg"
+temp_fit_script="_a.tcl"
+temp_result_name="_a.dat"
+
+
+def main(argv):
+ if len(argv)!=15:
+ print("Usage:%s <sample list> <mos1.evt> <mos2.evt> <pn.evt> <mos1.bkg> <mos2.bkg> <pn.bkg> <output file> <mos1.arf> <mos2.arf> <pn.arf> <mos1.rmf> <mos2.rmf> <pn.rmf> <fit script>"%(sys.argv[0]))
+ sys.exit()
+
+
+ sample_list_name=argv[0]
+ mos1_name=argv[1]
+ mos2_name=argv[2]
+ pn_name=argv[3]
+ mos1_bkg=argv[4]
+ mos2_bkg=argv[5]
+ pn_bkg=argv[6]
+ output_name=argv[7]
+ mos1_arf=argv[8]
+ mos2_arf=argv[9]
+ pn_arf=argv[10]
+ mos1_rmf=argv[11]
+ mos2_rmf=argv[12]
+ pn_rmf=argv[13]
+
+ fit_script_template=argv[14]
+
+# print fit_script_template
+
+
+ n=1
+ fit_script=open(temp_fit_script,'w');
+ grppha_comm=make_fit_script(open(fit_script_template),[temp_mos1_pi,temp_mos2_pi,temp_pn_pi],fit_script)
+ #fit_script.write("tclexit")
+
+
+ fit_script.write("tclexit\n")
+ fit_script.close()
+# print grppha_comm
+# sys.exit()
+ outfile=open(output_name,'w')
+ for i in open(sample_list_name):
+ regdesc=str.strip(i)
+ if not re.search('^#NULL',regdesc):
+ extract_mos(regdesc,mos1_name,temp_mos1_pi,mos1_bkg,temp_mos1_bg,mos1_rmf,mos1_arf,grppha_comm)
+ extract_mos(regdesc,mos2_name,temp_mos2_pi,mos2_bkg,temp_mos2_bg,mos2_rmf,mos2_arf,grppha_comm)
+ extract_pn(regdesc,pn_name,temp_pn_pi,pn_bkg,temp_pn_bg,pn_rmf,pn_arf,grppha_comm)
+ print("fitting...")
+ os.system("xspec %s >/dev/null"%temp_fit_script)
+ result=parse_result(open(temp_result_name))
+ outfile.write("%d %s %s"%(n,result,'\n'))
+ print("%d %s"%(n,result))
+ else:
+ outfile.write("\n")
+ print("%d"%(n))
+ n+=1
+ outfile.flush()
+
+
+if __name__=='__main__':
+ main(sys.argv[1:])