diff options
| -rwxr-xr-x | scripts/chandra_gensbpreg.sh | 173 |
1 files changed, 73 insertions, 100 deletions
diff --git a/scripts/chandra_gensbpreg.sh b/scripts/chandra_gensbpreg.sh index 9a2014e..5fc1a60 100755 --- a/scripts/chandra_gensbpreg.sh +++ b/scripts/chandra_gensbpreg.sh @@ -1,38 +1,43 @@ #!/bin/sh -## -## This script generate a series of regions for the extraction of -## radial surface brightness profile (SBP). -## -## Regions geneartion algorithm: -## (1) innermost 10 regions, we require a mininal of 5 pixel as well -## as 50 counts within the 0.7-7.0 keV range. -## (2) following regions: R_out = R_in * 1.2, and require STN > 1.5. -## -## Reference: -## [1] region generate algorithm ??? (TODO) -## -## Author: Zhenghao ZHU -## Created: ??? (TODO) -## -UPDATED="2015/06/03" -## ChangeLogs: -## v3.0, 2015/06/03, Aaron LI -## * Copy needed pfiles to current working directory, and -## set environment variable $PFILES to use these first. -## * Added missing punlearn -## * Removed the section of dmlist.par & dmextract.par deletion -## v2.1, 2015/02/13, Weitian LI -## * added '1' to denominators when calculate STN to avoid division by zero -## * added script description -## v2.0, 2013/03/06, Weitian LI -## * added the parameter `-t' to print STN results for testing -## +# +# This script generate a series of regions for the extraction of +# radial surface brightness profile (SBP). +# +# Regions geneartion algorithm: +# (1) innermost 10 regions, we require a mininal of 5 pixel as well +# as 50 counts within the 0.7-7.0 keV range. +# (2) following regions: R_out = R_in * 1.2, and require SNR > 1.5. +# SNR = ??? (TODO) +# +# Reference: +# [1] region generate algorithm ??? (TODO) +# +# Author: Zhenghao ZHU +# Created: ??? (TODO) +# +# Change logs: +# 2017-02-26, Weitian LI +# * Further simplify arguments handling +# * Remove test support (-t) for simplification +# * Rename 'stn' to 'SNR' +# * Add ds9 view +# v3.0, 2015/06/03, Weitian LI +# * Copy needed pfiles to current working directory, and +# set environment variable $PFILES to use these first. +# * Added missing punlearn +# * Removed the section of dmlist.par & dmextract.par deletion +# v2.1, 2015/02/13, Weitian LI +# * added '1' to denominators when calculate SNR to avoid division by zero +# * added script description +# v2.0, 2013/03/06, Weitian LI +# * added the parameter `-t' to print SNR results for testing # minimal counts CNT_MIN=50 # energy: 700-7000eV -- channel 49:479 +ERANGE=700:7000 CH_LOW=49 CH_HI=479 @@ -40,6 +45,26 @@ CH_HI=479 CH_BKG_LOW=651 CH_BKG_HI=822 +if [ $# -ne 4 ]; then + printf "usage:\n" + printf " `basename $0` <evt> <bkg_pi> <reg_in> <reg_out>\n" + exit 1 +fi + +EVT=$1 +BKGSPC=$2 +REG_IN=$3 +REG_OUT=$4 + +X=`\grep -i 'point' ${REG_IN} | head -n 1 | tr -d 'a-zA-Z() ' | awk -F',' '{ print $1 }'` +Y=`\grep -i 'point' ${REG_IN} | head -n 1 | tr -d 'a-zA-Z() ' | awk -F',' '{ print $2 }'` + +echo "EVT: ${EVT}" +echo "ERANGE: ${ERANGE}" +echo "Center: (${X},${Y})" +echo "BKGSPC: ${BKGSPC}" +echo "" + ## prepare parameter files (pfiles) {{{ CIAO_TOOLS="dmstat dmlist dmextract" @@ -54,97 +79,43 @@ export PFILES="./:${PFILES}" ## pfiles }}} -if [ $# -lt 6 ]; then - printf "usage:\n" - printf " `basename $0` <evt> <evt_e> <x> <y> <bkg_pi> <reg_out>\n" - printf " `basename $0` -t <evt> <evt_e> <x> <y> <bkg_pi> <rin1> ...\n" - printf "updated:\n" - printf " ${UPDATED}\n" - exit 1 -fi - -if [ "x$1" = "x-t" ] && [ $# -ge 7 ]; then - EVT=$2 - EVT_E=$3 - X=$4 - Y=$5 - BKGSPC=$6 - # process <rin> ... - printf "REG -- STN\n" - while [ ! -z "$7" ]; do - RIN=$7 - shift - ROUT=`echo "scale=2; $RIN * 1.2" | bc -l` - TMP_REG="pie($X,$Y,$RIN,$ROUT,0,360)" - TMP_SPC="_tmpspc.pi" - punlearn dmextract - dmextract infile="${EVT}[sky=${TMP_REG}][bin pi]" outfile=${TMP_SPC} wmap="[energy=300:12000][bin det=8]" clobber=yes - punlearn dmstat - INDEX_SRC=`dmstat "${TMP_SPC}[channel=${CH_BKG_LOW}:${CH_BKG_HI}][cols counts]" | \grep "sum:" | awk '{print $2}'` - INDEX_BKG=`dmstat "${BKGSPC}[channel=${CH_BKG_LOW}:${CH_BKG_HI}][cols counts]" | \grep "sum:" | awk '{print $2}'` - COUNT_SRC=`dmstat "${TMP_SPC}[channel=${CH_LOW}:${CH_HI}][cols counts]" | \grep "sum:" | awk '{print $2}'` - COUNT_BKG=`dmstat "${BKGSPC}[channel=${CH_LOW}:${CH_HI}][cols counts]" | \grep "sum:" | awk '{print $2}'` - # echo "CNT_SRC: ${COUNT_SRC}, IDX_SRC: ${INDEX_SRC}, CNT_BKG: ${COUNT_BKG}, IDX_BKG: ${INDEX_BKG}" - # exit - STN=`echo ${COUNT_SRC} ${INDEX_SRC} ${COUNT_BKG} ${INDEX_BKG} | awk '{ printf("%f",$1/$2/$3*$4) }'` - printf "${TMP_REG} -- ${STN}\n" - [ -e "${TMP_SPC}" ] && rm -f ${TMP_SPC} - done - # exit - exit 0 -fi - -EVT=$1 -EVT_E=$2 -X=$3 -Y=$4 -BKGSPC=$5 -REG_OUT=$6 - [ -f "${REG_OUT}" ] && mv -fv ${REG_OUT} ${REG_OUT}_bak -echo "EVT: ${EVT}" -echo "EVT_E: ${EVT_E}" -echo "X: ${X}" -echo "Y: ${Y}" -echo "BKGSPC: ${BKGSPC}" -echo "" - RIN=0 ROUT=0 CNTS=0 for i in `seq 1 10`; do - printf "gen reg #$i @ cnts:${CNT_MIN} ...\n" + printf "Generate region #$i @ cnts:${CNT_MIN} ...\n" ROUT=`expr $RIN + 5` TMP_REG="pie($X,$Y,$RIN,$ROUT,0,360)" punlearn dmlist - CNTS=`dmlist "${EVT_E}[sky=${TMP_REG}]" blocks | \grep 'EVENTS' | awk '{ print $8 }'` + CNTS=`dmlist "${EVT}[energy=${ERANGE}][sky=${TMP_REG}]" blocks | \grep 'EVENTS' | awk '{ print $8 }'` while [ `echo "$CNTS < $CNT_MIN" | bc -l` -eq 1 ]; do ROUT=`expr $ROUT + 1` TMP_REG="pie($X,$Y,$RIN,$ROUT,0,360)" punlearn dmlist - CNTS=`dmlist "${EVT_E}[sky=${TMP_REG}]" blocks | \grep 'EVENTS' | awk '{ print $8 }'` + CNTS=`dmlist "${EVT}[energy=${ERANGE}][sky=${TMP_REG}]" blocks | \grep 'EVENTS' | awk '{ print $8 }'` done TMP_REG="pie($X,$Y,$RIN,$ROUT,0,360)" echo "${TMP_REG}" >> ${REG_OUT} RIN=$ROUT done -# last reg width +# last region width #REG_W=`tail -n 1 ${REG_OUT} | tr -d '()' | awk -F',' '{ print $4-$3 }'` RIN=$ROUT ROUT=`echo "scale=2; $ROUT * 1.2" | bc -l` -STN=10 +SNR=10 i=11 -STN_FILE="sbp_stn.dat" -[ -e ${STN_FILE} ] && mv -fv ${STN_FILE} ${STN_FILE}_bak -while [ `echo "${STN} > 1.5" | bc -l` -eq 1 ]; do - printf "gen reg #$i \n" +SNR_FILE="sbp_snr.dat" +[ -e ${SNR_FILE} ] && mv -fv ${SNR_FILE} ${SNR_FILE}_bak +while [ `echo "${SNR} > 1.5" | bc -l` -eq 1 ]; do + printf "Generate SBP region #$i \n" i=`expr $i + 1` TMP_REG="pie($X,$Y,$RIN,$ROUT,0,360)" - - # next reg + + # next region RIN=$ROUT ROUT=`echo "scale=2; $ROUT * 1.2" | bc -l` TMP_SPC=_tmpspc.pi @@ -158,20 +129,22 @@ while [ `echo "${STN} > 1.5" | bc -l` -eq 1 ]; do COUNT_BKG=`dmstat "${BKGSPC}[channel=${CH_LOW}:${CH_HI}][cols counts]" | \grep "sum:" | awk '{print $2}'` #echo "CNT_SRC: ${COUNT_SRC}, IDX_SRC: ${INDEX_SRC}, CNT_BKG: ${COUNT_BKG}, IDX_BKG: ${INDEX_BKG}" - #exit 99 # Add '1' to the denominators to avoid division by zero. - STN=`echo ${COUNT_SRC} ${INDEX_SRC} ${COUNT_BKG} ${INDEX_BKG} | awk '{ printf("%f", ($1 / ($2 + 1)) / ($3 / ($4 + 1))) }'` + SNR=`echo ${COUNT_SRC} ${INDEX_SRC} ${COUNT_BKG} ${INDEX_BKG} | awk '{ printf("%f", ($1 / ($2 + 1)) / ($3 / ($4 + 1))) }'` punlearn dmlist - CNT=`dmlist "${EVT_E}[sky=${TMP_REG}]" blocks | \grep 'EVENTS' | awk '{ print $8 }'` + CNT=`dmlist "${EVT}[energy=${ERANGE}][sky=${TMP_REG}]" blocks | \grep 'EVENTS' | awk '{ print $8 }'` echo "CNT: ${CNT}" echo "CNT_MIN: ${CNT_MIN}" - if [ `echo "${CNT} < ${CNT_MIN}" | bc -l` -eq 1 ]; then + if [ `echo "${CNT} < ${CNT_MIN}" | bc -l` -eq 1 ]; then break fi - echo "STN: ${STN}" - echo "RIN: ${RIN}, ROUT: ${ROUT}" - echo "STN: ${STN}" >> ${STN_FILE} + echo "SNR: ${SNR}" + echo "RIN: ${RIN}, ROUT: ${ROUT}" + echo "SNR: ${SNR}" >> ${SNR_FILE} echo "${TMP_REG}" >> ${REG_OUT} done +printf "check SBP regions ...\n" +ds9 ${EVT} -regions format ciao -regions system physical \ + -regions ${REG_OUT} -cmap he -bin factor 4 |