Updates scripts aganist CIAO-4.6

5 files changed, 308 insertions, 28 deletions

diff --git a/README.md b/README.md
index 9c6cbdb..77425bc 100644
--- a/README.md
+++ b/README.md
@@ -4,6 +4,7 @@ README.md
 
 Install/Compile
 ---------------
+needs `libgsl-dev'
 (1) mass_profile
 $ cd mass_profile
 $ make clean
diff --git a/scripts/ciao_bkg_spectra_v4.sh b/scripts/ciao_bkg_spectra_v4.sh
index bc40151..b4dbbbd 100755
--- a/scripts/ciao_bkg_spectra_v4.sh
+++ b/scripts/ciao_bkg_spectra_v4.sh
@@ -31,6 +31,8 @@ export LC_COLLATE=C
 #     (through cmdline which similar to CIAO,
 #     and default filename match patterns)
 #   add simple `logging' function
+# v4.1, 2014/07/29, Weitian LI
+#   fix 'pbkfile' parameters for CIAO-4.6
 ###########################################################
 
 ## about, used in `usage' {{{
@@ -307,11 +309,18 @@ for reg_i in ${REGLIST}; do
     # NO background response files
     # NO background spectrum (generate by self)
     # NO spectrum grouping (group by self using `grppha')
+    # 'pbkfile' parameter deprecated in CIAO-4.6
+    if `pget specextract pbkfile >/dev/null 2>&1`; then
+        P_PBKFILE="pbkfile='${PBK}'"
+    else
+        P_PBKFILE=""
+    fi
+    #
     printf "use \`specextract' to generate spectra and response ...\n"
     punlearn specextract
     specextract infile="${EVT}[sky=region(${REG_TMP})]" \
         outroot=${LBKG_PI%.pi} bkgfile="" asp="@${ASOLIS}" \
-        pbkfile="${PBK}" mskfile="${MSK}" badpixfile="${BPIX}" \
+        ${P_PBKFILE} mskfile="${MSK}" badpixfile="${BPIX}" \
         weight=yes correct=no bkgresp=no \
         energy="0.3:11.0:0.01" channel="1:1024:1" \
         combine=no binarfwmap=2 \
diff --git a/scripts/ciao_deproj_spectra_v8.sh b/scripts/ciao_deproj_spectra_v8.sh
index 3818a44..abebc7c 100755
--- a/scripts/ciao_deproj_spectra_v8.sh
+++ b/scripts/ciao_deproj_spectra_v8.sh
@@ -39,6 +39,10 @@ export LC_COLLATE=C
 ## v8, 2012/08/14, LIweitiaNux
 ##   use `cmdline' args instead of `cfg file'
 ##   add `logging' function
+## v8.1, 2014/07/29, Weitian LI
+##   fix variable 'ABUND=grsa'
+## v8.2, 2014/07/29, Weitian LI
+##   fix 'pbkfile' parameters for CIAO-4.6
 ###########################################################
 
 ## about, used in `usage' {{{
@@ -88,6 +92,9 @@ DFT_BPIX_PAT="acis*repro*bpix?.fits"
 DFT_PBK_PAT="acis*pbk?.fits"
 # default `msk file pattern'
 DFT_MSK_PAT="acis*msk?.fits"
+
+## abundance standard
+ABUND="grsa"
 ## default parameters }}}
 
 ## error code {{{
@@ -416,10 +423,17 @@ for i in `seq ${LINES}`; do
     # NO background response files
     # NO background spectrum (generate by self)
     # NO spectrum grouping (group by self using `grppha')
+    # 'pbkfile' parameter deprecated in CIAO-4.6
+    if `pget specextract pbkfile >/dev/null 2>&1`; then
+        P_PBKFILE="pbkfile='${PBK}'"
+    else
+        P_PBKFILE=""
+    fi
+    #
     punlearn specextract
     specextract infile="${EVT}[sky=region(${REG_CIAO})]" \
         outroot="r${i}_${ROOTNAME}" bkgfile="" asp="@${ASOLIS}" \
-        pbkfile="${PBK}" mskfile="${MSK}" badpixfile="${BPIX}" \
+        ${P_PBKFILE} mskfile="${MSK}" badpixfile="${BPIX}" \
         weight=yes correct=no bkgresp=no \
         energy="0.3:11.0:0.01" channel="1:1024:1" \
         combine=no binarfwmap=2 \
