fix to 'ciao_expcorr.sh' & removed versions in filenames

1 files changed, 388 insertions, 0 deletions

diff --git a/scripts/ciao_expcorr.sh b/scripts/ciao_expcorr.sh
new file mode 100755
index 0000000..701ba01
--- /dev/null
+++ b/scripts/ciao_expcorr.sh
@@ -0,0 +1,388 @@
+#!/bin/sh
+#
+unalias -a
+export LC_COLLATE=C
+###########################################################
+## make `image' from `evt file'                          ##
+## make `spectral weight' using `make_instmap_weights'   ##
+## use `fluximage' to generating `exposure map'          ##
+## make `exposure-corrected' image                       ##
+## and extract `surface brighness profile'               ##
+##                                                       ##
+## NOTES:                                                ##
+## only ACIS-I (chip: 0-3) and ACIS-S (chip: 7) supported##
+## `merge_all' conflict with Heasoft's `pget' etc. tools ##
+##                                                       ##
+## LIweitiaNux <liweitianux@gmail.com>                   ##
+## August 16, 2012                                       ##
+###########################################################
+
+###########################################################
+## ChangeLogs:
+## v1.1, 2012-08-21, LIweitiaNux
+##   fix a bug with `sed'
+## v1.2, 2012-08-21, LIweitiaNux
+##   set `ardlib' before process `merge_all'
+## v2.0, 2014/07/29, Weitian LI
+##   `merge_all' deprecated, use `fluximage' if possible
+###########################################################
+
+## about, used in `usage' {{{
+VERSION="v1.2"
+UPDATE="2012-08-21"
+## about }}}
+
+## error code {{{
+ERR_USG=1
+ERR_DIR=11
+ERR_EVT=12
+ERR_BKG=13
+ERR_REG=14
+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> energy=<e_start:e_end:e_width> basedir=<base_dir> nh=<nH> z=<redshift> temp=<avg_temperature> abund=<avg_abund> [ logfile=<log_file> ]\n"
+        printf "\nversion:\n"
+        printf "${VERSION}, ${UPDATE}\n"
+        exit ${ERR_USG}
+        ;;
+esac
+## usage, help }}}
+
+## default parameters {{{
+# default `event file' which used to match `blanksky' files
+#DFT_EVT="_NOT_EXIST_"
+DFT_EVT="`ls evt2*_clean.fits 2> /dev/null`"
+# default dir which contains `asols, asol.lis, ...' files
+# DFT_BASEDIR="_NOT_EXIST_"
+DFT_BASEDIR=".."
+# default `radial region file' to extract surface brightness
+DFT_SBP_REG="_NOT_EXIST_"
+#DFT_SBP_REG="sbprofile.reg"
+# default `energy band' applying to `images' analysis
+# format: `E_START:E_END:E_WIDTH'
+DFT_ENERGY="700:7000:100"
+
+# default `log file'
+DFT_LOGFILE="expcorr_sbp_`date '+%Y%m%d'`.log"
+
+## howto find files in `basedir'
+# default `asol.lis pattern'
+DFT_ASOLIS_PAT="acis*asol?.lis"
+# default `bad pixel filename pattern'
+DFT_BPIX_PAT="acis*repro*bpix?.fits"
+# default `msk file pattern'
+DFT_MSK_PAT="acis*msk?.fits"
+## 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}"
+printf "## PATH: ${PATH}\n"
+## check CIAO }}}
+
+## parameters {{{
+# process cmdline args using `getopt_keyval'
+getopt_keyval "$@"
+
+## check log parameters {{{
+if [ ! -z "${log}" ]; then
+    LOGFILE="${log}"
+else
+    LOGFILE=${DFT_LOGFILE}
+fi
+printf "## use logfile: \`${LOGFILE}'\n"
+[ -e "${LOGFILE}" ] && mv -fv ${LOGFILE} ${LOGFILE}_bak
+TOLOG="tee -a ${LOGFILE}"
+echo "process script: `basename $0`" >> ${LOGFILE}
+echo "process date: `date`" >> ${LOGFILE}
+## log }}}
+
+# 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}
+fi
+printf "## use reg file(s): \`${SBP_REG}'\n" | ${TOLOG}
+## check given `energy' {{{
+if [ ! -z "${energy}" ]; then
+    ENERGY="${energy}"
+else
+    ENERGY="${DFT_ENERGY}"
+fi
+# split energy variable
+ENG_N=`echo ${ENERGY} | awk -F':' '{ print NF }'`
+if [ ${ENG_N} -eq 1 ]; then
+    E_START=`echo ${DFT_ENERGY} | awk -F':' '{ print $1 }'`
+    E_END=`echo ${DFT_ENERGY} | awk -F':' '{ print $2 }'`
+    E_WIDTH=${ENERGY}
+elif [ ${ENG_N} -eq 2 ]; then
+    E_START=`echo ${ENERGY} | awk -F':' '{ print $1 }'`
+    E_END=`echo ${ENERGY} | awk -F':' '{ print $2 }'`
+    E_WIDTH=`echo ${DFT_ENERGY} | awk -F':' '{ print $3 }'`
+elif [ ${ENG_N} -eq 3 ]; then
+    E_START=`echo ${ENERGY} | awk -F':' '{ print $1 }'`
+    E_END=`echo ${ENERGY} | awk -F':' '{ print $2 }'`
+    E_WIDTH=`echo ${ENERGY} | awk -F':' '{ print $3 }'`
+else
+    printf "ERROR: invalid energy: \`${ENERGY}'\n"
+    exit ${ERR_ENG}
+fi
+ENG_RANGE="${E_START}:${E_END}"
+printf "## use energy range: \`${ENG_RANGE}' eV\n" | ${TOLOG}
+printf "## use energy width: \`${E_WIDTH}' eV\n" | ${TOLOG}
+## parse energy }}}
+# parameters (nH, z, avg_temp, avg_abund) used to calculate {{{
+# `spectral weights' for `mkinstmap'
+# check given nH
+if [ -z "${nh}" ]; then
+    read -p "> value of nH: " N_H
+else
+    N_H=${nh}
+fi
+printf "## use nH: ${N_H}\n" | ${TOLOG}
+# check given redshift
+if [ -z "${z}" ]; then
+    read -p "> value of redshift: " REDSHIFT
+else
+    REDSHIFT=${z}
+fi
+printf "## use redshift: ${REDSHIFT}\n" | ${TOLOG}
+# check given temperature
+if [ -z "${temp}" ]; then
+    read -p "> object average temperature: " TEMP
+else
+    TEMP=${temp}
+fi
+printf "## use temperature: ${TEMP}\n" | ${TOLOG}
+# check given abundance
+if [ -z "${abund}" ]; then
+    read -p "> object average abundance: " ABUND
+else
+    ABUND=${abund}
+fi
+printf "## use abundance: ${ABUND}\n" | ${TOLOG}
+# `spectral weights' parameters }}}
+# 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}
+## parameters }}}
+
+## check files in `basedir' {{{
+# check asol files
+ASOLIS=`ls -1 ${BASEDIR}/${DFT_ASOLIS_PAT} | head -n 1`
+if [ -z "${ASOLIS}" ]; then
+    printf "ERROR: cannot find \"${DFT_ASOLIS_PAT}\" in dir \`${BASEDIR}'\n"
+    exit ${ERR_ASOL}
+fi
+printf "## use asolis: \`${ASOLIS}'\n"
+# check badpixel file
+BPIX=`ls -1 ${BASEDIR}/${DFT_BPIX_PAT} | head -n 1`
+if [ -z "${BPIX}" ]; then
+    printf "ERROR: cannot find \"${DFT_BPIX_PAT}\" in dir \`${BASEDIR}'\n"
+    exit ${ERR_BPIX}
+fi
+printf "## use badpixel: \`${BPIX}'\n" | ${TOLOG}
+# check msk file
+MSK=`ls -1 ${BASEDIR}/${DFT_MSK_PAT} | head -n 1`
+if [ -z "${MSK}" ]; then
+    printf "ERROR: cannot find \"${DFT_MSK_PAT}\" in dir \`${BASEDIR}'\n"
+    exit ${ERR_MSK}
+fi
+printf "## use msk file: \`${MSK}'\n" | ${TOLOG}
+## check files }}}
+
+## determine ACIS type {{{
+# consistent with `ciao_procevt'
+punlearn dmkeypar
+DETNAM=`dmkeypar ${EVT} DETNAM echo=yes`
+if echo ${DETNAM} | grep -q 'ACIS-0123'; then
+    printf "## \`DETNAM' (${DETNAM}) has chips 0123\n"
+    printf "## ACIS-I\n"
+    ACIS_TYPE="ACIS-I"
+    CCD="0:3"
+    NEW_DETNAM="ACIS-0123"
+    ROOTNAME="c0-3_e${E_START}-${E_END}"
+elif echo ${DETNAM} | grep -q 'ACIS-[0-6]*7'; then
+    printf "## \`DETNAM' (${DETNAM}) has chip 7\n"
+    printf "## ACIS-S\n"
+    ACIS_TYPE="ACIS-S"
+    CCD="7"
+    NEW_DETNAM="ACIS-7"
+    ROOTNAME="c7_e${E_START}-${E_END}"
+else
+    printf "ERROR: unknown detector type: ${DETNAM}\n"
+    exit ${ERR_DET}
+fi
+## ACIS type }}}
+
+
+## main process {{{
+## set `ardlib' at first
+printf "set \`ardlib' first ...\n"
+punlearn ardlib
+acis_set_ardlib badpixfile="${BPIX}"
+
+## generate `spectral weights' for `instrumental map' {{{
+printf "generate \`spectral weights' for making instrumental map ...\n"
+SPEC_WGT="instmap_weights.txt"
+# convert `eV' to `keV'
+EMIN=`echo ${E_START} | awk '{ print $1/1000 }'`
+EMAX=`echo ${E_END} | awk '{ print $1/1000 }'`
+EWIDTH=`echo ${E_WIDTH} | awk '{ print $1/1000 }'`
+punlearn make_instmap_weights
+make_instmap_weights outfile="${SPEC_WGT}" \
+    model="xswabs.gal*xsapec.ap" \
+    paramvals="gal.nh=${N_H};ap.kt=${TEMP};ap.abundanc=${ABUND};ap.redshift=${REDSHIFT}" \
+    emin="${EMIN}" emax="${EMAX}" ewidth="${EWIDTH}" \
+    abund="grsa" clobber=yes
+## spectral weights }}}
+
+## generate `skyfov'
+# XXX: omit `aspec', NOT provide `asol' file
+# otherwise the size of evt_img NOT match the `expmap'
+printf "generate skyfov ...\n"
+SKYFOV="_skyfov.fits"
+[ -e "${SKYFOV}" ] && rm -fv ${SKYFOV}
+punlearn skyfov
+skyfov infile="${EVT}" outfile="${SKYFOV}" aspect="@${ASOLIS}" clobber=yes
+
+## filter by energy band & make image
+printf "filter out events in energy band: \`${ENG_RANGE}' ...\n"
+EVT_E="evt2_${ROOTNAME}.fits"
+if [ ! -r "${EVT_E}" ]; then
+    punlearn dmcopy
+    dmcopy infile="${EVT}[energy=${E_START}:${E_END}]" outfile="${EVT_E}" clobber=yes
+fi
+
+printf "make image ...\n"
+IMG_ORIG="img_${ROOTNAME}.fits"
+punlearn dmcopy
+dmcopy infile="${EVT_E}[sky=region(${SKYFOV}[ccd_id=${CCD}])][bin sky=1]" \
+    outfile="${IMG_ORIG}" clobber=yes
+
+## modify `DETNAM' of image
+printf "modify keyword \`DETNAM' of image -> \`${NEW_DETNAM}'\n"
+punlearn dmhedit
+dmhedit infile="${IMG_ORIG}" filelist=none operation=add \
+    key=DETNAM value="${NEW_DETNAM}"
+
+## get `xygrid' for image
+punlearn get_sky_limits
+get_sky_limits image="${IMG_ORIG}"
+XYGRID=`pget get_sky_limits xygrid`
+printf "## get \`xygrid': \`${XYGRID}'\n" | ${TOLOG}
+
+EXPMAP="expmap_${ROOTNAME}.fits"
+IMG_EXPCORR="img_expcorr_${ROOTNAME}.fits"
+
+if `which merge_all >/dev/null 2>&1`; then
+    # merge_all available
+    printf "merge_all ...\n"
+    ## set `ardlib' again to make sure the matched bpix file specified
+    printf "set \`ardlib' again for \`merge_all' ...\n"
+    punlearn ardlib
+    acis_set_ardlib badpixfile="${BPIX}"
+    
+    ## XXX: `merge_all' needs `asol files' in working directory
+    printf "link asol files into currect dir (\`merge_all' needed) ...\n"
+    for f in `cat ${ASOLIS}`; do
+        ln -sv ${BASEDIR}/${f} .
+    done
+    
+    printf "use \`merge_all' to generate \`exposure map' ONLY ...\n"
+    punlearn merge_all
+    merge_all evtfile="${EVT_E}" asol="@${ASOLIS}" \
+        chip="${CCD}" xygrid="${XYGRID}" \
+        energy="${SPEC_WGT}" expmap="${EXPMAP}" \
+        dtffile="" refcoord="" merged="" expcorr="" \
+        clobber=yes
+    
+    ## apply exposure correction
+    printf "use \`dmimgcalc' to apply \`exposure correction' ...\n"
+    punlearn dmimgcalc
+    dmimgcalc infile="${IMG_ORIG}" infile2="${EXPMAP}" \
+        outfile="${IMG_EXPCORR}" operation=div clobber=yes
+else
+    ## `merge_all' deprecated and not available
+    ## use 'fluximage' to generate `exposure map' and apply exposure correction
+    printf "fluximage ...\n"
+    punlearn fluximage
+    fluximage infile="${EVT_E}" outroot="${ROOTNAME}" \
+        binsize=1 bands="${SPEC_WGT}" xygrid="${XYGRID}" \
+        asol="@${ASOLIS}" badpixfile="${BPIX}" \
+        maskfile="${MSK}" clobber=yes
+    ## make symbolic links
+    # clipped counts image
+    ln -svf ${ROOTNAME}*band*thresh.img ${IMG_ORIG%.fits}_thresh.fits
+    # clipped exposure map
+    ln -svf ${ROOTNAME}*band*thresh.expmap ${EXPMAP}
+    # exposure-corrected image
+    ln -svf ${ROOTNAME}*band*flux.img ${IMG_EXPCORR}
+fi
+
+## main }}}
+
+exit 0
+