chandra_update_xpeak2.sh: add to find X-peak in the cleaned image

1 files changed, 311 insertions, 0 deletions

diff --git a/scripts/chandra_update_xpeak2.sh b/scripts/chandra_update_xpeak2.sh
new file mode 100755
index 0000000..d8cab0c
--- /dev/null
+++ b/scripts/chandra_update_xpeak2.sh
@@ -0,0 +1,311 @@
+#!/bin/sh
+#
+unalias -a
+export LC_COLLATE=C
+##
+## Take the `X-ray centroid' coordinate from region `sbprofile.reg'
+## as the start location, then search for X-ray peak in a circle region
+## of radius ${SEARCH_RADIUS} pixels centered on the X-ray centroid.
+## After that the found X-ray peak is convert from physical coordinate to
+## WCS coordinate and added/updated to the INFO json file.
+##
+## Note: The image on which the X-ray peak is to be searched is generated
+## as following:
+##     (1) cleaned event file (c7, c0-3)
+##     (2) filter energy (${E_RANGE})
+##     (3) smooth with aconvolve
+##
+## Based on `chandra_update_xpeak.sh'
+##
+## Weitian LI <liweitianux@gmail.com>
+## Created: 2015-11-22
+##
+VERSION="v1.0"
+UPDATED="2015-11-20"
+##
+## ChangeLogs:
+##
+
+## error code {{{
+ERR_USG=1
+ERR_DIR=11
+ERR_EVT=12
+ERR_BKG=13
+ERR_REG=14
+ERR_INFO=15
+ERR_ASOL=21
+ERR_BPIX=22
+ERR_PBK=23
+ERR_MSK=24
+ERR_BKGTY=31
+ERR_SPEC=32
+ERR_DET=41
+ERR_ENG=42
+ERR_CIAO=100
+## error code }}}
+
+## usage, help {{{
+case "$1" in
+    -[hH]*|--[hH]*)
+        printf "usage:\n"
+        printf "    `basename $0` evt=<evt_file> reg=<sbp_reg> basedir=<base_dir> info=<INFO.json> conv=<YES|no> update=<YES|no>\n"
+        printf "\nversion:\n"
+        printf "${VERSION}, ${UPDATED}\n"
+        exit ${ERR_USG}
+        ;;
+esac
+## usage, help }}}
+
+## default parameters {{{
+# default circle radius within which to search the X-ray peak
+SEARCH_RADIUS=100
+# critical offset (in pixel)
+OFFSET_CRIC=20
+
+# energy range: 700-2000 eV
+E_RANGE="700:2000"
+# default `event file' which used to match `blanksky' files
+DFT_EVT="`\ls evt2*c7*_clean.fits evt2*c0-3*_clean.fits 2> /dev/null`"
+# default dir which contains `asols, asol.lis, ...' files
+DFT_BASEDIR=".."
+# default `radial region file' to extract surface brightness
+DFT_SBP_REG="`\ls sbprofile.reg rspec.reg 2> /dev/null | head -n 1`"
+
+## howto find files in `basedir'
+# default `asol.lis pattern'
+DFT_ASOLIS_PAT="acis*asol?.lis"
+# default INFO.json pattern
+DFT_INFO_PAT="*_INFO.json"
+
+# aconvolve parameters
+ACONV_KERNELSPEC="lib:gaus(2,3,1,3,3)"
+ACONV_METHOD="fft"
+## default parameters }}}
+
+## functions {{{
+# process commandline arguments
+# cmdline arg format: `KEY=VALUE'
+getopt_keyval() {
+    until [ -z "$1" ]
+    do
+        key=${1%%=*}                    # extract key
+        val=${1#*=}                     # extract value
+        keyval="${key}=\"${val}\""
+        echo "## getopt: eval '${keyval}'"
+        eval ${keyval}
+        shift                           # shift, process next one
+    done
+}
+## functions }}}
+
+## check CIAO init {{{
+if [ -z "${ASCDS_INSTALL}" ]; then
+    printf "ERROR: CIAO NOT initialized\n"
+    exit ${ERR_CIAO}
+fi
+
+## XXX: heasoft's `pget' etc. tools conflict with some CIAO tools
+printf "set \$PATH to avoid conflicts between HEAsoft and CIAO\n"
+export PATH="${ASCDS_BIN}:${ASCDS_CONTRIB}:${PATH}"
+## check CIAO }}}
+
+## parameters {{{
+# process cmdline args using `getopt_keyval'
+getopt_keyval "$@"
+
+# check given parameters
+# check evt file
+if [ -r "${evt}" ]; then
+    EVT=${evt}
+elif [ -r "${DFT_EVT}" ]; then
+    EVT=${DFT_EVT}
+else
+    read -p "clean evt2 file: " EVT
+    if [ ! -r "${EVT}" ]; then
+        printf "ERROR: cannot access given \`${EVT}' evt file\n"
+        exit ${ERR_EVT}
+    fi
+fi
+printf "## use evt file: \`${EVT}'\n" | ${TOLOG}
+# check given region file(s)
+if [ -r "${reg}" ]; then
+    SBP_REG="${reg}"
+elif [ -r "${DFT_SBP_REG}" ]; then
+    SBP_REG=${DFT_SBP_REG}
+else
+    read -p "> surface brighness radial region file: " SBP_REG
+    if [ ! -r "${SBP_REG}" ]; then
+        printf "ERROR: cannot access given \`${SBP_REG}' region file\n"
+        exit ${ERR_REG}
+    fi
+fi
+printf "## use reg file(s): \`${SBP_REG}'\n" | ${TOLOG}
+# check given dir
+if [ -d "${basedir}" ]; then
+    BASEDIR=${basedir}
+elif [ -d "${DFT_BASEDIR}" ]; then
+    BASEDIR=${DFT_BASEDIR}
+else
+    read -p "> basedir (contains asol files): " BASEDIR
+    if [ ! -d "${BASEDIR}" ]; then
+        printf "ERROR: given \`${BASEDIR}' NOT a directory\n"
+        exit ${ERR_DIR}
+    fi
+fi
+# remove the trailing '/'
+BASEDIR=`echo ${BASEDIR} | sed 's/\/*$//'`
+printf "## use basedir: \`${BASEDIR}'\n" | ${TOLOG}
+# check INFO.json file
+if [ ! -z "${info}" ] && [ -r "${BASEDIR}/${info}" ]; then
+    INFO_JSON="${info}"
+elif [ "`\ls ${BASEDIR}/${DFT_INFO_PAT} | wc -l`" -eq 1 ]; then
+    INFO_JSON=`( cd ${BASEDIR} && \ls ${DFT_INFO_PAT} )`
+else
+    read -p "> info json file: " INFO_JSON
+    if ! [ -r "${BASEDIR}/${INFO_JSON}" ]; then
+        printf "ERROR: cannot access given \`${BASEDIR}/${INFO_JSON}' file\n"
+        exit ${ERR_INFO}
+    fi
+fi
+INFO_JSON=`readlink -f ${BASEDIR}/${INFO_JSON}`
+printf "## use info json file: \`${INFO_JSON}'\n"
+# update flag: whether to update xcentroid in the info.json file
+if [ ! -z "${update}" ]; then
+    case "${update}" in
+        [nN][oO]|[fF]*)
+            F_UPDATE="NO"
+            ;;
+        *)
+            F_UPDATE="YES"
+            ;;
+    esac
+else
+    F_UPDATE="YES"
+fi
+
+# convolve (optional)
+if [ -z "${conv}" ]; then
+    CONV="YES"
+else
+    case "${conv}" in
+        [nN]*)
+            CONV="NO"
+            printf "## Do NOT apply \`aconvolve' !\n"
+            ;;
+        *)
+            CONV="YES"
+            printf "## apply \`aconvolve'\n"
+            ;;
+    esac
+fi
+## parameters }}}
+
+## prepare parameter files (pfiles) {{{
+CIAO_TOOLS="dmcopy dmstat aconvolve dmcoords"
+
+# 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 process {{{
+# Use previously generated `skyfov'
+SKYFOV=`\ls *skyfov*.fits 2>/dev/null | head -n 1`
+
+# Extract chip(s) from evt filename
+CHIP=`echo "${EVT}" | sed 's/^.*_c\(7\|0-3\)_.*$/\1/' | tr '-' ':'`
+
+# generate image
+IMG="img_c`echo ${CHIP} | tr ':' '-'`_e`echo ${E_RANGE} | tr ':' '-'`.fits"
+if [ -r "${IMG}" ]; then
+    printf "use previously generated image: \`${IMG}'\n"
+else
+    printf "generate image: \`${IMG}' ...\n"
+    punlearn dmcopy
+    dmcopy infile="${EVT_DEFLARE}[sky=region(${SKYFOV}[ccd_id=${CHIP}])][energy=${E_RANGE}][bin sky=::1]" outfile="${IMG}" clobber=yes
+fi
+
+# aconvolve
+if [ "${CONV}" = "YES" ]; then
+    IMG_ACONV="${IMG%.fits}_aconv.fits"
+    #if [ -r "${IMG_ACONV}" ]; then
+    #    printf "use previously convolved image: \`${IMG_ACONV}'\n"
+    #else
+        printf "\`aconvolve' to smooth img: \`${IMG_ACONV}' ...\n"
+        printf "## aconvolve: kernelspec=\"${ACONV_KERNELSPEC}\" method=\"${ACONV_METHOD}\"\n"
+        punlearn aconvolve
+        aconvolve infile="${IMG}" outfile="${IMG_ACONV}" kernelspec="${ACONV_KERNELSPEC}" method="${ACONV_METHOD}" clobber=yes
+    #fi
+else
+    IMG_ACONV=${IMG}
+fi
+
+# Get X-ray centroid coordinate from sbp region
+printf "get (x,y) from ${SBP_REG}\n"
+CNTRD_X=`\grep -iE '(pie|annulus)' ${SBP_REG} | head -n 1 | awk -F',' '{ print $1 }' | tr -d 'a-zA-Z() '`
+CNTRD_Y=`\grep -iE '(pie|annulus)' ${SBP_REG} | head -n 1 | awk -F',' '{ print $2 }' | tr -d 'a-zA-Z() '`
+
+SEARCH_REGION="circle(${CNTRD_X},${CNTRD_Y},${SEARCH_RADIUS})"
+# dmstat to find the maximum location
+printf "\`dmstat' to find the peak ...\n"
+punlearn dmstat
+dmstat infile="${IMG_ACONV}[sky=${SEARCH_REGION}]" verbose=0
+PEAK_X=`pget dmstat out_max_loc | cut -d',' -f1`
+PEAK_Y=`pget dmstat out_max_loc | cut -d',' -f2`
+
+# asolis
+ASOLIS=`( cd ${BASEDIR} && \ls ${DFT_ASOLIS_PAT} 2> /dev/null )`
+
+# Use "dmcoords" to convert (x,y) to (ra,dec)
+printf "\`dmcoords' to convert (x,y) to (ra,dec) ...\n"
+punlearn dmcoords
+dmcoords infile="${EVT}" asolfile="@${BASEDIR}/${ASOLIS}" option=sky x=${PEAK_X} y=${PEAK_Y}
+PEAK_RA=`pget dmcoords ra`
+PEAK_DEC=`pget dmcoords dec`
+
+# Calculate the offset between peak and centroid coordinate
+OFFSET=`echo "scale=5; sqrt((${PEAK_X}-${CNTRD_X})^2 + (${PEAK_Y}-${CNTRD_Y})^2)" | bc -l`
+
+printf "## X-ray centroid (x,y): (${CNTRD_X},${CNTRD_Y})\n"
+printf "## X-ray peak (x,y): (${PEAK_X},${PEAK_Y})\n"
+printf "## X-ray peak (ra,dec): (${PEAK_RA},${PEAK_DEC})\n"
+printf "## Offset (pixel): ${OFFSET}\n"
+if [ `echo "${OFFSET} > ${OFFSET_CRIC}" | bc -l` -eq 1 ]; then
+    printf "*** WARNING: large offset (> ${OFFSET_CRIC}) ***\n"
+fi
+
+# Output X-ray peak coordinate to a region file
+PEAK_PHY_REG="peak2_phy.reg"
+[ -e "${PEAK_PHY_REG}" ] && mv -f ${PEAK_PHY_REG} ${PEAK_PHY_REG}_bak
+echo "point(${PEAK_X},${PEAK_Y})" > ${PEAK_PHY_REG}
+PEAK_WCS_REG="peak2_wcs.reg"
+[ -e "${PEAK_WCS_REG}" ] && mv -f ${PEAK_WCS_REG} ${PEAK_WCS_REG}_bak
+echo "point(${PEAK_RA},${PEAK_DEC})" > ${PEAK_WCS_REG}
+
+if [ "${F_UPDATE}" = "YES" ]; then
+    cp -f ${INFO_JSON} ${INFO_JSON}_bak
+    printf "update/add X-ray peak coordinate to info.json ...\n"
+    if \grep -qE 'XPEAK2_(RA|DEC)' ${INFO_JSON}; then
+        printf "update ...\n"
+        sed -i'' "s/XPEAK2_RA.*$/XPEAK2_RA\":\ \"${PEAK_RA}\",/" ${INFO_JSON}
+        sed -i'' "s/XPEAK2_DEC.*$/XPEAK2_DEC\":\ \"${PEAK_DEC}\",/" ${INFO_JSON}
+        sed -i'' "s/XPEAK2_XCNTRD_dist.*$/XPEAK2_XCNTRD_dist\ (pix)\":\ \"${OFFSET}\",/" ${INFO_JSON}
+    else
+        printf "add ...\n"
+        sed -i'' "/\"Dec\.\"/ a\
+\ \ \ \ \"XPEAK2_XCNTRD_dist\ (pix)\": \"${OFFSET}\"," ${INFO_JSON}
+        sed -i'' "/\"Dec\.\"/ a\
+\ \ \ \ \"XPEAK2_DEC\": \"${PEAK_DEC}\"," ${INFO_JSON}
+        sed -i'' "/\"Dec\.\"/ a\
+\ \ \ \ \"XPEAK2_RA\": \"${PEAK_RA}\"," ${INFO_JSON}
+    fi
+fi
+## main }}}
+
+exit 0
+