diff --git a/scripts/ciao_expcorr_only.sh b/scripts/ciao_expcorr_only.sh
index 526540d..1321382 100755
--- a/scripts/ciao_expcorr_only.sh
+++ b/scripts/ciao_expcorr_only.sh
@@ -5,7 +5,7 @@ export LC_COLLATE=C
 ###########################################################
 ## make `image' from `evt file'                          ##
 ## make `spectral weight' using `make_instmap_weights'   ##
-## use `merge_all' to generating `exposure map'          ##
+## use `fluximage' to generating `exposure map'          ##
 ## make `exposure-corrected' image                       ##
 ## and extract `surface brighness profile'               ##
 ##                                                       ##
@@ -23,6 +23,8 @@ export LC_COLLATE=C
 ##   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' {{{
@@ -81,6 +83,8 @@ DFT_LOGFILE="expcorr_sbp_`date '+%Y%m%d'`.log"
 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 {{{
@@ -233,11 +237,6 @@ if [ -z "${ASOLIS}" ]; then
     exit ${ERR_ASOL}
 fi
 printf "## use asolis: \`${ASOLIS}'\n"
-## 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
 # check badpixel file
 BPIX=`ls -1 ${BASEDIR}/${DFT_BPIX_PAT} | head -n 1`
 if [ -z "${BPIX}" ]; then
@@ -245,6 +244,13 @@ if [ -z "${BPIX}" ]; then
     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 {{{
@@ -328,29 +334,46 @@ get_sky_limits image="${IMG_ORIG}"
 XYGRID=`pget get_sky_limits xygrid`
 printf "## get \`xygrid': \`${XYGRID}'\n" | ${TOLOG}
 
-## for simplicity, use `merge_all' to generating `exposure map' {{{
-## 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}"
-
-printf "use \`merge_all' to generate \`exposure map' ONLY ...\n"
 EXPMAP="expmap_${ROOTNAME}.fits"
-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
-## `merge_all' }}}
-
-## apply exposure correction
-printf "use \`dmimgcalc' to apply \`exposure correction' ...\n"
 IMG_EXPCORR="img_expcorr_${ROOTNAME}.fits"
-punlearn dmimgcalc
-dmimgcalc infile="${IMG_ORIG}" infile2="${EXPMAP}" \
-    outfile="${IMG_EXPCORR}" operation=div clobber=yes
 
+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 verbose=2
+fi
 
 ## main }}}
 
diff --git a/scripts/ciao_expcorr_sbp_v4.sh b/scripts/ciao_expcorr_sbp_v4.sh
new file mode 100755
index 0000000..861c6dd
--- /dev/null
+++ b/scripts/ciao_expcorr_sbp_v4.sh
@@ -0,0 +1,233 @@
+#!/bin/sh
+unalias -a
+export LC_COLLATE=C
+
+#####################################################################
+## make exposure map and exposure-corrected image
+## extract surface brightness profile (revoke 'ciao_sbp_v3.1.sh')
+##
+## ChangeLogs:
+## v4, 2013/10/12, LIweitiaNux
+##   split out the 'generate regions' parts -> 'ciao_genreg_v1.sh'
+## v3, 2013/05/03, LIweitiaNux
+##   add parameter 'ds9' to check the centroid and regions
+#####################################################################
+
+SCRIPT_PATH=`readlink -f $0`
+SCRIPT_DIR=`dirname ${SCRIPT_PATH}`
+EXPCORR_SCRIPT="ciao_expcorr_only.sh"
+EXTRACT_SBP_SCRIPT="ciao_sbp_v3.1.sh"
+
+## about, used in `usage' {{{
+VERSION="v4"
+UPDATE="2013-10-12"
+## about }}}
+
+## err 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 {{{
+case "$1" in
+    -[hH]*|--[hH]*)
+        printf "usage:\n"
+        printf "    `basename $0` evt=<evt2_clean> sbp_reg=<sbprofile.reg> nh=<nh> z=<redshift> temp=<avg_temp> abund=<avg_abund> cellreg=<celld_reg> expcorr=<yes|no>\n"
+        printf "\nversion:\n"
+        printf "${VERSION}, ${UPDATE}\n"
+        exit ${ERR_USG}
+        ;;
+esac
+## usage }}}
+
+## default parameters {{{ 
+## clean evt2 file
+DFT_EVT=`ls evt2*_clean.fits 2> /dev/null`
+## the repro dir
+DFT_BASEDIR=".."
+
+## xcentroid region
+DFT_SBP_REG="sbprofile.reg"
+
+## log file
+DFT_LOGFILE="expcorr_sbp_`date '+%Y%m%d'`.log"
+
+## cell region
+DFT_CELL_REG=`ls celld*.reg 2> /dev/null`
+
+## background spectra
+DFT_BKGD=`ls bkgcorr_bl*.pi 2> /dev/null`
+## default parameters }}}
+
+## functions {{{
+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 and set path to solve conflit with heasoft {{{
+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&heasoft}}}
+
+## parameters {{{
+getopt_keyval "$@"
+
+##log file
+LOGFILE="not_exist"
+
+##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
+
+##check ori sbp region
+if [ -r "${sbp_reg}" ]; then
+    SBP_REG="${sbp_reg}"
+elif [ -r "${DFT_SBP_REG}" ] ;then
+    SBP_REG="${DFT_SBP_REG}" 
+else
+    read -p "sbp region file: " SBP_REG
+    if [ ! -r "${SBP_REG}" ]; then
+        printf "ERROR: cannot access given \`${SBP}' sbp region file\n"
+        exit ${ERR_REG}
+    fi
+fi
+
+## nh z temp abund
+if [ -z "${nh}" ]; then
+    read -p "> value of nH: " N_H
+else
+    N_H=${nh}
+fi
+if [ -z "${z}" ]; then
+    read -p "> value of redshift: " REDSHIFT
+else
+    REDSHIFT=${z}
+fi
+if [ -z "${temp}" ]; then
+    read -p "> object average temperature: " TEMP
+else
+    TEMP=${temp}
+fi
+if [ -z "${abund}" ]; then
+    read -p "> object average abundance: " ABUND
+else
+    ABUND=${abund}
+fi
+# check give basedir
+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/\/*$//'`
+
+# point source region
+if [ -r "${cellreg}" ]; then
+    CELL_REG=${cellreg}
+elif [ -r "${DFT_CELL_REG}" ]; then
+    CELL_REG=${DFT_CELL_REG}
+else
+    read -p ">point source region file: " CELL_REG
+    if [ ! -d ${CELL_REG} ] ; then
+        printf " ERROR no point source region\n"
+        exit ${ERR_REG}
+    fi
+fi
+
+## expcorr: flag to determine whether to process expcorr
+if [ ! -z "${expcorr}" ]; then
+    case "${expcorr}" in
+        [nN][oO]|[fF]*)
+            F_EXPCORR="NO"
+            ;;
+        *)
+            F_EXPCORR="YES"
+            ;;
+    esac
+else
+    F_EXPCORR="YES"
+fi
+## parameters }}}
+
+if [ "${F_EXPCORR}" = "NO" ]; then
+    printf "################################\n"
+    printf "### SKIP EXPOSURE CORRECTION ###\n"
+    printf "################################\n\n"
+else
+    printf "======== EXPOSURE CORRECTION =======\n"
+    CMD="${SCRIPT_DIR}/${EXPCORR_SCRIPT} evt=${EVT} basedir=${BASEDIR} nh=${N_H} z=${REDSHIFT} temp=${TEMP} abund=${ABUND}"
+    printf "CMD: ${CMD}\n"
+    ${SCRIPT_DIR}/${EXPCORR_SCRIPT} evt=${EVT} basedir=${BASEDIR} nh=${N_H} z=${REDSHIFT} temp=${TEMP} abund=${ABUND}
+    printf "======== EXPOSURE CORRECTION FINISHED =======\n\n"
+fi
+
+EXPMAP=`ls expmap*e700-7000*fits 2> /dev/null`
+EVT_E=`ls evt*e700-7000*fits 2> /dev/null`
+
+printf "======== EXTRACT SBP =======\n"
+CMD="${SCRIPT_DIR}/${EXTRACT_SBP_SCRIPT} evt_e=${EVT_E} reg=${SBP_REG} expmap=${EXPMAP} cellreg=${CELL_REG}"
+printf "CMD: ${CMD}\n"
+${SCRIPT_DIR}/${EXTRACT_SBP_SCRIPT} evt_e=${EVT_E} reg=${SBP_REG} expmap=${EXPMAP} cellreg=${CELL_REG}
+printf "======== EXTRACT SBP FINISHED =======\n\n"
+
+#generate a cfg file for specextract
+[ -e "spc_fit.cfg" ] && mv -fv spc_fit.cfg spc_fit.cfg_bak
+cat > spc_fit.cfg << _EOF_
+nh     ${N_H}
+z      ${REDSHIFT}
+basedir ../..
+bkgd ../${BKGD}
+_EOF_
+
+#link for sbp2.dat radius.dat sbp3.dat
+ln -svf radius_sbp.txt radius.dat
+ln -svf flux_sbp.txt sbp2.dat
+ln -svf sbprofile.txt sbp3.dat
+
+exit 0
+