diff options
Diffstat (limited to 'scripts/ciao_calc_ct.sh')
| -rwxr-xr-x | scripts/ciao_calc_ct.sh | 151 |
1 files changed, 94 insertions, 57 deletions
diff --git a/scripts/ciao_calc_ct.sh b/scripts/ciao_calc_ct.sh index 5f76ed5..58bf566 100755 --- a/scripts/ciao_calc_ct.sh +++ b/scripts/ciao_calc_ct.sh @@ -3,30 +3,34 @@ unalias -a export LC_COLLATE=C ########################################################### -## based on `ciao_r500avgt' ## -## for calculating the `cooling time ' ## +## based on `ciao_r500avgt.sh' ## +## for calculating the `cooling time' ## ## within (0.-0.048 r500) region ## ## ## -## Junhua Gu ## -## August 22, 2012 ## +## Author: Junhua GU ## +## Created: 2012/08/22 ## ## ## -## LIweitiaNux ## +## Modified by: Weitian LI ## ## 2013/04/28 ## ########################################################### - -########################################################### -## ChangeLogs -## v1.1, 2014/06/18 -## 'cooling_time2.sh' -> 'ciao_calc_ct.sh' -## v1.2, 2015/05/27, Weitian LI +## +VERSION="v2.0" +UPDATE="2015/06/03" +## +## Changelogs: +## v2.0, 2015/06/03, Aaron LI +## * Copy needed pfiles to current working directory, and +## set environment variable $PFILES to use these first. +## * Replace 'grep' with '\grep', 'ls' with '\ls' +## * replaced 'grppha' with 'dmgroup' to group spectra +## (dmgroup will add history to fits file, while grppha NOT) +## * ds9 colormap changed from 'sls' to 'he' +## v1.2, 2015/05/27, Aaron LI ## update 'DFT_ARF' & 'DFT_RMF' to find '*.arf' & '*.rmf' files ## (specextract only use .arf & .rmf extensions since revision 2014-12) -########################################################### - -## about, used in `usage' {{{ -VERSION="v1.2" -UPDATE="2015-05-27" -## about }}} +## v1.1, 2014/06/18 +## 'cooling_time2.sh' -> 'ciao_calc_ct.sh' +## ## error code {{{ ERR_USG=1 @@ -49,17 +53,15 @@ ERR_UNI=61 case "$1" in -[hH]*|--[hH]*) printf "usage:\n" - printf " `basename $0` basedir=<repro_dir> evt=<evt2_clean> r500=<r500_kpc> regin=<input_reg> regout=<output_reg> bkgd=<blank_evt | lbkg_reg | bkg_spec> nh=<nH> z=<redshift> arf=<arf_file> rmf=<rmf_file> [ grpcmd=<grppha_cmd> log=<log_file> ]\n" + printf " `basename $0` basedir=<repro_dir> evt=<evt2_clean> r500=<r500_kpc> regin=<input_reg> regout=<output_reg> bkgd=<blank_evt | lbkg_reg | bkg_spec> nh=<nH> z=<redshift> arf=<arf_file> rmf=<rmf_file> [ grouptype=<NUM_CTS|BIN> grouptypeval=<number> binspec=<binspec> log=<log_file> ]\n" printf "\nversion:\n" - printf "${VERSION}, ${UPDATE}\n" + printf " ${VERSION}, ${UPDATED}\n" exit ${ERR_USG} ;; esac ## usage, help }}} ## comology calculator {{{ -## XXX: MODIFY THIS TO YOUR OWN CASE -## and make sure this `calc' is executable ## NOTES: use `$HOME' instead of `~' in path BASE_PATH=`dirname $0` COSCALC="`which cosmo_calc calc_distance 2>/dev/null | head -n 1`" @@ -72,9 +74,9 @@ fi ## default parameters {{{ # default `event file' which used to match `blanksky' files #DFT_EVT="_NOT_EXIST_" -DFT_EVT="`ls evt2*_clean.fits 2> /dev/null`" +DFT_EVT="`\ls evt2*_clean.fits 2> /dev/null`" # default `bkgd', use `bkgcorr_blanksky*' corrected bkg spectrum -DFT_BKGD="`ls bkgcorr_blanksky_*.pi 2> /dev/null`" +DFT_BKGD="`\ls bkgcorr_blanksky_*.pi 2> /dev/null`" # default basedir DFT_BASEDIR="../.." # default info_json pattern @@ -85,14 +87,14 @@ DFT_REG_IN="rspec.reg" # default output region file (0.1-0.5 r500 region) DFT_REG_OUT="cooling.reg" # default ARF/RMF, the one of the outmost region -DFT_ARF="`ls r1_*.warf r1_*.arf 2> /dev/null`" -DFT_RMF="`ls r1_*.wrmf r1_*.rmf 2> /dev/null`" +DFT_ARF="`\ls r1_*.warf r1_*.arf 2> /dev/null`" +DFT_RMF="`\ls r1_*.wrmf r1_*.rmf 2> /dev/null`" -# default `group command' for `grppha' -#DFT_GRP_CMD="group 1 128 2 129 256 4 257 512 8 513 1024 16" -DFT_GRP_CMD="group min 20" -# default `log file' -DFT_LOGFILE="cooling_`date '+%Y%m%d'`.log" +# default parameters for 'dmgroup' +DFT_GROUPTYPE="NUM_CTS" +DFT_GROUPTYPEVAL="20" +#DFT_GROUPTYPE="BIN" +DFT_BINSPEC="1:128:2,129:256:4,257:512:8,513:1024:16" INNER=0.0 OUTER=0.048 @@ -102,6 +104,9 @@ XSPEC_SCRIPT="xspec_cooling.xcm" # deproj xspec script, generated by `deproj_spectra' # from which get `nh' and `redshift' XSPEC_DEPROJ="xspec_deproj.xcm" + +# default `log file' +DFT_LOGFILE="cooling_`date '+%Y%m%d'`.log" ## default parameters }}} ## functions {{{ @@ -126,7 +131,7 @@ CH_LOW=651 CH_HI=822 pb_flux() { punlearn dmstat - COUNTS=`dmstat "$1[channel=${CH_LOW}:${CH_HI}][cols COUNTS]" | grep -i 'sum:' | awk '{ print $2 }'` + COUNTS=`dmstat "$1[channel=${CH_LOW}:${CH_HI}][cols COUNTS]" | \grep -i 'sum:' | awk '{ print $2 }'` punlearn dmkeypar EXPTIME=`dmkeypar $1 EXPOSURE echo=yes` BACK=`dmkeypar $1 BACKSCAL echo=yes` @@ -179,7 +184,7 @@ if [ ! -r "${XSPEC_DEPROJ}" ]; then read -p "> value of redshift: " REDSHIFT else # get `nh' and `redshift' from xspec script - LN=`grep -n 'projct\*wabs\*apec' ${XSPEC_DEPROJ} | tail -n 1 | cut -d':' -f1` + LN=`\grep -n 'projct\*wabs\*apec' ${XSPEC_DEPROJ} | tail -n 1 | cut -d':' -f1` # calc the line number of which contains `nh' LN_NH=`expr ${LN} + 4` NH_XCM=`head -n ${LN_NH} ${XSPEC_DEPROJ} | tail -n 1 | awk '{ print $1 }'` @@ -215,8 +220,8 @@ fi # json file if [ ! -z "${json}" ] && [ -r "${BASEDIR}/${json}" ]; then JSON_FILE="${BASEDIR}/${json}" -elif [ `ls ${BASEDIR}/${DFT_JSON_PAT} 2> /dev/null | wc -l` -eq 1 ]; then - JSON_FILE=`ls ${BASEDIR}/${DFT_JSON_PAT}` +elif [ `\ls ${BASEDIR}/${DFT_JSON_PAT} 2> /dev/null | wc -l` -eq 1 ]; then + JSON_FILE=`\ls ${BASEDIR}/${DFT_JSON_PAT}` else read -p "> JSON_file: " JSON_FILE if [ ! -r "${JSON_FILE}" ]; then @@ -227,7 +232,7 @@ fi printf "## use json_file: \`${JSON_FILE}'\n" | ${TOLOG} # process `r500' {{{ -R500_RAW=`grep '"R500.*kpc' ${JSON_FILE} | sed 's/.*"R500.*":\ //' | sed 's/\ *,$//'` +R500_RAW=`\grep '"R500.*kpc' ${JSON_FILE} | sed 's/.*"R500.*":\ //' | sed 's/\ *,$//'` if [ ! -z "${r500}" ]; then R500_RAW=${r500} fi @@ -247,7 +252,7 @@ case "${R500_UNI}" in ;; *) printf "## units in \`kpc', convert to \`Chandra pixel'\n" | ${TOLOG} - KPC_PER_PIX=`${COSCALC} ${REDSHIFT} | grep 'kpc.*pix' | tr -d 'a-zA-Z_#=(),:/ '` + KPC_PER_PIX=`${COSCALC} ${REDSHIFT} | \grep 'kpc.*pix' | tr -d 'a-zA-Z_#=(),:/ '` # convert scientific notation for `bc' KPC_PER_PIX_B=`echo ${KPC_PER_PIX} | sed 's/[eE]/\*10\^/' | sed 's/+//'` printf "## calculated \`kpc/pixel': ${KPC_PER_PIX_B}\n" @@ -299,7 +304,7 @@ printf "## set output regfile (0.1-0.5 r500 region): \`${REG_OUT}'\n" | ${TOLOG} # get center position from `regin' # only consider `pie' or `annulus'-shaped region -TMP_REG=`grep -iE '(pie|annulus)' ${REG_IN} | head -n 1` +TMP_REG=`\grep -iE '(pie|annulus)' ${REG_IN} | head -n 1` XC=`echo ${TMP_REG} | tr -d ' ' | awk -F'[(),]' '{ print $2 }'` YC=`echo ${TMP_REG} | tr -d ' ' | awk -F'[(),]' '{ print $3 }'` printf "## get center coord: (${XC},${YC})\n" | ${TOLOG} @@ -321,7 +326,7 @@ printf "## use bkgd: \`${BKGD}'\n" | ${TOLOG} # determine bkg type: blanksky, lbkg_reg, bkg_spec ? # according to file type first: text / FITS # if FITS, then get values of `HDUCLAS1' and `OBJECT' -if file -bL ${BKGD} | grep -qi 'text'; then +if file -bL ${BKGD} | \grep -qi 'text'; then printf "## given \`${BKGD}' is a \`text file'\n" printf "## use it as local bkg region file\n" printf "## use *LOCAL BKG SPEC*\n" | ${TOLOG} @@ -329,9 +334,10 @@ if file -bL ${BKGD} | grep -qi 'text'; then USE_LBKG_REG=YES USE_BLANKSKY=NO USE_BKG_SPEC=NO -elif file -bL ${BKGD} | grep -qi 'FITS'; then +elif file -bL ${BKGD} | \grep -qi 'FITS'; then printf "## given \`${BKGD}' is a \`FITS file'\n" # get FITS header keyword + punlearn dmkeypar HDUCLAS1=`dmkeypar ${BKGD} HDUCLAS1 echo=yes` if [ "${HDUCLAS1}" = "EVENTS" ]; then # event file @@ -400,21 +406,49 @@ fi printf "## use RMF: \`${RMF}'\n" | ${TOLOG} # arf & rmf }}} -# check given `grpcmd' -if [ ! -z "${grpcmd}" ]; then - GRP_CMD="${grpcmd}" +# check given dmgroup parameters: grouptype, grouptypeval, binspec +if [ -z "${grouptype}" ]; then + GROUPTYPE="${DFT_GROUPTYPE}" +elif [ "x${grouptype}" = "xNUM_CTS" ] || [ "x${grouptype}" = "xBIN" ]; then + GROUPTYPE="${grouptype}" +else + printf "ERROR: given grouptype \`${grouptype}' invalid.\n" + exit ${ERR_GRPTYPE} +fi +printf "## use grouptype: \`${GROUPTYPE}'\n" | ${TOLOG} +if [ ! -z "${grouptypeval}" ]; then + GROUPTYPEVAL="${grouptypeval}" else - GRP_CMD="${DFT_GRP_CMD}" + GROUPTYPEVAL="${DFT_GROUPTYPEVAL}" fi -printf "## use grppha cmd: \`${GRP_CMD}'\n" | ${TOLOG} +printf "## use grouptypeval: \`${GROUPTYPEVAL}'\n" | ${TOLOG} +if [ ! -z "${binspec}" ]; then + BINSPEC="${binspec}" +else + BINSPEC="${DFT_BINSPEC}" +fi +printf "## use binspec: \`${BINSPEC}'\n" | ${TOLOG} ## parameters }}} -################################################## -#### main -## D_A -#D_A_CM=`${COSCALC} ${REDSHIFT} | grep '^d_a_cm' | awk '{ print $2 }'` -D_A_CM=`${COSCALC} ${REDSHIFT} | grep -i 'd_a.*cm' | awk -F'=' '{ print $2 }' | awk '{ print $1 }'` +## prepare parameter files (pfiles) {{{ +CIAO_TOOLS="dmstat dmkeypar dmhedit dmextract dmlist dmgroup" + +# Copy necessary pfiles for localized usage +for tool in ${CIAO_TOOLS}; do + pfile=`paccess ${tool}` + [ -n "${pfile}" ] && punlearn ${tool} && cp -Lvf ${pfile} . +done + +# Modify environment variable 'PFILES' to use local pfiles first +export PFILES="./:${PFILES}" +## pfiles }}} + +### main ### + +## D_A {{{ +D_A_CM=`${COSCALC} ${REDSHIFT} | \grep -i 'd_a.*cm' | awk -F'=' '{ print $2 }' | awk '{ print $1 }'` printf "D_A_CM(${REDSHIFT})= ${D_A_CM}\n" +## D_A }}} ## region related {{{ ## generate the needed region file @@ -427,11 +461,11 @@ _EOF_ ## open the evt file to verify or modify printf "## check the generated pie region ...\n" printf "## if modified, save with the same name \`${REG_OUT}' (overwrite)\n" -ds9 ${EVT} -regions ${REG_OUT} -cmap sls -bin factor 4 +ds9 ${EVT} -regions ${REG_OUT} -cmap he -bin factor 4 ## check the (modified) region (pie region end angle) printf "check the above region (for pie region end angle) ...\n" -INVALID=`grep -i 'pie' ${REG_OUT} | awk -F'[,()]' '$7 > 360'` +INVALID=`\grep -i 'pie' ${REG_OUT} | awk -F'[,()]' '$7 > 360'` if [ "x${INVALID}" != "x" ]; then printf "*** WARNING: there are pie regions' END_ANGLE > 360\n" | ${TOLOG} printf "*** will to fix ...\n" @@ -452,7 +486,7 @@ fi ## generate spectrum {{{ # check counts punlearn dmlist -CNT_RC=`dmlist infile="${EVT}[sky=region(${REG_OUT})][energy=700:7000]" opt=block | grep 'EVENTS' | awk '{ print $8 }'` +CNT_RC=`dmlist infile="${EVT}[sky=region(${REG_OUT})][energy=700:7000]" opt=block | \grep 'EVENTS' | awk '{ print $8 }'` if [ ${CNT_RC} -lt 500 ]; then F_WC=true WC="LOW_COUNTS" @@ -468,8 +502,11 @@ dmextract infile="${EVT}[sky=region(${REG_OUT})][bin PI]" \ outfile="${AVGT_SPEC}" wmap="[bin det=8]" clobber=yes # group spectrum printf "group object spectrum ...\n" -grppha infile="${AVGT_SPEC}" outfile="${AVGT_SPEC_GRP}" \ - comm="${GRP_CMD} & exit" clobber=yes > /dev/null +punlearn dmgroup +dmgroup infile="${AVGT_SPEC}" outfile="${AVGT_SPEC_GRP}" \ + grouptype="${GROUPTYPE}" grouptypeval=${GROUPTYPEVAL} \ + binspec="${BINSPEC}" xcolumn="CHANNEL" ycolumn="COUNTS" \ + clobber=yes # background printf "generate the background spectrum ...\n" @@ -585,9 +622,9 @@ else printf "invoking XSPEC to calculate cooling time ...\n" xspec - ${XSPEC_SCRIPT} | tee ${CT_RES} - OBS_ID=`grep '"Obs.*ID' ${JSON_FILE} | awk -F':' '{ print $2 }' | tr -d ' ,'` - OBJ_NAME=`grep '"Source\ Name' ${JSON_FILE} | awk -F':' '{ print $2 }' | sed -e 's/\ *"//' -e 's/"\ *,$//'` - CT=`grep -i '^Cooling_time' ${CT_RES} | awk '{ print $2 }'` + OBS_ID=`\grep '"Obs.*ID' ${JSON_FILE} | awk -F':' '{ print $2 }' | tr -d ' ,'` + OBJ_NAME=`\grep '"Source\ Name' ${JSON_FILE} | awk -F':' '{ print $2 }' | sed -e 's/\ *"//' -e 's/"\ *,$//'` + CT=`\grep -i '^Cooling_time' ${CT_RES} | awk '{ print $2 }'` printf "\n" | tee -a ${CT_RES} printf "# OBS_ID,OBJ_NAME,CT_gyr\n" | tee -a ${CT_RES} |