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.

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:])