diff options
| -rwxr-xr-x | astro/spectrum/sfr_xmm.py | 197 |
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:]) |