diff options
| -rw-r--r-- | xspec/xspec_bkgcorr_v2.tcl | 370 | ||||
| -rw-r--r-- | xspec/xspec_coolfunc_v2.tcl | 246 |
2 files changed, 616 insertions, 0 deletions
diff --git a/xspec/xspec_bkgcorr_v2.tcl b/xspec/xspec_bkgcorr_v2.tcl new file mode 100644 index 0000000..61a09f1 --- /dev/null +++ b/xspec/xspec_bkgcorr_v2.tcl @@ -0,0 +1,370 @@ +##################################################################### +## XSPEC Tcl script +## +## Task: +## analysis chandra acis background components +## output model fitting results, then generate +## fake spectra to correct the blanksky background +## +## LIweitiaNux <liweitianux@gmail.com> +## August 1, 2012 +## +## NOTES: needs XSPEC v12.x +## +## ChangeLogs: +## v2, 2012/08/14, LIweitiaNux +## skip calc `error' of params +## change variable names +## improve the final comparison graph +## replace `trim' with `regexp' +##################################################################### + +## about {{{ +set NAME "xspec_bkgcorr_v2.tcl" +set VERSION "v2, 2012-08-14" +## about }}} + +## set basic variables {{{ +## tclout: create tcl var from current state +set tcl_date [ exec date ] +tclout filename 1 +set spec_name [ exec basename $xspec_tclout ] +# `trim' is not reliable: `lbkg_grp.pi' --> `lbkg_gr' +# set spec_rootname "[ string trimright $spec_name {.pi} ]" +regexp -- {(.+)\.(pi)} $spec_name match spec_rootname fext +# get `rootname' from `spec_name' +regexp -- {(.+)_(bin|grp)([0-9]*)\.(pi|fits)} $spec_name match rootname grptype grpval fext +puts "## rootname: $rootname" +puts "## grptype: $grptype" +puts "## grpval: $grpval" +puts "## file extension: $fext" +# backgrnd name +tclout backgrnd 1 +set back_name $xspec_tclout +# set back_rootname "[ string trimright $back_name {.pi} ]" +regexp -- {(.+)\.(pi)} $back_name match back_rootname fext +# regexp -- {(.+)\.(pi)} $spec_name match spec_rootname fext +# `specscal': scale factor to create fake spectrum, +# use a large factor (i.e. 10) to get a good statistic for final spectrum +set specscal 10.0 +## basic variables }}} + +## save current xspec fitting results {{{ +set xspec_outroot "xspec_${rootname}_fit" +# `save all', to save a xcm script +if {[ file exists ${xspec_outroot}.xcm ]} { + exec mv -fv ${xspec_outroot}.xcm ${xspec_outroot}.xcm_bak +} +save all "${xspec_outroot}.xcm" +# writefits: (tcl scripts, v12.x) +#if {[ file exists ${xspec_outroot}.fits ]} { +# exec mv -fv ${xspec_outroot}.fits ${xspec_outroot}.fits_bak +#} +#writefits "${xspec_outroot}.fits" +## save & writefits }}} + +## set output file {{{ +set out_fn "bkgcorr_${rootname}.log" +if {[ file exists ${out_fn} ]} { + exec mv -fv ${out_fn} ${out_fn}_bak +} +set out_fd [ open ${out_fn} w ] +## output file }}} + +## header, output basic process info {{{ +puts $out_fd "# process by script: ${NAME}" +puts $out_fd "# version: ${VERSION}" +puts $out_fd "#" +puts $out_fd "# process date: $tcl_date" +puts $out_fd "#" +puts $out_fd "" +## header }}} + +## save basic info about current fitting {{{ +# files in use +puts $out_fd "data $spec_name" +tclout response 1 +set rmf $xspec_tclout +puts $out_fd "response $rmf" +tclout arf 1 +set arf $xspec_tclout +puts $out_fd "arf $arf" +puts $out_fd "backgrnd $back_name" +puts $out_fd "" +# exposure time, backscale +tclout expos 1 s +scan $xspec_tclout "%f" expos_src +tclout expos 1 b +scan $xspec_tclout "%f" expos_bkg +puts $out_fd "# src/bkg exptime: $expos_src/$expos_bkg" +tclout backscal 1 s +scan $xspec_tclout "%f" backscal_src +tclout backscal 1 b +scan $xspec_tclout "%f" backscal_bkg +puts $out_fd "# src/bkg backscal: $backscal_src/$backscal_bkg" +puts $out_fd "" +# model, and basic fitting results +tclout model +puts $out_fd "model $xspec_tclout" +tclout weight +puts $out_fd "# weight: $xspec_tclout" +tclout stat +scan $xspec_tclout "%f" chisq +tclout dof +scan $xspec_tclout "%d" dof +set rechisq_cmd { format "%.5f" [ expr { $chisq / $dof } ] } +set rechisq [ eval $rechisq_cmd ] +puts $out_fd "# chisq/dof = $chisq/$dof = $rechisq" +tclout noticed 1 +puts $out_fd "# noticed channel: $xspec_tclout" +tclout noticed energy 1 +puts $out_fd "# noticed energy (keV): $xspec_tclout" +puts $out_fd "" +## basic fitting info }}} + +## model fitting results {{{ +puts $out_fd "# fitting results of each components" +puts $out_fd "# errors are emitted for simplicity" +puts $out_fd "#" +puts $out_fd "# name value sigma_err" +puts $out_fd "" +# apec_1, Galactic X-ray background, soft-soft +tclout param 1 +scan $xspec_tclout "%f" temp_a1 +puts $out_fd "temp_a1: $temp_a1" +tclout param 2 +scan $xspec_tclout "%f" abund_a1 +puts $out_fd "abund_a1: $abund_a1" +tclout param 3 +scan $xspec_tclout "%f" redshift_a1 +puts $out_fd "redshift_a1: $redshift_a1" +tclout param 4 +scan $xspec_tclout "%f" norm_a1 +tclout sigma 4 +set norm_sigma_a1 $xspec_tclout +puts $out_fd "norm_a1: $norm_a1 $norm_sigma_a1" +puts $out_fd "" +# apec_2, Galactic X-ray Background, soft +tclout param 5 +scan $xspec_tclout "%f" temp_a2 +puts $out_fd "temp_a2: $temp_a2" +tclout param 6 +scan $xspec_tclout "%f" abund_a2 +puts $out_fd "abund_a2: $abund_a2" +tclout param 7 +scan $xspec_tclout "%f" redshift_a2 +puts $out_fd "redshift_a2: $redshift_a2" +tclout param 8 +scan $xspec_tclout "%f" norm_a2 +tclout sigma 8 +set norm_sigma_a2 $xspec_tclout +puts $out_fd "norm_a2: $norm_a2 $norm_sigma_a2" +puts $out_fd "" +# wabs +tclout param 9 +scan $xspec_tclout "%f" wabs +puts $out_fd "wabs: $wabs" +puts $out_fd "" +# powerlaw, Cosmological X-ray Background, hard +tclout param 10 +scan $xspec_tclout "%f" gamma_po +puts $out_fd "gamma_po: $gamma_po" +tclout param 11 +scan $xspec_tclout "%f" norm_po +tclout sigma 11 +set norm_sigma_po $xspec_tclout +puts $out_fd "norm_po: $norm_po $norm_sigma_po" +puts $out_fd "" +# apec_3, gas emission from object source +tclout param 12 +scan $xspec_tclout "%f" temp_gas +tclout sigma 12 +set temp_sigma_gas $xspec_tclout +puts $out_fd "temp_gas: $temp_gas $temp_sigma_gas" +tclout param 13 +scan $xspec_tclout "%f" abund_gas +tclout sigma 13 +set abund_sigma_gas $xspec_tclout +puts $out_fd "abund_gas: $abund_gas $abund_sigma_gas" +tclout param 14 +scan $xspec_tclout "%f" redshift_gas +puts $out_fd "redshift_gas: $redshift_gas" +tclout param 15 +scan $xspec_tclout "%f" norm_gas +tclout sigma 15 +set norm_sigma_gas $xspec_tclout +puts $out_fd "norm_gas: $norm_gas $norm_sigma_gas" +puts $out_fd "" +## fittins results }}} +# end output results +close $out_fd + +## prepare data for fake spectrum {{{ +# see `specscal' in top `basic variables' section +# apec_1 +if {$norm_a1 < 0} { + set pm_a1 "-" +} else { + set pm_a1 "+" +} +set pars_a1 "$temp_a1 & $abund_a1 & $redshift_a1 & [ expr abs($norm_a1)*$specscal ]" +set model_a1 "apec" +# apec_2 +if {$norm_a2 < 0} { + set pm_a2 "-" +} else { + set pm_a2 "+" +} +set pars_a2 "$temp_a2 & $abund_a2 & $redshift_a2 & [ expr abs($norm_a2)*$specscal ]" +set model_a2 "apec" +# powerlaw +if {$norm_po < 0} { + set pm_po "-" +} else { + set pm_po "+" +} +set pars_po "$wabs & $gamma_po & [ expr abs($norm_po)*$specscal ]" +set model_po "wabs*powerlaw" +# gas emssion +# norm of `gas component' cannot be negative +set pars_gas "$wabs & $temp_gas & $abund_gas & $redshift_gas & [ expr abs($norm_gas)*$specscal ]" +set model_gas "wabs*apec" +## prepare data }}} + +## functions to load model, fake spectrum {{{ +proc tcl_model {model_str pars args} { + model none + model $model_str & $pars & /* +} +proc tcl_fakeit {rmf arf fakedata exptime args} { + data none + fakeit none & $rmf & $arf & y & & $fakedata & $exptime & /* +} +## functions }}} + +## set fake spectrum and check previous files {{{ +set fake_a1 "fake_${rootname}_a1.pi" +set fake_a2 "fake_${rootname}_a2.pi" +set fake_po "fake_${rootname}_po.pi" +set fake_gas "fake_${rootname}_gas.pi" +if {[ file exists $fake_a1 ]} { + exec mv -fv $fake_a1 ${fake_a1}_bak +} +if {[ file exists $fake_a2 ]} { + exec mv -fv $fake_a2 ${fake_a2}_bak +} +if {[ file exists $fake_po ]} { + exec mv -fv $fake_po ${fake_po}_bak +} +if {[ file exists $fake_gas ]} { + exec mv -fv $fake_gas ${fake_gas}_bak +} +## fake spectrum }}} + +## generate fake spectrum {{{ +# blanksky, apec_1 +tcl_model $model_a1 $pars_a1 +tcl_fakeit $rmf $arf $fake_a1 $expos_bkg +# blanksky, apec_2 +tcl_model $model_a2 $pars_a2 +tcl_fakeit $rmf $arf $fake_a2 $expos_bkg +# blanksky, wabs*powerlaw +tcl_model $model_po $pars_po +tcl_fakeit $rmf $arf $fake_po $expos_bkg +# src gas, wabs*apec +tcl_model $model_gas $pars_gas +tcl_fakeit $rmf $arf $fake_gas $expos_src +## fake spectrum }}} + +## background correction {{{ +## blanksky {{{ +set bbkg_expr "${specscal}*${back_name} ${pm_a1}${fake_a1} ${pm_a2}${fake_a2} ${pm_po}${fake_po}" +set bbkg_outf "bkgcorr_blanksky_${rootname}.pi" +set bbkg_back [ expr { $backscal_bkg * $specscal } ] +set bbkg_cmm1 "corrected background spectrum, blanksky based" +set bbkg_cmm2 "norm*${specscal}, backscal*${specscal}, properr=no" +set bbkg_mathpha "mathpha expr=\"${bbkg_expr}\" units=C outfil=${bbkg_outf} exposure=${expos_bkg} areascal=% ncomments=2 comment1=\"${bbkg_cmm1}\" comment2=\"${bbkg_cmm2}\" properr=no backscal=${bbkg_back}" +## blanksky }}} + +## localbkg {{{ +set lbkg_expr "${specscal}*${spec_name} -${fake_gas}" +set lbkg_outf "bkgcorr_localbkg_${rootname}.pi" +set lbkg_back [ expr { $backscal_src * $specscal } ] +set lbkg_cmm1 "corrected background spectrum, local background based" +set lbkg_cmm2 "norm*${specscal}, backscal*${specscal}, properr=no" +set lbkg_mathpha "mathpha expr=\"${lbkg_expr}\" units=C outfil=${lbkg_outf} exposure=${expos_src} areascal=% ncomments=2 comment1=\"${lbkg_cmm1}\" comment2=\"${lbkg_cmm2}\" properr=no backscal=${lbkg_back}" +## local }}} + +## XXX: cannot figure out how to `eval' these cmd in Tcl. !!! +## generate a shell script to exec mathpha {{{ +if {[ file exists ${bbkg_outf} ]} { + exec mv -fv ${bbkg_outf} ${bbkg_outf}_bak +} +if {[ file exists ${lbkg_outf} ]} { + exec mv -fv ${lbkg_outf} ${lbkg_outf}_bak +} +set mathpha_fn "_mathpha.sh" +if {[ file exists ${mathpha_fn} ]} { + exec rm -fv ${mathpha_fn} +} +set mathpha_fd [ open ${mathpha_fn} w ] +puts $mathpha_fd $bbkg_mathpha +puts $mathpha_fd $lbkg_mathpha +close $mathpha_fd +# exec +exec sh ${mathpha_fn} +#exec rm -rf ${mathpha_sh} +## mathpha }}} +## background correction }}} + +## load the corrected spectra, save a picture {{{ +set bkgcorr_img "bkgcorr_${rootname}.ps" +if {[ file exists ${bkgcorr_img} ]} { + exec mv -fv ${bkgcorr_img} ${bkgcorr_img}_bak +} + +model none +data none + +# corrected spectrum (blanksky based) +data 1:1 ${bbkg_outf} +resp 1:1 ${rmf} +arf 1:1 ${arf} +# corrected spectrum (local bkg based) +data 2:2 ${lbkg_outf} +resp 1:2 ${rmf} +arf 1:2 ${arf} +# original blanksky spectrum +data 3:3 ${back_name} +resp 1:3 ${rmf} +arf 1:3 ${arf} +# original local spectrum +data 4:4 ${spec_name} +resp 1:4 ${rmf} +arf 1:4 ${arf} + +ignore bad +ignore **:0.0-0.3,11.0-** + +setplot energy + +# rebin for plot +setplot rebin 20 20 1 +setplot rebin 20 20 2 +setplot rebin 20 20 3 +setplot rebin 10 20 4 + +# label +setplot command LABEL T "Corrected & Original Background Spectrum Comparison" +setplot command LABEL F "corr_blank(BLACK), corr_local(RED), orig_blank(GREEN), orig_local(BLUE)" + +setplot device ${bkgcorr_img}/cps +plot ldata +setplot device /xw +plot +## bkgcorr picture }}} + +## end analysis +# tclexit + +# vim: set ts=8 sw=4 tw=0 fenc=utf-8 ft=tcl: # diff --git a/xspec/xspec_coolfunc_v2.tcl b/xspec/xspec_coolfunc_v2.tcl new file mode 100644 index 0000000..0bfcbf5 --- /dev/null +++ b/xspec/xspec_coolfunc_v2.tcl @@ -0,0 +1,246 @@ +##################################################################### +## XSPEC Tcl script +## +## Task: +## To generate `Temperature profile' +## and `cooling function data file' +## according to `deprojection spectra analysis' results +## +## Requirements: +## deprojection spectra analysis, +## spectra loaded *in order* (in order to correctly calc *radius*) +## use XSPEC model `projct*wabs*apec' (AS IT IS) +## because this script cannot handle other *parameter numbers* +## +## NOTES: +## *solar abundance std* assumed to be `GRSA' (1998) +## if error < 0.01, then assume error = 0.01 +## if fabs(error) < 1e-7, then set error as `NULL' in output +## (this may be caused by frozen or tied params) +## +## Output: +## tcl_xspec_saveall.xcm (save all for the current state) +## tcl_fitted_params.log (log the fitted parameters) +## tcl_temp_profile.qdp (temp profile QDP file for graphing) +## tcl_temp_profile.txt (temp profile data file for later fit) +## tcl_coolfunc_dat.txt (coolfunc data file for later fit) +## +## LIweitiaNux <liweitianux@gmail.com> +## August 12, 2012 +## +## v1.1, 2012/08/12 +## improve description +## modify error treatment +## v2.0, 2012/08/14, LIweitiaNux +## add `output a QDP file for temperature profile' +## cancel `calc errors of abundance' +## improve error treatment (add warning) +##################################################################### + +## about {{{ +set NAME "xspec_coolfunc_v2.tcl" +set VERSION "v2, 2012-08-14" +## about }}} + +## ask `normalization' used to calculate cooling function +puts -nonewline "LY> Input the NORMALIZATION used to calc cooling function: " +set norm [ gets stdin ] + +## basic variables {{{ +## record process date +set DATE [ exec date ] +set ERR_FLAG "FALSE" +## basic vars }}} + +## xspec settings {{{ +# solar abundance std +# use `grsa' (1998) +set abund_std "grsa" + +# Return TCL results for XSPEC commands. +set xs_return_result 1 + +# Keep going until fit converges. +query yes +## xspec settings }}} + +## xspec, save current results {{{ +set fn_save "tcl_xspec_saveall.xcm" +if {[ file exists $fn_save ]} { + exec mv -fv $fn_save ${fn_save}_bak +} +save all $fn_save +## xspec save all }}} + +## xspec log, log current fitted params {{{ +set fn_log "tcl_fitted_params.log" +if {[ file exists $fn_log ]} { + exec mv -fv $fn_log ${fn_log}_bak +} +log $fn_log +show param +log none +## xspec log }}} + +## set output file {{{ +set tpro_qdp_fn "tcl_temp_profile.qdp" +set tpro_fn "tcl_temp_profile.txt" +set cool_fn "tcl_coolfunc_dat.txt" +# tmp file used to calc `cooling function' +set tmpxspec_fn "_tcl_xspec_tmp.xcm" +set tmpres_fn "_tcl_coolres_tmp.txt" +# check file status +if {[ file exists $tpro_qdp_fn ]} { + exec mv -fv $tpro_qdp_fn ${tpro_qdp_fn}_bak +} +if {[ file exists $tpro_fn ]} { + exec mv -fv $tpro_fn ${tpro_fn}_bak +} +if {[ file exists $cool_fn ]} { + exec mv -fv $cool_fn ${cool_fn}_bak +} +set tpro_qdp_fd [ open $tpro_qdp_fn w ] +set tpro_fd [ open $tpro_fn w ] +set cool_fd [ open $cool_fn w ] +## output files }}} + +## datasets, number of data group +tclout datasets +set datasets $xspec_tclout + +## get `nH' and `redshift' at first +tclout param 4 +scan $xspec_tclout "%f" nh +tclout param 7 +scan $xspec_tclout "%f" redshift + +## QDP file header {{{ +puts $tpro_qdp_fd "! generated by: ${NAME}" +puts $tpro_qdp_fd "! created date: ${DATE}" +puts $tpro_qdp_fd "SKIP SINGLE" +puts $tpro_qdp_fd "READ SERR 1 2" +puts $tpro_qdp_fd "" +puts $tpro_qdp_fd "LABEL T Temperature Profile" +puts $tpro_qdp_fd "LABEL F temperature error at 68% confidence level" +puts $tpro_qdp_fd "LABEL X Radius (pixel)" +puts $tpro_qdp_fd "LABEL Y Temperature (keV)" +puts $tpro_qdp_fd "" +puts $tpro_qdp_fd "! radius(pixel) radius_err temperature(keV) temp_err(60%)" +## QDP }}} + +## output temperature profile and abundance +for {set i 1} {$i <= ${datasets}} {incr i} { + # get the MID-radius of current spectra + tclout xflt $i + scan $xspec_tclout "%f %f" holder r_out + if {$i == 1} { + # the innermost region + set r_in 0.0 + } else { + set j [ expr {$i - 1} ] + tclout xflt $j + scan $xspec_tclout "%f %f" holder r_in + } + set r [ expr {($r_in + $r_out) / 2.0} ] + set r_err [ expr {($r_out - $r_in) / 2.0} ] + puts "radius: $r, radius_err: $r_err ($r_in, $r_out)" + # output values: temperature and abundance + # determine the param number of temperature and abundance + # temperature + set temp_pn [ expr {8 * $i - 3} ] + tclout param $temp_pn + scan $xspec_tclout "%f" temp_val + error 1.0 $temp_pn + tclout error $temp_pn + scan $xspec_tclout "%f %f" temp_val_l temp_val_u + set temp_err [ expr {($temp_val_u - $temp_val_l) / 2.0} ] + # error treatment + if {abs($temp_err) < 1.0e-7} { + puts "*** WARNING: WRONG error values" + set temp_err "NULL" + set ERR_FLAG "TRUE" + } elseif {$temp_err < 0.01} { + set temp_err 0.01 + } + puts "temperature (p${temp_pn}): $temp_val, $temp_err" + # abundance + set abund_pn [ expr {8 * $i - 2} ] + tclout param $abund_pn + scan $xspec_tclout "%f" abund_val + #error 1.0 $abund_pn + #tclout error $abund_pn + #scan $xspec_tclout "%f %f" abund_val_l abund_val_u + #set abund_err [ expr {($abund_val_u - $abund_val_l) / 2.0} ] + #if {abs($abund_err) < 1.0e-7} { + # puts "*** WARNING: WRONG error values" + # set abund_err "NULL" + # set ERR_FLAG "TRUE" + #} elseif {$abund_err < 0.01} { + # set abund_err 0.01 + #} + #puts "abundance (p${abund_pn}): $abund_val $abund_err" + # output `temp_profile' related results + puts $tpro_qdp_fd "$r $r_err $temp_val $temp_err" + puts $tpro_fd "$r $r_err $temp_val $temp_err" + + ## generate a xspec script to calc cooling function {{{ + # remove previous files + if {[ file exists $tmpxspec_fn ]} { + exec rm -fv $tmpxspec_fn + } + if {[ file exists $tmpres_fn ]} { + exec rm -fv $tmpres_fn + } + # open script file + set tmpxspec_fd [ open $tmpxspec_fn w ] + puts $tmpxspec_fd "set tmpres_fn \"$tmpres_fn\"" + puts $tmpxspec_fd "set tmpres_fd \[ open \$tmpres_fn w \]" + puts $tmpxspec_fd "query yes" + puts $tmpxspec_fd "statistic chi" + puts $tmpxspec_fd "method leven 10 0.01" + puts $tmpxspec_fd "abund ${abund_std}" + puts $tmpxspec_fd "xsect bcmc" + puts $tmpxspec_fd "cosmo 70 0 0.73" + puts $tmpxspec_fd "xset delta 0.01" + puts $tmpxspec_fd "systematic 0" + puts $tmpxspec_fd "dummyrsp 0.3 11.0 1024" + puts $tmpxspec_fd "model wabs*apec & ${nh} & ${temp_val} & ${abund_val} & ${redshift} & $norm & /*" + puts $tmpxspec_fd "flux 0.7 7.0" + puts $tmpxspec_fd "tclout flux 1" + puts $tmpxspec_fd "scan \$xspec_tclout \"%f %f %f %f\" holder holder holder res_flux" + puts $tmpxspec_fd "puts \$res_flux" + puts $tmpxspec_fd "puts \$tmpres_fd \$res_flux" + puts $tmpxspec_fd "close \$tmpres_fd" + puts $tmpxspec_fd "tclexit" + close $tmpxspec_fd + ## xspec script }}} + ## calc cooling function and output results {{{ + # exec XSPEC to calculate + exec xspec - $tmpxspec_fn + # read results from files + set tmpres_fd [ open $tmpres_fn r ] + set tmpline [ gets $tmpres_fd ] + scan $tmpline "%f" coolf_val + puts "cooling function: $coolf_val" + close $tmpres_fd + puts $cool_fd "$r $coolf_val" + ## cooling function }}} +} + +## remove tmp files +if {[ file exists $tmpres_fn ]} { + exec rm -fv $tmpres_fn +} + +## close file +close $tpro_qdp_fd +close $tpro_fd +close $cool_fd + +## check `ERR_FLAG', print WARNING +if {[ string equal $ERR_FLAG "TRUE" ]} { + puts "*** WARNING: there are WRONG error values" +} + +# EOF +# vim: set ts=8 sw=4 tw=0 fenc= ft=tcl: # |