From 92c9b0121001d5cab4127f4ae5ae95f43958513a Mon Sep 17 00:00:00 2001 From: Weitian LI Date: Thu, 30 Oct 2014 18:55:16 +0800 Subject: fix to 'ciao_expcorr.sh' & removed versions in filenames --- scripts/ciao_bkg_spectra.sh | 425 +++++++++++++++++++++++ scripts/ciao_bkg_spectra_v4.sh | 425 ----------------------- scripts/ciao_blanksky.sh | 242 ++++++++++++++ scripts/ciao_blanksky_v4.sh | 242 -------------- scripts/ciao_deproj_spectra.sh | 684 ++++++++++++++++++++++++++++++++++++++ scripts/ciao_deproj_spectra_v8.sh | 684 -------------------------------------- scripts/ciao_expcorr.sh | 388 +++++++++++++++++++++ scripts/ciao_expcorr_only.sh | 388 --------------------- scripts/ciao_expcorr_sbp.sh | 236 +++++++++++++ scripts/ciao_expcorr_sbp_v4.sh | 233 ------------- scripts/ciao_genregs.sh | 277 +++++++++++++++ scripts/ciao_genregs_v1.sh | 277 --------------- scripts/ciao_procevt.sh | 185 +++++++++++ scripts/ciao_procevt_v2.sh | 182 ---------- scripts/ciao_r500avgt.sh | 539 ++++++++++++++++++++++++++++++ scripts/ciao_r500avgt_v3.sh | 539 ------------------------------ scripts/ciao_sbp.sh | 529 +++++++++++++++++++++++++++++ scripts/ciao_sbp_v3.1.sh | 529 ----------------------------- 18 files changed, 3505 insertions(+), 3499 deletions(-) create mode 100755 scripts/ciao_bkg_spectra.sh delete mode 100755 scripts/ciao_bkg_spectra_v4.sh create mode 100755 scripts/ciao_blanksky.sh delete mode 100755 scripts/ciao_blanksky_v4.sh create mode 100755 scripts/ciao_deproj_spectra.sh delete mode 100755 scripts/ciao_deproj_spectra_v8.sh create mode 100755 scripts/ciao_expcorr.sh delete mode 100755 scripts/ciao_expcorr_only.sh create mode 100755 scripts/ciao_expcorr_sbp.sh delete mode 100755 scripts/ciao_expcorr_sbp_v4.sh create mode 100755 scripts/ciao_genregs.sh delete mode 100755 scripts/ciao_genregs_v1.sh create mode 100755 scripts/ciao_procevt.sh delete mode 100755 scripts/ciao_procevt_v2.sh create mode 100755 scripts/ciao_r500avgt.sh delete mode 100755 scripts/ciao_r500avgt_v3.sh create mode 100755 scripts/ciao_sbp.sh delete mode 100755 scripts/ciao_sbp_v3.1.sh (limited to 'scripts') diff --git a/scripts/ciao_bkg_spectra.sh b/scripts/ciao_bkg_spectra.sh new file mode 100755 index 0000000..b4dbbbd --- /dev/null +++ b/scripts/ciao_bkg_spectra.sh @@ -0,0 +1,425 @@ +#!/bin/sh - +# +trap date INT +unalias -a +export GREP_OPTIONS="" +export LC_COLLATE=C +########################################################### +## extract background spectra from src and blanksky ## +## renormalization the blank spectrum ## +## ## +## Ref: Chandra spectrum analysis ## +## http://cxc.harvard.edu/ciao/threads/extended/ ## +## Ref: specextract ## +## http://cxc.harvard.edu/ciao/ahelp/specextract.html ## +## Ref: CIAO v4.4 region bugs ## +## http://cxc.harvard.edu/ciao/bugs/regions.html#bug-12187 +## ## +## LIweitiaNux, July 24, 2012 ## +########################################################### + +########################################################### +# ChangeLogs: +# v3, 2012/08/09 +# fix `scientific notation' for `bc' +# change `spec group' method to `min 15' +# v4, 2012/08/13 +# add `clobber=yes' +# improve error code +# improve cmdline arguements +# provide a flexible way to pass parameters +# (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' {{{ +VERSION="v4" +UPDATE="2012-08-14" +## about }}} + +## usage, help {{{ +case "$1" in + -[hH]*|--[hH]*) + printf "usage:\n" + printf " `basename $0` evt= reg= blank= basedir= nh= z= [ grpcmd= log= ]\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`" +# default `blanksky file' +#DFT_BLANK="_NOT_EXIST_" +DFT_BLANK="`ls blanksky*.fits`" +# default dir which contains `asols, asol.lis, ...' files +#DFT_BASEDIR="_NOT_EXIST_" +DFT_BASEDIR=".." +# 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="bkg_spectra_`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 `pbk file pattern' +DFT_PBK_PAT="acis*pbk?.fits" +# default `msk file pattern' +DFT_MSK_PAT="acis*msk?.fits" +## default parameters }}} + +## 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 +## error code }}} + +## 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 +} + +## background renormalization (BACKSCAL) {{{ +# renorm background according to particle background +# energy range: 9.5-12.0 keV (channel: 651-822) +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 }'` + punlearn dmkeypar + EXPTIME=`dmkeypar $1 EXPOSURE echo=yes` + BACK=`dmkeypar $1 BACKSCAL echo=yes` + # fix `scientific notation' bug for `bc' + EXPTIME_B=`echo ${EXPTIME} | sed 's/[eE]/\*10\^/' | sed 's/+//'` + BACK_B=`echo "( ${BACK} )" | sed 's/[eE]/\*10\^/' | sed 's/+//'` + PB_FLUX=`echo "scale = 16; ${COUNTS} / ${EXPTIME_B} / ${BACK_B}" | bc -l` + echo ${PB_FLUX} +} + +bkg_renorm() { + # $1: src spectrum, $2: back spectrum + PBFLUX_SRC=`pb_flux $1` + PBFLUX_BKG=`pb_flux $2` + BACK_OLD=`dmkeypar $2 BACKSCAL echo=yes` + BACK_OLD_B=`echo "( ${BACK_OLD} )" | sed 's/[eE]/\*10\^/' | sed 's/+//'` + BACK_NEW=`echo "scale = 16; ${BACK_OLD_B} * ${PBFLUX_BKG} / ${PBFLUX_SRC}" | bc -l` + printf "\`$2': BACKSCAL:\n" + printf " ${BACK_OLD} --> ${BACK_NEW}\n" + punlearn dmhedit + dmhedit infile=$2 filelist=none operation=add \ + key=BACKSCAL value=${BACK_NEW} comment="old value: ${BACK_OLD}" +} +## bkg renorm }}} +## functions end }}} + +## 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 [ -z "${reg}" ]; then + read -p "> selected local bkg region file: " REGLIST +else + REGLIST="${reg}" +fi +REGLIST=`echo ${REGLIST} | tr ',' ' '` # use *space* to separate +printf "## use reg file(s): \`${REGLIST}'\n" | ${TOLOG} +# check given blanksky +if [ -r "${blank}" ]; then + BLANK=${blank} +elif [ -r "${DFT_BLANK}" ]; then + BLANK=${DFT_BLANK} +else + read -p "> matched blanksky evtfile: " BLANK + if [ ! -r "${BLANK}" ]; then + printf "ERROR: cannot acces given \`${BLANK}' blanksky file\n" + exit ${ERR_BKG} + fi +fi +# 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 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 given `grpcmd' +if [ ! -z "${grpcmd}" ]; then + GRP_CMD="${grpcmd}" +else + GRP_CMD="${DFT_GRP_CMD}" +fi +printf "## use grppha cmd: \`${GRP_CMD}'\n" | ${TOLOG} +## parameters }}} + +## check needed files {{{ +# check reg file(s) +printf "check accessibility of reg file(s) ...\n" +for reg_f in ${REGLIST}; do + if [ ! -r "${reg_f}" ]; then + printf "ERROR: file \`${reg_f}' NOT accessiable\n" + exit ${ERR_REG} + fi +done +# check the validity of *pie* regions +printf "check pie reg validity ...\n" +INVALID=`cat ${REGLIST} | grep -i 'pie' | awk -F, '{ print $6 }' | tr -d ')' | awk '$1 > 360'` +if [ "x${INVALID}" != "x" ]; then + printf "WARNING: some pie region's END_ANGLE > 360\n" | ${TOLOG} + printf " CIAO v4.4 tools may run into trouble\n" +fi + +# check files in `basedir' +printf "check needed files in basedir \`${BASEDIR}' ...\n" +# check asolis 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" | ${TOLOG} +# 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 pbk file +PBK=`ls -1 ${BASEDIR}/${DFT_PBK_PAT} | head -n 1` +if [ -z "${PBK}" ]; then + printf "ERROR: cannot find \"${DFT_PBK_PAT}\" in dir \`${BASEDIR}'\n" + exit ${ERR_PBK} +fi +printf "## use pbk: \`${PBK}'\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: \`${MSK}'\n" | ${TOLOG} +## check files }}} + +## main part {{{ +## use 'for' loop to process every region file +for reg_i in ${REGLIST}; do + printf "\n==============================\n" + printf "PROCESS REGION fle \`${reg_i}' ...\n" + + REG_TMP="_tmp.reg" + [ -f "${REG_TMP}" ] && rm -fv ${REG_TMP} # remove tmp files + cp -fv ${reg_i} ${REG_TMP} + # check the validity of *pie* regions {{{ + INVALID=`grep -i 'pie' ${REG_TMP} | awk -F, '{ print $6 }' | tr -d ')' | awk '$1 > 360'` + if [ "x${INVALID}" != "x" ]; then + printf "WARNING: fix for *pie* region in file \`${reg_i}'\n" + cat ${REG_TMP} + A_OLD=`echo ${INVALID} | sed 's/\./\\\./'` + A_NEW=`echo ${INVALID}-360 | bc -l | sed 's/\./\\\./'` + sed -i'' "s/${A_OLD}\ *)/${A_NEW})/" ${REG_TMP} + printf " --> " + cat ${REG_TMP} + fi + ## check pie region }}} + + LBKG_PI="${reg_i%.reg}.pi" + ## use `specextract' to extract local bkg spectrum {{{ + # NOTE: set `binarfwmap=2' to save the time for generating `ARF' + # I have tested that this bin factor has little impact on the results. + # 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}" \ + ${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 \ + grouptype=NONE binspec=NONE \ + verbose=2 clobber=yes + # `specextract' }}} + + ## generate the blanksky bkg spectrum {{{ + printf "generate the blanksky bkg spectrum ...\n" + BBKG_PI="blanksky_${LBKG_PI}" + punlearn dmextract + dmextract infile="${BLANK}[sky=region(${REG_TMP})][bin pi]" \ + outfile=${BBKG_PI} wmap="[bin det=8]" clobber=yes + ## blanksky bkg spectrum }}} + + ## bkg renormalization {{{ + printf "Renormalize background ...\n" + bkg_renorm ${LBKG_PI} ${BBKG_PI} + ## bkg renorm }}} + + ## group spectrum {{{ + printf "group spectrum using \`grppha'\n" + LBKG_GRP_PI="${LBKG_PI%.pi}_grp.pi" + grppha infile="${LBKG_PI}" outfile="${LBKG_GRP_PI}" \ + comm="${GRP_CMD} & exit" clobber=yes > /dev/null + ## group spectra }}} + + ## generate a script for XSPEC {{{ + XSPEC_XCM="xspec_${LBKG_PI%.pi}_model.xcm" + if [ -e ${XSPEC_XCM} ]; then + mv -fv ${XSPEC_XCM} ${XSPEC_XCM}_bak + fi + cat >> ${XSPEC_XCM} << _EOF_ +## xspec script +## analysis chandra acis background components +## xspec model: apec+apec+wabs*(pow+apec) +## +## generated by script \``basename $0`' +## `date` +## NOTES: needs XSPEC v12.x + +# settings +statistic chi +#weight churazov +abund grsa +query yes + +# data +data ${LBKG_GRP_PI} +response ${LBKG_PI%.pi}.wrmf +arf ${LBKG_PI%.pi}.warf +backgrnd ${BBKG_PI} + +# fitting range +ignore bad +ignore 0.0-0.4,8.0-** + +# plot related +setplot energy + +method leven 10 0.01 +xsect bcmc +cosmo 70 0 0.73 +xset delta 0.01 +systematic 0 + +# model +model apec + apec + wabs(powerlaw + apec) + 0.08 -0.01 0.008 0.008 64 64 + 1 -0.001 0 0 5 5 + 0 -0.01 -0.999 -0.999 10 10 + 0.0 0.01 -1 0 0 1 + 0.2 -0.01 0.008 0.008 64 64 + 1 -0.001 0 0 5 5 + 0 -0.01 -0.999 -0.999 10 10 + 0.0 0.01 -1 0 0 1 + ${N_H} -0.001 0 0 100000 1e+06 + 1.4 -0.01 -3 -2 9 10 + 0.0 0.01 -1 0 0 1 + 1.0 0.01 0.008 0.008 64 64 + 0.4 0.001 0 0 5 5 + ${REDSHIFT} -0.01 -0.999 -0.999 10 10 + 0.0 0.01 0 0 1e+24 1e+24 + +freeze 1 2 3 +freeze 5 6 7 +freeze 9 10 14 +thaw 12 13 +_EOF_ + ## XSPEC script }}} + +done # end 'for', `specextract' +## main part }}} + +# clean +printf "clean ...\n" +rm -f ${REG_TMP} + +printf "DONE\n" +########################################################### + +# vim: set ts=8 sw=4 tw=0 fenc=utf-8 ft=sh: # diff --git a/scripts/ciao_bkg_spectra_v4.sh b/scripts/ciao_bkg_spectra_v4.sh deleted file mode 100755 index b4dbbbd..0000000 --- a/scripts/ciao_bkg_spectra_v4.sh +++ /dev/null @@ -1,425 +0,0 @@ -#!/bin/sh - -# -trap date INT -unalias -a -export GREP_OPTIONS="" -export LC_COLLATE=C -########################################################### -## extract background spectra from src and blanksky ## -## renormalization the blank spectrum ## -## ## -## Ref: Chandra spectrum analysis ## -## http://cxc.harvard.edu/ciao/threads/extended/ ## -## Ref: specextract ## -## http://cxc.harvard.edu/ciao/ahelp/specextract.html ## -## Ref: CIAO v4.4 region bugs ## -## http://cxc.harvard.edu/ciao/bugs/regions.html#bug-12187 -## ## -## LIweitiaNux, July 24, 2012 ## -########################################################### - -########################################################### -# ChangeLogs: -# v3, 2012/08/09 -# fix `scientific notation' for `bc' -# change `spec group' method to `min 15' -# v4, 2012/08/13 -# add `clobber=yes' -# improve error code -# improve cmdline arguements -# provide a flexible way to pass parameters -# (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' {{{ -VERSION="v4" -UPDATE="2012-08-14" -## about }}} - -## usage, help {{{ -case "$1" in - -[hH]*|--[hH]*) - printf "usage:\n" - printf " `basename $0` evt= reg= blank= basedir= nh= z= [ grpcmd= log= ]\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`" -# default `blanksky file' -#DFT_BLANK="_NOT_EXIST_" -DFT_BLANK="`ls blanksky*.fits`" -# default dir which contains `asols, asol.lis, ...' files -#DFT_BASEDIR="_NOT_EXIST_" -DFT_BASEDIR=".." -# 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="bkg_spectra_`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 `pbk file pattern' -DFT_PBK_PAT="acis*pbk?.fits" -# default `msk file pattern' -DFT_MSK_PAT="acis*msk?.fits" -## default parameters }}} - -## 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 -## error code }}} - -## 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 -} - -## background renormalization (BACKSCAL) {{{ -# renorm background according to particle background -# energy range: 9.5-12.0 keV (channel: 651-822) -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 }'` - punlearn dmkeypar - EXPTIME=`dmkeypar $1 EXPOSURE echo=yes` - BACK=`dmkeypar $1 BACKSCAL echo=yes` - # fix `scientific notation' bug for `bc' - EXPTIME_B=`echo ${EXPTIME} | sed 's/[eE]/\*10\^/' | sed 's/+//'` - BACK_B=`echo "( ${BACK} )" | sed 's/[eE]/\*10\^/' | sed 's/+//'` - PB_FLUX=`echo "scale = 16; ${COUNTS} / ${EXPTIME_B} / ${BACK_B}" | bc -l` - echo ${PB_FLUX} -} - -bkg_renorm() { - # $1: src spectrum, $2: back spectrum - PBFLUX_SRC=`pb_flux $1` - PBFLUX_BKG=`pb_flux $2` - BACK_OLD=`dmkeypar $2 BACKSCAL echo=yes` - BACK_OLD_B=`echo "( ${BACK_OLD} )" | sed 's/[eE]/\*10\^/' | sed 's/+//'` - BACK_NEW=`echo "scale = 16; ${BACK_OLD_B} * ${PBFLUX_BKG} / ${PBFLUX_SRC}" | bc -l` - printf "\`$2': BACKSCAL:\n" - printf " ${BACK_OLD} --> ${BACK_NEW}\n" - punlearn dmhedit - dmhedit infile=$2 filelist=none operation=add \ - key=BACKSCAL value=${BACK_NEW} comment="old value: ${BACK_OLD}" -} -## bkg renorm }}} -## functions end }}} - -## 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 [ -z "${reg}" ]; then - read -p "> selected local bkg region file: " REGLIST -else - REGLIST="${reg}" -fi -REGLIST=`echo ${REGLIST} | tr ',' ' '` # use *space* to separate -printf "## use reg file(s): \`${REGLIST}'\n" | ${TOLOG} -# check given blanksky -if [ -r "${blank}" ]; then - BLANK=${blank} -elif [ -r "${DFT_BLANK}" ]; then - BLANK=${DFT_BLANK} -else - read -p "> matched blanksky evtfile: " BLANK - if [ ! -r "${BLANK}" ]; then - printf "ERROR: cannot acces given \`${BLANK}' blanksky file\n" - exit ${ERR_BKG} - fi -fi -# 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 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 given `grpcmd' -if [ ! -z "${grpcmd}" ]; then - GRP_CMD="${grpcmd}" -else - GRP_CMD="${DFT_GRP_CMD}" -fi -printf "## use grppha cmd: \`${GRP_CMD}'\n" | ${TOLOG} -## parameters }}} - -## check needed files {{{ -# check reg file(s) -printf "check accessibility of reg file(s) ...\n" -for reg_f in ${REGLIST}; do - if [ ! -r "${reg_f}" ]; then - printf "ERROR: file \`${reg_f}' NOT accessiable\n" - exit ${ERR_REG} - fi -done -# check the validity of *pie* regions -printf "check pie reg validity ...\n" -INVALID=`cat ${REGLIST} | grep -i 'pie' | awk -F, '{ print $6 }' | tr -d ')' | awk '$1 > 360'` -if [ "x${INVALID}" != "x" ]; then - printf "WARNING: some pie region's END_ANGLE > 360\n" | ${TOLOG} - printf " CIAO v4.4 tools may run into trouble\n" -fi - -# check files in `basedir' -printf "check needed files in basedir \`${BASEDIR}' ...\n" -# check asolis 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" | ${TOLOG} -# 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 pbk file -PBK=`ls -1 ${BASEDIR}/${DFT_PBK_PAT} | head -n 1` -if [ -z "${PBK}" ]; then - printf "ERROR: cannot find \"${DFT_PBK_PAT}\" in dir \`${BASEDIR}'\n" - exit ${ERR_PBK} -fi -printf "## use pbk: \`${PBK}'\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: \`${MSK}'\n" | ${TOLOG} -## check files }}} - -## main part {{{ -## use 'for' loop to process every region file -for reg_i in ${REGLIST}; do - printf "\n==============================\n" - printf "PROCESS REGION fle \`${reg_i}' ...\n" - - REG_TMP="_tmp.reg" - [ -f "${REG_TMP}" ] && rm -fv ${REG_TMP} # remove tmp files - cp -fv ${reg_i} ${REG_TMP} - # check the validity of *pie* regions {{{ - INVALID=`grep -i 'pie' ${REG_TMP} | awk -F, '{ print $6 }' | tr -d ')' | awk '$1 > 360'` - if [ "x${INVALID}" != "x" ]; then - printf "WARNING: fix for *pie* region in file \`${reg_i}'\n" - cat ${REG_TMP} - A_OLD=`echo ${INVALID} | sed 's/\./\\\./'` - A_NEW=`echo ${INVALID}-360 | bc -l | sed 's/\./\\\./'` - sed -i'' "s/${A_OLD}\ *)/${A_NEW})/" ${REG_TMP} - printf " --> " - cat ${REG_TMP} - fi - ## check pie region }}} - - LBKG_PI="${reg_i%.reg}.pi" - ## use `specextract' to extract local bkg spectrum {{{ - # NOTE: set `binarfwmap=2' to save the time for generating `ARF' - # I have tested that this bin factor has little impact on the results. - # 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}" \ - ${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 \ - grouptype=NONE binspec=NONE \ - verbose=2 clobber=yes - # `specextract' }}} - - ## generate the blanksky bkg spectrum {{{ - printf "generate the blanksky bkg spectrum ...\n" - BBKG_PI="blanksky_${LBKG_PI}" - punlearn dmextract - dmextract infile="${BLANK}[sky=region(${REG_TMP})][bin pi]" \ - outfile=${BBKG_PI} wmap="[bin det=8]" clobber=yes - ## blanksky bkg spectrum }}} - - ## bkg renormalization {{{ - printf "Renormalize background ...\n" - bkg_renorm ${LBKG_PI} ${BBKG_PI} - ## bkg renorm }}} - - ## group spectrum {{{ - printf "group spectrum using \`grppha'\n" - LBKG_GRP_PI="${LBKG_PI%.pi}_grp.pi" - grppha infile="${LBKG_PI}" outfile="${LBKG_GRP_PI}" \ - comm="${GRP_CMD} & exit" clobber=yes > /dev/null - ## group spectra }}} - - ## generate a script for XSPEC {{{ - XSPEC_XCM="xspec_${LBKG_PI%.pi}_model.xcm" - if [ -e ${XSPEC_XCM} ]; then - mv -fv ${XSPEC_XCM} ${XSPEC_XCM}_bak - fi - cat >> ${XSPEC_XCM} << _EOF_ -## xspec script -## analysis chandra acis background components -## xspec model: apec+apec+wabs*(pow+apec) -## -## generated by script \``basename $0`' -## `date` -## NOTES: needs XSPEC v12.x - -# settings -statistic chi -#weight churazov -abund grsa -query yes - -# data -data ${LBKG_GRP_PI} -response ${LBKG_PI%.pi}.wrmf -arf ${LBKG_PI%.pi}.warf -backgrnd ${BBKG_PI} - -# fitting range -ignore bad -ignore 0.0-0.4,8.0-** - -# plot related -setplot energy - -method leven 10 0.01 -xsect bcmc -cosmo 70 0 0.73 -xset delta 0.01 -systematic 0 - -# model -model apec + apec + wabs(powerlaw + apec) - 0.08 -0.01 0.008 0.008 64 64 - 1 -0.001 0 0 5 5 - 0 -0.01 -0.999 -0.999 10 10 - 0.0 0.01 -1 0 0 1 - 0.2 -0.01 0.008 0.008 64 64 - 1 -0.001 0 0 5 5 - 0 -0.01 -0.999 -0.999 10 10 - 0.0 0.01 -1 0 0 1 - ${N_H} -0.001 0 0 100000 1e+06 - 1.4 -0.01 -3 -2 9 10 - 0.0 0.01 -1 0 0 1 - 1.0 0.01 0.008 0.008 64 64 - 0.4 0.001 0 0 5 5 - ${REDSHIFT} -0.01 -0.999 -0.999 10 10 - 0.0 0.01 0 0 1e+24 1e+24 - -freeze 1 2 3 -freeze 5 6 7 -freeze 9 10 14 -thaw 12 13 -_EOF_ - ## XSPEC script }}} - -done # end 'for', `specextract' -## main part }}} - -# clean -printf "clean ...\n" -rm -f ${REG_TMP} - -printf "DONE\n" -########################################################### - -# vim: set ts=8 sw=4 tw=0 fenc=utf-8 ft=sh: # diff --git a/scripts/ciao_blanksky.sh b/scripts/ciao_blanksky.sh new file mode 100755 index 0000000..da1d6f6 --- /dev/null +++ b/scripts/ciao_blanksky.sh @@ -0,0 +1,242 @@ +#!/bin/sh +# +unalias -a +export LC_COLLATE=C +########################################################### +## process `blanksky' bkg files for spectra analysis ## +## http://cxc.harvard.edu/ciao/threads/acisbackground/ ## +## ## +## for extracting the spectrum of the background ## +## for determining the background components ## +## ## +## inputs: `evt2_clean', `asol files' ## +## output: `blanksky_c7.fits' for ACIS-S, ## +## `blanksky-c0-3.fits' for ACIS-I ## +## ## +## LIweitiaNux, January 11, 2011 ## +########################################################### + +########################################################### +## ChangeLogs: +## v2: 2012/08/01 +## add ACIS-I support (chips 0-3) +## v3: 2012/08/06, Junhua Gu +## pass parameters by cmd line param +## v4: 2012/08/13, LIweitiaNux +## add `clobber=yes' parameters to CIAO tools +## improve `commandline arguements' +## add `default parameters' +########################################################### + +## about, used in `usage' {{{ +VERSION="v4" +UPDATE="2012-08-14" +## about }}} + +## usage, help {{{ +case "$1" in + -[hH]*|--[hH]*) + printf "usage:\n" + printf " `basename $0` evt= basedir= [ outfile= ]\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`" +# default dir which contains `asols, asol.lis, ...' files +# DFT_BASEDIR="_NOT_EXIST_" +DFT_BASEDIR=".." + +## howto find files in `basedir' +# default `asol.lis pattern' +DFT_ASOLIS_PAT="acis*asol?.lis" +## default parameters }}} + +## 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 +## error code }}} + +## 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 +} + +# reprocess blanksky evt with matched gainfile +blank_regain() { + mv $1 ${1%.fits}_ungain.fits + punlearn acis_process_events + acis_process_events infile="${1%.fits}_ungain.fits" \ + outfile="$1" \ + acaofffile=NONE stop="none" doevtgrade=no \ + apply_cti=yes apply_tgain=no \ + calculate_pi=yes pix_adj=NONE \ + gainfile="$2" \ + eventdef="{s:ccd_id,s:node_id,i:expno,s:chip,s:tdet,f:det,f:sky,s:phas,l:pha,l:pha_ro,f:energy,l:pi,s:fltgrade,s:grade,x:status}" \ + clobber=yes + rm -fv ${1%.fits}_ungain.fits +} +## functions end }}} + +## 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 "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" +# 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" +## 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 files }}} +#exit 0 # test script + +#### main start {{{ +printf "look up coresponding background file ...\n" +BKG_LKP="`acis_bkgrnd_lookup ${EVT}`" +AS_NUM=`echo ${BKG_LKP} | tr ' ' '\n' | \grep 'acis7sD' | wc -l` +AI_NUM=`echo ${BKG_LKP} | tr ' ' '\n' | \grep 'acis[0123]iD' | wc -l` +## determine detector type: ACIS-S / ACIS-I {{{ +if [ ${AS_NUM} -eq 1 ]; then + printf "## ACIS-S, chip: 7\n" + BKG_ROOT="blanksky_c7" + cp -v ${BKG_LKP} ${BKG_ROOT}_orig.fits +elif [ ${AI_NUM} -eq 4 ]; then + printf "## ACIS-I, chip: 0-3\n" + BKG_ROOT="blanksky_c0-3" + AI_FILES="" + for bf in ${BKG_LKP}; do + cp -v ${bf} . + AI_FILES="${AI_FILES},`basename ${bf}`" + done + AI_FILES=${AI_FILES#,} # remove the first ',' + printf "## ACIS-I blanksky files to merge:\n" + printf "## \`${AI_FILES}'\n" + printf "\`dmmerge' to merge the above blanksky files ...\n" + # merge 4 chips blanksky evt files + punlearn dmmerge + dmmerge "${AI_FILES}" ${BKG_ROOT}_orig.fits clobber=yes + rm -fv `echo ${AI_FILES} | tr ',' ' '` # remove original files +else + printf "## ERROR: UNKNOW blanksky files:\n" + printf "## ${BKG_ORIG}\n" + exit ${ERR_BKG} +fi +## determine ACIS type }}} + +## check 'DATMODE' {{{ +## filter blanksky files (status=0) for `VFAINT' observations +DATA_MODE="`dmkeypar ${EVT} DATAMODE echo=yes`" +printf "## DATAMODE: ${DATA_MODE}\n" +if [ "${DATA_MODE}" = "VFAINT" ]; then + mv -fv ${BKG_ROOT}_orig.fits ${BKG_ROOT}_tmp.fits + printf "apply \`status=0' to filter blanksky file ...\n" + punlearn dmcopy + dmcopy "${BKG_ROOT}_tmp.fits[status=0]" ${BKG_ROOT}_orig.fits clobber=yes + rm -fv ${BKG_ROOT}_tmp.fits +fi +## DATAMODE, status=0 }}} + +## check `GAINFILE' of blanksky and evt2 file {{{ +## if NOT match, then reprocess blanksky +GAINFILE_EVT="`dmkeypar ${EVT} GAINFILE echo=yes`" +GAINFILE_BG="`dmkeypar "${BKG_ROOT}_orig.fits" GAINFILE echo=yes`" +if ! [ "${GAINFILE_EVT}" = "${GAINFILE_BG}" ]; then + printf "WARNING: GAINFILE NOT match.\n" + printf "event: ${GAINFILE_EVT}\n" + printf "blank: ${GAINFILE_BG}\n" + printf "reprocess blanksky with evt gainfile ...\n" + # reprocess blanksky using matched evt GAINFILE + GAINFILE="$CALDB/data/chandra/acis/det_gain/`basename ${GAINFILE_EVT}`" + printf "GAINFILE: ${GAINFILE}\n" + blank_regain "${BKG_ROOT}_orig.fits" ${GAINFILE} +fi +## check & match GAINFILE }}} + +printf "add the PNT header keywords ... " +EVT_HEADER="_evt_header.par" +EVT_PNT="_evt_pnt.par" +punlearn dmmakepar +dmmakepar ${EVT} ${EVT_HEADER} clobber=yes +grep -i '_pnt' ${EVT_HEADER} > ${EVT_PNT} +punlearn dmreadpar +dmreadpar ${EVT_PNT} "${BKG_ROOT}_orig.fits[EVENTS]" clobber=yes +printf "DONE\n" + +printf "reproject the background ...\n" +punlearn reproject_events +reproject_events infile=${BKG_ROOT}_orig.fits \ + outfile=${BKG_ROOT}.fits match=${EVT} \ + aspect="@${ASOLIS}" random=0 clobber=yes + +# rename output file if specified +if ! [ -z "${outfile}" ]; then + mv -fv ${BKG_ROOT}.fits ${outfile} +fi +## main end }}} + +# clean +printf "\nclean ...\n" +rm -fv ${BKG_ROOT}_orig.fits ${EVT_PNT} + +printf "\nFINISHED\n" + +# vim: set ts=8 sw=4 tw=0 fenc=utf-8 ft=sh # diff --git a/scripts/ciao_blanksky_v4.sh b/scripts/ciao_blanksky_v4.sh deleted file mode 100755 index da1d6f6..0000000 --- a/scripts/ciao_blanksky_v4.sh +++ /dev/null @@ -1,242 +0,0 @@ -#!/bin/sh -# -unalias -a -export LC_COLLATE=C -########################################################### -## process `blanksky' bkg files for spectra analysis ## -## http://cxc.harvard.edu/ciao/threads/acisbackground/ ## -## ## -## for extracting the spectrum of the background ## -## for determining the background components ## -## ## -## inputs: `evt2_clean', `asol files' ## -## output: `blanksky_c7.fits' for ACIS-S, ## -## `blanksky-c0-3.fits' for ACIS-I ## -## ## -## LIweitiaNux, January 11, 2011 ## -########################################################### - -########################################################### -## ChangeLogs: -## v2: 2012/08/01 -## add ACIS-I support (chips 0-3) -## v3: 2012/08/06, Junhua Gu -## pass parameters by cmd line param -## v4: 2012/08/13, LIweitiaNux -## add `clobber=yes' parameters to CIAO tools -## improve `commandline arguements' -## add `default parameters' -########################################################### - -## about, used in `usage' {{{ -VERSION="v4" -UPDATE="2012-08-14" -## about }}} - -## usage, help {{{ -case "$1" in - -[hH]*|--[hH]*) - printf "usage:\n" - printf " `basename $0` evt= basedir= [ outfile= ]\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`" -# default dir which contains `asols, asol.lis, ...' files -# DFT_BASEDIR="_NOT_EXIST_" -DFT_BASEDIR=".." - -## howto find files in `basedir' -# default `asol.lis pattern' -DFT_ASOLIS_PAT="acis*asol?.lis" -## default parameters }}} - -## 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 -## error code }}} - -## 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 -} - -# reprocess blanksky evt with matched gainfile -blank_regain() { - mv $1 ${1%.fits}_ungain.fits - punlearn acis_process_events - acis_process_events infile="${1%.fits}_ungain.fits" \ - outfile="$1" \ - acaofffile=NONE stop="none" doevtgrade=no \ - apply_cti=yes apply_tgain=no \ - calculate_pi=yes pix_adj=NONE \ - gainfile="$2" \ - eventdef="{s:ccd_id,s:node_id,i:expno,s:chip,s:tdet,f:det,f:sky,s:phas,l:pha,l:pha_ro,f:energy,l:pi,s:fltgrade,s:grade,x:status}" \ - clobber=yes - rm -fv ${1%.fits}_ungain.fits -} -## functions end }}} - -## 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 "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" -# 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" -## 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 files }}} -#exit 0 # test script - -#### main start {{{ -printf "look up coresponding background file ...\n" -BKG_LKP="`acis_bkgrnd_lookup ${EVT}`" -AS_NUM=`echo ${BKG_LKP} | tr ' ' '\n' | \grep 'acis7sD' | wc -l` -AI_NUM=`echo ${BKG_LKP} | tr ' ' '\n' | \grep 'acis[0123]iD' | wc -l` -## determine detector type: ACIS-S / ACIS-I {{{ -if [ ${AS_NUM} -eq 1 ]; then - printf "## ACIS-S, chip: 7\n" - BKG_ROOT="blanksky_c7" - cp -v ${BKG_LKP} ${BKG_ROOT}_orig.fits -elif [ ${AI_NUM} -eq 4 ]; then - printf "## ACIS-I, chip: 0-3\n" - BKG_ROOT="blanksky_c0-3" - AI_FILES="" - for bf in ${BKG_LKP}; do - cp -v ${bf} . - AI_FILES="${AI_FILES},`basename ${bf}`" - done - AI_FILES=${AI_FILES#,} # remove the first ',' - printf "## ACIS-I blanksky files to merge:\n" - printf "## \`${AI_FILES}'\n" - printf "\`dmmerge' to merge the above blanksky files ...\n" - # merge 4 chips blanksky evt files - punlearn dmmerge - dmmerge "${AI_FILES}" ${BKG_ROOT}_orig.fits clobber=yes - rm -fv `echo ${AI_FILES} | tr ',' ' '` # remove original files -else - printf "## ERROR: UNKNOW blanksky files:\n" - printf "## ${BKG_ORIG}\n" - exit ${ERR_BKG} -fi -## determine ACIS type }}} - -## check 'DATMODE' {{{ -## filter blanksky files (status=0) for `VFAINT' observations -DATA_MODE="`dmkeypar ${EVT} DATAMODE echo=yes`" -printf "## DATAMODE: ${DATA_MODE}\n" -if [ "${DATA_MODE}" = "VFAINT" ]; then - mv -fv ${BKG_ROOT}_orig.fits ${BKG_ROOT}_tmp.fits - printf "apply \`status=0' to filter blanksky file ...\n" - punlearn dmcopy - dmcopy "${BKG_ROOT}_tmp.fits[status=0]" ${BKG_ROOT}_orig.fits clobber=yes - rm -fv ${BKG_ROOT}_tmp.fits -fi -## DATAMODE, status=0 }}} - -## check `GAINFILE' of blanksky and evt2 file {{{ -## if NOT match, then reprocess blanksky -GAINFILE_EVT="`dmkeypar ${EVT} GAINFILE echo=yes`" -GAINFILE_BG="`dmkeypar "${BKG_ROOT}_orig.fits" GAINFILE echo=yes`" -if ! [ "${GAINFILE_EVT}" = "${GAINFILE_BG}" ]; then - printf "WARNING: GAINFILE NOT match.\n" - printf "event: ${GAINFILE_EVT}\n" - printf "blank: ${GAINFILE_BG}\n" - printf "reprocess blanksky with evt gainfile ...\n" - # reprocess blanksky using matched evt GAINFILE - GAINFILE="$CALDB/data/chandra/acis/det_gain/`basename ${GAINFILE_EVT}`" - printf "GAINFILE: ${GAINFILE}\n" - blank_regain "${BKG_ROOT}_orig.fits" ${GAINFILE} -fi -## check & match GAINFILE }}} - -printf "add the PNT header keywords ... " -EVT_HEADER="_evt_header.par" -EVT_PNT="_evt_pnt.par" -punlearn dmmakepar -dmmakepar ${EVT} ${EVT_HEADER} clobber=yes -grep -i '_pnt' ${EVT_HEADER} > ${EVT_PNT} -punlearn dmreadpar -dmreadpar ${EVT_PNT} "${BKG_ROOT}_orig.fits[EVENTS]" clobber=yes -printf "DONE\n" - -printf "reproject the background ...\n" -punlearn reproject_events -reproject_events infile=${BKG_ROOT}_orig.fits \ - outfile=${BKG_ROOT}.fits match=${EVT} \ - aspect="@${ASOLIS}" random=0 clobber=yes - -# rename output file if specified -if ! [ -z "${outfile}" ]; then - mv -fv ${BKG_ROOT}.fits ${outfile} -fi -## main end }}} - -# clean -printf "\nclean ...\n" -rm -fv ${BKG_ROOT}_orig.fits ${EVT_PNT} - -printf "\nFINISHED\n" - -# vim: set ts=8 sw=4 tw=0 fenc=utf-8 ft=sh # diff --git a/scripts/ciao_deproj_spectra.sh b/scripts/ciao_deproj_spectra.sh new file mode 100755 index 0000000..abebc7c --- /dev/null +++ b/scripts/ciao_deproj_spectra.sh @@ -0,0 +1,684 @@ +#!/bin/sh +# +trap date INT +unalias -a +export GREP_OPTIONS="" +export LC_COLLATE=C +# +########################################################### +## generate `src' and `bkg' spectra as well as ## +## RMF/ARFs using `specextract' ## +## for radial spectra analysis ## +## to get temperature/abundance profile ## +## XSPEC model `projct' for deprojection analysis ## +## ## +## Ref: Chandra spectrum analysis ## +## http://cxc.harvard.edu/ciao/threads/extended/ ## +## Ref: specextract ## +## http://cxc.harvard.edu/ciao/ahelp/specextract.html ## +## Ref: CIAO v4.4 region bugs ## +## http://cxc.harvard.edu/ciao/bugs/regions.html#bug-12187 +## ## +## LIweitiaNux, July 24, 2012 ## +########################################################### + +########################################################### +## ChangeLogs: +## v5, 2012/08/05 +## XFLT0005 not modified as pie end angle +## add `background renormalization' +## v6, 2012/08/08, Gu Junhua +## Modified to using config file to pass parameters +## Use grppha to rebin the spectrum +## v7, 2012/08/10, LIweitiaNux +## account blanksky, local bkg, specified bkg +## change name to `ciao_deproj_spectra_v*.sh' +## add `error status' +## Imporve ${CFG_FILE} +## Imporve comments +## 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' {{{ +VERSION="v8" +UPDATE="2012-08-14" +## about }}} + +## usage, help {{{ +case "$1" in + -[hH]*|--[hH]*) + printf "usage:\n" + printf " `basename $0` evt= reg= bkgd= basedir= nh= z= [ grpcmd= log= ]\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`" +# default `radial region file' +#DFT_REG_IN="_NOT_EXIST_" +DFT_REG_IN="rspec.reg" +# default dir which contains `asols, asol.lis, ...' files +# DFT_BASEDIR=".." +DFT_BASEDIR="_NOT_EXIST_" +# 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 1 128 4 129 256 8 257 512 16 513 1024 32" +DFT_GRP_CMD="group min 20" +# default `log file' +DFT_LOGFILE="deproj_spectra_`date '+%Y%m%d'`.log" + +# default output xspec scripts +XSPEC_DEPROJ="xspec_deproj.xcm" +XSPEC_PROJTD="xspec_projected.xcm" + +## 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 `pbk file pattern' +DFT_PBK_PAT="acis*pbk?.fits" +# default `msk file pattern' +DFT_MSK_PAT="acis*msk?.fits" + +## abundance standard +ABUND="grsa" +## default parameters }}} + +## 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 +## error code }}} + +## 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 +} + +## background renormalization (BACKSCAL) {{{ +# renorm background according to particle background +# energy range: 9.5-12.0 keV (channel: 651-822) +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 }'` + punlearn dmkeypar + EXPTIME=`dmkeypar $1 EXPOSURE echo=yes` + BACK=`dmkeypar $1 BACKSCAL echo=yes` + # fix `scientific notation' bug for `bc' + EXPTIME_B=`echo ${EXPTIME} | sed 's/[eE]/\*10\^/' | sed 's/+//'` + BACK_B=`echo "( ${BACK} )" | sed 's/[eE]/\*10\^/' | sed 's/+//'` + PB_FLUX=`echo "scale = 16; ${COUNTS} / ${EXPTIME_B} / ${BACK_B}" | bc -l` + echo ${PB_FLUX} +} + +bkg_renorm() { + # $1: src spectrum, $2: back spectrum + PBFLUX_SRC=`pb_flux $1` + PBFLUX_BKG=`pb_flux $2` + BACK_OLD=`dmkeypar $2 BACKSCAL echo=yes` + BACK_OLD_B=`echo "( ${BACK_OLD} )" | sed 's/[eE]/\*10\^/' | sed 's/+//'` + BACK_NEW=`echo "scale = 16; ${BACK_OLD_B} * ${PBFLUX_BKG} / ${PBFLUX_SRC}" | bc -l` + printf "\`$2': BACKSCAL:\n" + printf " ${BACK_OLD} --> ${BACK_NEW}\n" + punlearn dmhedit + dmhedit infile=$2 filelist=none operation=add \ + key=BACKSCAL value=${BACK_NEW} comment="old value: ${BACK_OLD}" +} +## bkg renorm }}} +## functions end }}} + +## 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 + REG_IN="${reg}" +elif [ -r "${DFT_REG_IN}" ]; then + REG_IN=${DFT_REG_IN} +else + read -p "> radial spec region file: " REG_IN + if [ ! -r "${REG_IN}" ]; then + printf "ERROR: cannot access given \`${REG_IN}' region file\n" + exit ${ERR_REG} + fi +fi +printf "## use radial reg: \`${REG_IN}'\n" | ${TOLOG} +# check given bkgd, determine background {{{ +if [ -z "${bkgd}" ]; then + read -p "> background (blanksky_evt | lbkg_reg | bkg_spec): " BKGD +else + BKGD=${bkgd} +fi +if [ ! -r "${BKGD}" ]; then + printf "ERROR: cannot access given \`${BKGD}'\n" + exit ${ERR_BKG} +fi +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 + printf "## given \`${BKGD}' is a \`text file'\n" + printf "## use it as local bkg region file\n" + printf "## use *LOCAL BKG SPEC*\n" | ${TOLOG} + # just set flags, extract spectrum later + USE_LBKG_REG=YES + USE_BLANKSKY=NO + USE_BKG_SPEC=NO +elif file -bL ${BKGD} | grep -qi 'FITS'; then + printf "## given \`${BKGD}' is a \`FITS file'\n" + # get FITS header keyword + HDUCLAS1=`dmkeypar ${BKGD} HDUCLAS1 echo=yes` + if [ "${HDUCLAS1}" = "EVENTS" ]; then + # event file + printf "## given file is \`event'\n" + # check if `blanksky' or `stowed bkg' + BKG_OBJ=`dmkeypar ${BKGD} OBJECT echo=yes` + if [ "${BKG_OBJ}" = "BACKGROUND DATASET" ] || [ "${BKG_OBJ}" = "ACIS STOWED" ]; then + # valid bkg evt file + printf "## given FITS file is a valid bkgrnd file\n" + printf "## use *BLANKSKY*\n" | ${TOLOG} + USE_BLANKSKY=YES + USE_LBKG_REG=NO + USE_BKG_SPEC=NO + # specify `BLANKSKY' + BLANKSKY=${BKGD} + else + # invalid bkg evt file + printf "ERROR: invalid bkg evt file given\n" + exit ${ERR_BKGTY} + fi + elif [ "${HDUCLAS1}" = "SPECTRUM" ]; then + # spectrum file + printf "## given file is \`spectrum'\n" + printf "## use *BKG SPECTRUM*\n" | ${TOLOG} + USE_BKG_SPEC=YES + USE_BLANKSKY=NO + USE_LBKG_REG=NO + # specify `BKG_SPEC' + BKG_SPEC=${BKGD} + else + # other type + printf "ERROR: other type FITS given\n" + exit ${ERR_BKGTY} + fi +else + printf "ERROR: given \`${BKGD}' type UNKNOWN\n" + exit ${ERR_BKGTY} +fi +# bkgd }}} +# 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 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 given `grpcmd' +if [ ! -z "${grpcmd}" ]; then + GRP_CMD="${grpcmd}" +else + GRP_CMD="${DFT_GRP_CMD}" +fi +printf "## use grppha cmd: \`${GRP_CMD}'\n" | ${TOLOG} +# rootname for output files +[ "x${ROOTNAME}" = "x" ] && ROOTNAME="${REG_IN%.reg}" +printf "## use rootname: \`${ROOTNAME}'\n" | ${TOLOG} +## parameters }}} + +## check needed files {{{ +# check the validity of *pie* regions +printf "check pie reg validity ...\n" +INVALID=`cat ${REG_IN} | grep -i 'pie' | awk -F, '{ print $6 }' | tr -d ')' | awk '$1 > 360'` +if [ "x${INVALID}" != "x" ]; then + printf "WARNING: some pie region's END_ANGLE > 360\n" | ${TOLOG} + printf " CIAO v4.4 tools may run into trouble\n" +fi + +# check files in `basedir' +printf "check needed files in basedir \`${BASEDIR}' ...\n" +# check asolis 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" | ${TOLOG} +# 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 pbk file +PBK=`ls -1 ${BASEDIR}/${DFT_PBK_PAT} | head -n 1` +if [ -z "${PBK}" ]; then + printf "ERROR: cannot find \"${DFT_PBK_PAT}\" in dir \`${BASEDIR}'\n" + exit ${ERR_PBK} +fi +printf "## use pbk: \`${PBK}'\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: \`${MSK}'\n" | ${TOLOG} +## check files }}} + +## process local background {{{ +if [ "${USE_LBKG_REG}" = "YES" ]; then + BKG_EVT=${EVT} + LBKG_REG=${BKGD} + LBKG_REG_CIAO="_tmp_${LBKG_REG}" + cp -fv ${LBKG_REG} ${LBKG_REG_CIAO} + ## check background region (CIAO v4.4 bug) {{{ + printf "check background region ...\n" + INVALID=`grep -i 'pie' ${LBKG_REG_CIAO} | awk -F, '{ print $6 }' | tr -d ')' | awk '$1 > 360'` + if [ "x${INVALID}" != "x" ]; then + printf "WARNING: fix for pie region:\n" + cat ${LBKG_REG_CIAO} + for angle in ${INVALID}; do + A_OLD=`echo ${angle} | sed 's/\./\\\./'` + A_NEW=`echo ${angle}-360 | bc -l | sed 's/\./\\\./'` + sed -i'' "s/${A_OLD}\ *)/${A_NEW})/" ${LBKG_REG_CIAO} + done + printf " -->" + cat ${LBKG_REG_CIAO} + printf "======================\n" + fi + ## background region }}} + ## extract local background spectrum + printf "extract local background spectrum ...\n" + BKG_SPEC="${LBKG_REG%.reg}.pi" + punlearn dmextract + dmextract infile="${BKG_EVT}[sky=region(${LBKG_REG_CIAO})][bin pi]" \ + outfile=${BKG_SPEC} wmap="[bin det=8]" clobber=yes + rm -fv ${LBKG_REG_CIAO} + printf "renormalizing the spectrum later ...\n" +fi +## local bkg }}} + +## modify the region file, remove the commented and blank lines {{{ +REG_NEW="_new.reg" +REG_TMP="_tmp.reg" +[ -f ${REG_NEW} ] && rm -fv ${REG_NEW} +[ -f ${REG_TMP} ] && rm -fv ${REG_TMP} +cat ${REG_IN} | sed 's/#.*$//' | grep -Ev '^\s*$' > ${REG_NEW} +## REG_IN }}} + +## `specextract' to extract spectrum {{{ +LINES="`wc -l ${REG_NEW} | cut -d' ' -f1`" +printf "\n======================================\n" +printf "TOTAL *${LINES}* regions to process ......\n" +for i in `seq ${LINES}`; do + printf "\n==============================\n" + printf ">>> PROCESS REGION ${i} ...\n" + + ## generate corresponding `region' file + rm -f ${REG_TMP} ${REG_CIAO} 2> /dev/null + head -n ${i} ${REG_NEW} | tail -n 1 > ${REG_TMP} + REG_CIAO="${REG_TMP%.reg}_ciao.reg" + cp -fv ${REG_TMP} ${REG_CIAO} + ## check the validity of *pie* regions {{{ + INVALID=`grep -i 'pie' ${REG_TMP} | awk -F, '{ print $6 }' | tr -d ')' | awk '$1 > 360'` + if [ "x${INVALID}" != "x" ]; then + printf "WARNING: fix for pie region:\n" + cat ${REG_CIAO} + A_OLD=`echo ${INVALID} | sed 's/\./\\\./'` + A_NEW=`echo ${INVALID}-360 | bc -l | sed 's/\./\\\./'` + sed -i'' "s/${A_OLD}\ *)/${A_NEW})/" ${REG_CIAO} + printf " -->\n" + cat ${REG_CIAO} + fi + # check pie region }}} + + ## use `specextract' to extract spectra {{{ + # NOTE: set `binarfwmap=2' to save the time for generating `ARF' + # I have tested that this bin factor has little impact on the results. + # 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}" \ + ${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 \ + grouptype=NONE binspec=NONE \ + clobber=yes verbose=2 + ## specextract }}} + + RSPEC_PI="r${i}_${ROOTNAME}.pi" + RSPEC_BKG_PI="${RSPEC_PI%.pi}_bkg.pi" + ## background spectrum {{{ + ## generate the blanksky bkg spectrum by self + if [ "${USE_BLANKSKY}" = "YES" ]; then + # use blanksky as background file + printf "extract blanksky bkg spectrum ...\n" + punlearn dmextract + dmextract infile="${BLANKSKY}[sky=region(${REG_CIAO})][bin pi]" \ + outfile=${RSPEC_BKG_PI} wmap="[bin det=8]" clobber=yes + elif [ "${USE_LBKG_REG}" = "YES" ] || [ "${USE_BKG_SPEC}" = "YES" ]; then + # use *local background* or specified background spectrum + cp -fv ${BKG_SPEC} ${RSPEC_BKG_PI} + fi + ## background }}} + + ## bkg renormalization {{{ + printf "Renormalize background ...\n" + bkg_renorm ${RSPEC_PI} ${RSPEC_BKG_PI} + ## bkg renorm }}} + + ## group spectrum {{{ + printf "group spectrum \`${RSPEC_PI}' using \`grppha'\n" + RSPEC_GRP_PI="${RSPEC_PI%.pi}_grp.pi" + grppha infile="${RSPEC_PI}" outfile="${RSPEC_GRP_PI}" \ + comm="${GRP_CMD} & exit" clobber=yes > /dev/null + ## group }}} + + ## `XFLT####' keywords for XSPEC model `projct' {{{ + printf "update file headers ...\n" + punlearn dmhedit + if grep -qi 'pie' ${REG_TMP}; then + R_OUT="`awk -F',' '{ print $4 }' ${REG_TMP} | tr -d ')'`" + A_BEGIN="`awk -F',' '{ print $5 }' ${REG_TMP} | tr -d ')'`" + A_END="`awk -F',' '{ print $6 }' ${REG_TMP} | tr -d ')'`" + # RSPEC_PI + dmhedit infile=${RSPEC_PI} filelist=NONE operation=add \ + key="XFLT0001" value=${R_OUT} + dmhedit infile=${RSPEC_PI} filelist=NONE operation=add \ + key="XFLT0002" value=${R_OUT} + dmhedit infile=${RSPEC_PI} filelist=NONE operation=add \ + key="XFLT0003" value=0 + dmhedit infile=${RSPEC_PI} filelist=NONE operation=add \ + key="XFLT0004" value=${A_BEGIN} + dmhedit infile=${RSPEC_PI} filelist=NONE operation=add \ + key="XFLT0005" value=${A_END} + # RSPEC_GRP_PI + dmhedit infile=${RSPEC_GRP_PI} filelist=NONE operation=add \ + key="XFLT0001" value=${R_OUT} + dmhedit infile=${RSPEC_GRP_PI} filelist=NONE operation=add \ + key="XFLT0002" value=${R_OUT} + dmhedit infile=${RSPEC_GRP_PI} filelist=NONE operation=add \ + key="XFLT0003" value=0 + dmhedit infile=${RSPEC_GRP_PI} filelist=NONE operation=add \ + key="XFLT0004" value=${A_BEGIN} + dmhedit infile=${RSPEC_GRP_PI} filelist=NONE operation=add \ + key="XFLT0005" value=${A_END} + elif grep -qi 'annulus' ${REG_TMP}; then + R_OUT="`awk -F',' '{ print $4 }' ${REG_TMP} | tr -d ')'`" + # RSPEC_PI + dmhedit infile=${RSPEC_PI} filelist=NONE operation=add \ + key="XFLT0001" value=${R_OUT} + dmhedit infile=${RSPEC_PI} filelist=NONE operation=add \ + key="XFLT0002" value=${R_OUT} + dmhedit infile=${RSPEC_PI} filelist=NONE operation=add \ + key="XFLT0003" value=0 + # RSPEC_GRP_PI + dmhedit infile=${RSPEC_GRP_PI} filelist=NONE operation=add \ + key="XFLT0001" value=${R_OUT} + dmhedit infile=${RSPEC_GRP_PI} filelist=NONE operation=add \ + key="XFLT0002" value=${R_OUT} + dmhedit infile=${RSPEC_GRP_PI} filelist=NONE operation=add \ + key="XFLT0003" value=0 + else + printf "*** WARNING: region file NOT MATCH!!\n" + fi + ## `XFLT####' }}} + +done # end *for*, `specextract' +## `specextract' }}} + +## clean +printf "clean ...\n" +rm -f ${REG_TMP} ${REG_CIAO} 2> /dev/null + + +########################################################### +## generate a script file for XSPEC ## +########################################################### +printf "generate scripts for XSPEC ...\n" +## xspec script (deproj) {{{ +printf "XSPEC script for deprojection analysis\n" +[ -e ${XSPEC_DEPROJ} ] && rm -fv ${XSPEC_DEPROJ} + +cat >> ${XSPEC_DEPROJ} << _EOF_ +# +# XSPEC +# radial spectra (deprojection analysis) +# model projct*wabs*apec +# +# generated by script: \``basename $0`' +# `date` + +statistic chi + +# load data +_EOF_ + +for i in `seq ${LINES}`; do + RSPEC="r${i}_${ROOTNAME}" + cat >> ${XSPEC_DEPROJ} << _EOF_ +data ${i}:${i} ${RSPEC}_grp.pi +response 1:${i} ${RSPEC}.wrmf +arf 1:${i} ${RSPEC}.warf +backgrnd ${i} ${RSPEC}_bkg.pi + +_EOF_ +done + +cat >> ${XSPEC_DEPROJ} << _EOF_ + +# filter needed energy range +ignore bad +ignore **:0.0-0.7,7.0-** + +method leven 1000 0.01 +# change abundance standard +abund ${ABUND} +xsect bcmc +cosmo 70 0 0.73 +xset delta 0.01 +systematic 0 +# auto answer +query yes + +# plot related +setplot energy + +# model to use +model projct*wabs*apec + 0 + 0 + 0 + ${N_H} -0.001 0 0 100000 1e+06 + 1 0.01 0.008 0.008 64 64 + 0.5 0.001 0 0 5 5 + ${REDSHIFT} -0.01 -0.999 -0.999 10 10 + 1 0.01 0 0 1e+24 1e+24 +_EOF_ + +INPUT_TIMES=`expr ${LINES} - 1` +for i in `seq ${INPUT_TIMES}`; do + cat >> ${XSPEC_DEPROJ} << _EOF_ += 1 += 2 += 3 += 4 + 1 0.01 0.008 0.008 64 64 + 0.5 0.001 0 0 5 5 += 7 + 1 0.01 0 0 1e+24 1e+24 +_EOF_ +done +## xspec script }}} + +########################################################### +## xspec script (projected) {{{ +printf "XSPEC script for projected analysis\n" +[ -e ${XSPEC_PROJTD} ] && rm -fv ${XSPEC_PROJTD} + +cat >> ${XSPEC_PROJTD} << _EOF_ +# +# XSPEC +# radial spectra (projected analysis) +# model wabs*apec +# +# generated by script: \``basename $0`' +# `date` + +statistic chi + +# load data +_EOF_ + +for i in `seq ${LINES}`; do + RSPEC="r${i}_${ROOTNAME}" + cat >> ${XSPEC_PROJTD} << _EOF_ +data ${i}:${i} ${RSPEC}_grp.pi +response 1:${i} ${RSPEC}.wrmf +arf 1:${i} ${RSPEC}.warf +backgrnd ${i} ${RSPEC}_bkg.pi + +_EOF_ +done + +cat >> ${XSPEC_PROJTD} << _EOF_ + +# filter needed energy range +ignore bad +ignore **:0.0-0.7,7.0-** + +method leven 1000 0.01 +# change abundance standard +abund ${ABUND} +xsect bcmc +cosmo 70 0 0.73 +xset delta 0.01 +systematic 0 +# auto answer +query yes + +# plot related +setplot energy + +# model to use +model wabs*apec + ${N_H} -0.001 0 0 100000 1e+06 + 1 0.01 0.008 0.008 64 64 + 0.5 0.001 0 0 5 5 + ${REDSHIFT} -0.01 -0.999 -0.999 10 10 + 1 0.01 0 0 1e+24 1e+24 +_EOF_ + +INPUT_TIMES=`expr ${LINES} - 1` +for i in `seq ${INPUT_TIMES}`; do + cat >> ${XSPEC_PROJTD} << _EOF_ += 1 + 1 0.01 0.008 0.008 64 64 + 0.5 0.001 0 0 5 5 += 4 + 1 0.01 0 0 1e+24 1e+24 +_EOF_ +done + +printf "DONE\n" +## xspec script }}} +########################################################### + +printf "ALL FINISHED\n" + +# vim: set ts=8 sw=4 tw=0 fenc=utf-8 ft=sh: # diff --git a/scripts/ciao_deproj_spectra_v8.sh b/scripts/ciao_deproj_spectra_v8.sh deleted file mode 100755 index abebc7c..0000000 --- a/scripts/ciao_deproj_spectra_v8.sh +++ /dev/null @@ -1,684 +0,0 @@ -#!/bin/sh -# -trap date INT -unalias -a -export GREP_OPTIONS="" -export LC_COLLATE=C -# -########################################################### -## generate `src' and `bkg' spectra as well as ## -## RMF/ARFs using `specextract' ## -## for radial spectra analysis ## -## to get temperature/abundance profile ## -## XSPEC model `projct' for deprojection analysis ## -## ## -## Ref: Chandra spectrum analysis ## -## http://cxc.harvard.edu/ciao/threads/extended/ ## -## Ref: specextract ## -## http://cxc.harvard.edu/ciao/ahelp/specextract.html ## -## Ref: CIAO v4.4 region bugs ## -## http://cxc.harvard.edu/ciao/bugs/regions.html#bug-12187 -## ## -## LIweitiaNux, July 24, 2012 ## -########################################################### - -########################################################### -## ChangeLogs: -## v5, 2012/08/05 -## XFLT0005 not modified as pie end angle -## add `background renormalization' -## v6, 2012/08/08, Gu Junhua -## Modified to using config file to pass parameters -## Use grppha to rebin the spectrum -## v7, 2012/08/10, LIweitiaNux -## account blanksky, local bkg, specified bkg -## change name to `ciao_deproj_spectra_v*.sh' -## add `error status' -## Imporve ${CFG_FILE} -## Imporve comments -## 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' {{{ -VERSION="v8" -UPDATE="2012-08-14" -## about }}} - -## usage, help {{{ -case "$1" in - -[hH]*|--[hH]*) - printf "usage:\n" - printf " `basename $0` evt= reg= bkgd= basedir= nh= z= [ grpcmd= log= ]\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`" -# default `radial region file' -#DFT_REG_IN="_NOT_EXIST_" -DFT_REG_IN="rspec.reg" -# default dir which contains `asols, asol.lis, ...' files -# DFT_BASEDIR=".." -DFT_BASEDIR="_NOT_EXIST_" -# 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 1 128 4 129 256 8 257 512 16 513 1024 32" -DFT_GRP_CMD="group min 20" -# default `log file' -DFT_LOGFILE="deproj_spectra_`date '+%Y%m%d'`.log" - -# default output xspec scripts -XSPEC_DEPROJ="xspec_deproj.xcm" -XSPEC_PROJTD="xspec_projected.xcm" - -## 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 `pbk file pattern' -DFT_PBK_PAT="acis*pbk?.fits" -# default `msk file pattern' -DFT_MSK_PAT="acis*msk?.fits" - -## abundance standard -ABUND="grsa" -## default parameters }}} - -## 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 -## error code }}} - -## 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 -} - -## background renormalization (BACKSCAL) {{{ -# renorm background according to particle background -# energy range: 9.5-12.0 keV (channel: 651-822) -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 }'` - punlearn dmkeypar - EXPTIME=`dmkeypar $1 EXPOSURE echo=yes` - BACK=`dmkeypar $1 BACKSCAL echo=yes` - # fix `scientific notation' bug for `bc' - EXPTIME_B=`echo ${EXPTIME} | sed 's/[eE]/\*10\^/' | sed 's/+//'` - BACK_B=`echo "( ${BACK} )" | sed 's/[eE]/\*10\^/' | sed 's/+//'` - PB_FLUX=`echo "scale = 16; ${COUNTS} / ${EXPTIME_B} / ${BACK_B}" | bc -l` - echo ${PB_FLUX} -} - -bkg_renorm() { - # $1: src spectrum, $2: back spectrum - PBFLUX_SRC=`pb_flux $1` - PBFLUX_BKG=`pb_flux $2` - BACK_OLD=`dmkeypar $2 BACKSCAL echo=yes` - BACK_OLD_B=`echo "( ${BACK_OLD} )" | sed 's/[eE]/\*10\^/' | sed 's/+//'` - BACK_NEW=`echo "scale = 16; ${BACK_OLD_B} * ${PBFLUX_BKG} / ${PBFLUX_SRC}" | bc -l` - printf "\`$2': BACKSCAL:\n" - printf " ${BACK_OLD} --> ${BACK_NEW}\n" - punlearn dmhedit - dmhedit infile=$2 filelist=none operation=add \ - key=BACKSCAL value=${BACK_NEW} comment="old value: ${BACK_OLD}" -} -## bkg renorm }}} -## functions end }}} - -## 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 - REG_IN="${reg}" -elif [ -r "${DFT_REG_IN}" ]; then - REG_IN=${DFT_REG_IN} -else - read -p "> radial spec region file: " REG_IN - if [ ! -r "${REG_IN}" ]; then - printf "ERROR: cannot access given \`${REG_IN}' region file\n" - exit ${ERR_REG} - fi -fi -printf "## use radial reg: \`${REG_IN}'\n" | ${TOLOG} -# check given bkgd, determine background {{{ -if [ -z "${bkgd}" ]; then - read -p "> background (blanksky_evt | lbkg_reg | bkg_spec): " BKGD -else - BKGD=${bkgd} -fi -if [ ! -r "${BKGD}" ]; then - printf "ERROR: cannot access given \`${BKGD}'\n" - exit ${ERR_BKG} -fi -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 - printf "## given \`${BKGD}' is a \`text file'\n" - printf "## use it as local bkg region file\n" - printf "## use *LOCAL BKG SPEC*\n" | ${TOLOG} - # just set flags, extract spectrum later - USE_LBKG_REG=YES - USE_BLANKSKY=NO - USE_BKG_SPEC=NO -elif file -bL ${BKGD} | grep -qi 'FITS'; then - printf "## given \`${BKGD}' is a \`FITS file'\n" - # get FITS header keyword - HDUCLAS1=`dmkeypar ${BKGD} HDUCLAS1 echo=yes` - if [ "${HDUCLAS1}" = "EVENTS" ]; then - # event file - printf "## given file is \`event'\n" - # check if `blanksky' or `stowed bkg' - BKG_OBJ=`dmkeypar ${BKGD} OBJECT echo=yes` - if [ "${BKG_OBJ}" = "BACKGROUND DATASET" ] || [ "${BKG_OBJ}" = "ACIS STOWED" ]; then - # valid bkg evt file - printf "## given FITS file is a valid bkgrnd file\n" - printf "## use *BLANKSKY*\n" | ${TOLOG} - USE_BLANKSKY=YES - USE_LBKG_REG=NO - USE_BKG_SPEC=NO - # specify `BLANKSKY' - BLANKSKY=${BKGD} - else - # invalid bkg evt file - printf "ERROR: invalid bkg evt file given\n" - exit ${ERR_BKGTY} - fi - elif [ "${HDUCLAS1}" = "SPECTRUM" ]; then - # spectrum file - printf "## given file is \`spectrum'\n" - printf "## use *BKG SPECTRUM*\n" | ${TOLOG} - USE_BKG_SPEC=YES - USE_BLANKSKY=NO - USE_LBKG_REG=NO - # specify `BKG_SPEC' - BKG_SPEC=${BKGD} - else - # other type - printf "ERROR: other type FITS given\n" - exit ${ERR_BKGTY} - fi -else - printf "ERROR: given \`${BKGD}' type UNKNOWN\n" - exit ${ERR_BKGTY} -fi -# bkgd }}} -# 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 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 given `grpcmd' -if [ ! -z "${grpcmd}" ]; then - GRP_CMD="${grpcmd}" -else - GRP_CMD="${DFT_GRP_CMD}" -fi -printf "## use grppha cmd: \`${GRP_CMD}'\n" | ${TOLOG} -# rootname for output files -[ "x${ROOTNAME}" = "x" ] && ROOTNAME="${REG_IN%.reg}" -printf "## use rootname: \`${ROOTNAME}'\n" | ${TOLOG} -## parameters }}} - -## check needed files {{{ -# check the validity of *pie* regions -printf "check pie reg validity ...\n" -INVALID=`cat ${REG_IN} | grep -i 'pie' | awk -F, '{ print $6 }' | tr -d ')' | awk '$1 > 360'` -if [ "x${INVALID}" != "x" ]; then - printf "WARNING: some pie region's END_ANGLE > 360\n" | ${TOLOG} - printf " CIAO v4.4 tools may run into trouble\n" -fi - -# check files in `basedir' -printf "check needed files in basedir \`${BASEDIR}' ...\n" -# check asolis 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" | ${TOLOG} -# 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 pbk file -PBK=`ls -1 ${BASEDIR}/${DFT_PBK_PAT} | head -n 1` -if [ -z "${PBK}" ]; then - printf "ERROR: cannot find \"${DFT_PBK_PAT}\" in dir \`${BASEDIR}'\n" - exit ${ERR_PBK} -fi -printf "## use pbk: \`${PBK}'\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: \`${MSK}'\n" | ${TOLOG} -## check files }}} - -## process local background {{{ -if [ "${USE_LBKG_REG}" = "YES" ]; then - BKG_EVT=${EVT} - LBKG_REG=${BKGD} - LBKG_REG_CIAO="_tmp_${LBKG_REG}" - cp -fv ${LBKG_REG} ${LBKG_REG_CIAO} - ## check background region (CIAO v4.4 bug) {{{ - printf "check background region ...\n" - INVALID=`grep -i 'pie' ${LBKG_REG_CIAO} | awk -F, '{ print $6 }' | tr -d ')' | awk '$1 > 360'` - if [ "x${INVALID}" != "x" ]; then - printf "WARNING: fix for pie region:\n" - cat ${LBKG_REG_CIAO} - for angle in ${INVALID}; do - A_OLD=`echo ${angle} | sed 's/\./\\\./'` - A_NEW=`echo ${angle}-360 | bc -l | sed 's/\./\\\./'` - sed -i'' "s/${A_OLD}\ *)/${A_NEW})/" ${LBKG_REG_CIAO} - done - printf " -->" - cat ${LBKG_REG_CIAO} - printf "======================\n" - fi - ## background region }}} - ## extract local background spectrum - printf "extract local background spectrum ...\n" - BKG_SPEC="${LBKG_REG%.reg}.pi" - punlearn dmextract - dmextract infile="${BKG_EVT}[sky=region(${LBKG_REG_CIAO})][bin pi]" \ - outfile=${BKG_SPEC} wmap="[bin det=8]" clobber=yes - rm -fv ${LBKG_REG_CIAO} - printf "renormalizing the spectrum later ...\n" -fi -## local bkg }}} - -## modify the region file, remove the commented and blank lines {{{ -REG_NEW="_new.reg" -REG_TMP="_tmp.reg" -[ -f ${REG_NEW} ] && rm -fv ${REG_NEW} -[ -f ${REG_TMP} ] && rm -fv ${REG_TMP} -cat ${REG_IN} | sed 's/#.*$//' | grep -Ev '^\s*$' > ${REG_NEW} -## REG_IN }}} - -## `specextract' to extract spectrum {{{ -LINES="`wc -l ${REG_NEW} | cut -d' ' -f1`" -printf "\n======================================\n" -printf "TOTAL *${LINES}* regions to process ......\n" -for i in `seq ${LINES}`; do - printf "\n==============================\n" - printf ">>> PROCESS REGION ${i} ...\n" - - ## generate corresponding `region' file - rm -f ${REG_TMP} ${REG_CIAO} 2> /dev/null - head -n ${i} ${REG_NEW} | tail -n 1 > ${REG_TMP} - REG_CIAO="${REG_TMP%.reg}_ciao.reg" - cp -fv ${REG_TMP} ${REG_CIAO} - ## check the validity of *pie* regions {{{ - INVALID=`grep -i 'pie' ${REG_TMP} | awk -F, '{ print $6 }' | tr -d ')' | awk '$1 > 360'` - if [ "x${INVALID}" != "x" ]; then - printf "WARNING: fix for pie region:\n" - cat ${REG_CIAO} - A_OLD=`echo ${INVALID} | sed 's/\./\\\./'` - A_NEW=`echo ${INVALID}-360 | bc -l | sed 's/\./\\\./'` - sed -i'' "s/${A_OLD}\ *)/${A_NEW})/" ${REG_CIAO} - printf " -->\n" - cat ${REG_CIAO} - fi - # check pie region }}} - - ## use `specextract' to extract spectra {{{ - # NOTE: set `binarfwmap=2' to save the time for generating `ARF' - # I have tested that this bin factor has little impact on the results. - # 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}" \ - ${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 \ - grouptype=NONE binspec=NONE \ - clobber=yes verbose=2 - ## specextract }}} - - RSPEC_PI="r${i}_${ROOTNAME}.pi" - RSPEC_BKG_PI="${RSPEC_PI%.pi}_bkg.pi" - ## background spectrum {{{ - ## generate the blanksky bkg spectrum by self - if [ "${USE_BLANKSKY}" = "YES" ]; then - # use blanksky as background file - printf "extract blanksky bkg spectrum ...\n" - punlearn dmextract - dmextract infile="${BLANKSKY}[sky=region(${REG_CIAO})][bin pi]" \ - outfile=${RSPEC_BKG_PI} wmap="[bin det=8]" clobber=yes - elif [ "${USE_LBKG_REG}" = "YES" ] || [ "${USE_BKG_SPEC}" = "YES" ]; then - # use *local background* or specified background spectrum - cp -fv ${BKG_SPEC} ${RSPEC_BKG_PI} - fi - ## background }}} - - ## bkg renormalization {{{ - printf "Renormalize background ...\n" - bkg_renorm ${RSPEC_PI} ${RSPEC_BKG_PI} - ## bkg renorm }}} - - ## group spectrum {{{ - printf "group spectrum \`${RSPEC_PI}' using \`grppha'\n" - RSPEC_GRP_PI="${RSPEC_PI%.pi}_grp.pi" - grppha infile="${RSPEC_PI}" outfile="${RSPEC_GRP_PI}" \ - comm="${GRP_CMD} & exit" clobber=yes > /dev/null - ## group }}} - - ## `XFLT####' keywords for XSPEC model `projct' {{{ - printf "update file headers ...\n" - punlearn dmhedit - if grep -qi 'pie' ${REG_TMP}; then - R_OUT="`awk -F',' '{ print $4 }' ${REG_TMP} | tr -d ')'`" - A_BEGIN="`awk -F',' '{ print $5 }' ${REG_TMP} | tr -d ')'`" - A_END="`awk -F',' '{ print $6 }' ${REG_TMP} | tr -d ')'`" - # RSPEC_PI - dmhedit infile=${RSPEC_PI} filelist=NONE operation=add \ - key="XFLT0001" value=${R_OUT} - dmhedit infile=${RSPEC_PI} filelist=NONE operation=add \ - key="XFLT0002" value=${R_OUT} - dmhedit infile=${RSPEC_PI} filelist=NONE operation=add \ - key="XFLT0003" value=0 - dmhedit infile=${RSPEC_PI} filelist=NONE operation=add \ - key="XFLT0004" value=${A_BEGIN} - dmhedit infile=${RSPEC_PI} filelist=NONE operation=add \ - key="XFLT0005" value=${A_END} - # RSPEC_GRP_PI - dmhedit infile=${RSPEC_GRP_PI} filelist=NONE operation=add \ - key="XFLT0001" value=${R_OUT} - dmhedit infile=${RSPEC_GRP_PI} filelist=NONE operation=add \ - key="XFLT0002" value=${R_OUT} - dmhedit infile=${RSPEC_GRP_PI} filelist=NONE operation=add \ - key="XFLT0003" value=0 - dmhedit infile=${RSPEC_GRP_PI} filelist=NONE operation=add \ - key="XFLT0004" value=${A_BEGIN} - dmhedit infile=${RSPEC_GRP_PI} filelist=NONE operation=add \ - key="XFLT0005" value=${A_END} - elif grep -qi 'annulus' ${REG_TMP}; then - R_OUT="`awk -F',' '{ print $4 }' ${REG_TMP} | tr -d ')'`" - # RSPEC_PI - dmhedit infile=${RSPEC_PI} filelist=NONE operation=add \ - key="XFLT0001" value=${R_OUT} - dmhedit infile=${RSPEC_PI} filelist=NONE operation=add \ - key="XFLT0002" value=${R_OUT} - dmhedit infile=${RSPEC_PI} filelist=NONE operation=add \ - key="XFLT0003" value=0 - # RSPEC_GRP_PI - dmhedit infile=${RSPEC_GRP_PI} filelist=NONE operation=add \ - key="XFLT0001" value=${R_OUT} - dmhedit infile=${RSPEC_GRP_PI} filelist=NONE operation=add \ - key="XFLT0002" value=${R_OUT} - dmhedit infile=${RSPEC_GRP_PI} filelist=NONE operation=add \ - key="XFLT0003" value=0 - else - printf "*** WARNING: region file NOT MATCH!!\n" - fi - ## `XFLT####' }}} - -done # end *for*, `specextract' -## `specextract' }}} - -## clean -printf "clean ...\n" -rm -f ${REG_TMP} ${REG_CIAO} 2> /dev/null - - -########################################################### -## generate a script file for XSPEC ## -########################################################### -printf "generate scripts for XSPEC ...\n" -## xspec script (deproj) {{{ -printf "XSPEC script for deprojection analysis\n" -[ -e ${XSPEC_DEPROJ} ] && rm -fv ${XSPEC_DEPROJ} - -cat >> ${XSPEC_DEPROJ} << _EOF_ -# -# XSPEC -# radial spectra (deprojection analysis) -# model projct*wabs*apec -# -# generated by script: \``basename $0`' -# `date` - -statistic chi - -# load data -_EOF_ - -for i in `seq ${LINES}`; do - RSPEC="r${i}_${ROOTNAME}" - cat >> ${XSPEC_DEPROJ} << _EOF_ -data ${i}:${i} ${RSPEC}_grp.pi -response 1:${i} ${RSPEC}.wrmf -arf 1:${i} ${RSPEC}.warf -backgrnd ${i} ${RSPEC}_bkg.pi - -_EOF_ -done - -cat >> ${XSPEC_DEPROJ} << _EOF_ - -# filter needed energy range -ignore bad -ignore **:0.0-0.7,7.0-** - -method leven 1000 0.01 -# change abundance standard -abund ${ABUND} -xsect bcmc -cosmo 70 0 0.73 -xset delta 0.01 -systematic 0 -# auto answer -query yes - -# plot related -setplot energy - -# model to use -model projct*wabs*apec - 0 - 0 - 0 - ${N_H} -0.001 0 0 100000 1e+06 - 1 0.01 0.008 0.008 64 64 - 0.5 0.001 0 0 5 5 - ${REDSHIFT} -0.01 -0.999 -0.999 10 10 - 1 0.01 0 0 1e+24 1e+24 -_EOF_ - -INPUT_TIMES=`expr ${LINES} - 1` -for i in `seq ${INPUT_TIMES}`; do - cat >> ${XSPEC_DEPROJ} << _EOF_ -= 1 -= 2 -= 3 -= 4 - 1 0.01 0.008 0.008 64 64 - 0.5 0.001 0 0 5 5 -= 7 - 1 0.01 0 0 1e+24 1e+24 -_EOF_ -done -## xspec script }}} - -########################################################### -## xspec script (projected) {{{ -printf "XSPEC script for projected analysis\n" -[ -e ${XSPEC_PROJTD} ] && rm -fv ${XSPEC_PROJTD} - -cat >> ${XSPEC_PROJTD} << _EOF_ -# -# XSPEC -# radial spectra (projected analysis) -# model wabs*apec -# -# generated by script: \``basename $0`' -# `date` - -statistic chi - -# load data -_EOF_ - -for i in `seq ${LINES}`; do - RSPEC="r${i}_${ROOTNAME}" - cat >> ${XSPEC_PROJTD} << _EOF_ -data ${i}:${i} ${RSPEC}_grp.pi -response 1:${i} ${RSPEC}.wrmf -arf 1:${i} ${RSPEC}.warf -backgrnd ${i} ${RSPEC}_bkg.pi - -_EOF_ -done - -cat >> ${XSPEC_PROJTD} << _EOF_ - -# filter needed energy range -ignore bad -ignore **:0.0-0.7,7.0-** - -method leven 1000 0.01 -# change abundance standard -abund ${ABUND} -xsect bcmc -cosmo 70 0 0.73 -xset delta 0.01 -systematic 0 -# auto answer -query yes - -# plot related -setplot energy - -# model to use -model wabs*apec - ${N_H} -0.001 0 0 100000 1e+06 - 1 0.01 0.008 0.008 64 64 - 0.5 0.001 0 0 5 5 - ${REDSHIFT} -0.01 -0.999 -0.999 10 10 - 1 0.01 0 0 1e+24 1e+24 -_EOF_ - -INPUT_TIMES=`expr ${LINES} - 1` -for i in `seq ${INPUT_TIMES}`; do - cat >> ${XSPEC_PROJTD} << _EOF_ -= 1 - 1 0.01 0.008 0.008 64 64 - 0.5 0.001 0 0 5 5 -= 4 - 1 0.01 0 0 1e+24 1e+24 -_EOF_ -done - -printf "DONE\n" -## xspec script }}} -########################################################### - -printf "ALL FINISHED\n" - -# vim: set ts=8 sw=4 tw=0 fenc=utf-8 ft=sh: # 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 ## +## 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= energy= basedir= nh= z= temp= abund= [ logfile= ]\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 + diff --git a/scripts/ciao_expcorr_only.sh b/scripts/ciao_expcorr_only.sh deleted file mode 100755 index f4d390b..0000000 --- a/scripts/ciao_expcorr_only.sh +++ /dev/null @@ -1,388 +0,0 @@ -#!/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 ## -## 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= energy= basedir= nh= z= temp= abund= [ logfile= ]\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}" 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 -sv "${ROOTNAME}*band*thresh.img" "${IMG_ORIG%.fits}_thresh.fits" - # clipped exposure map - ln -sv "${ROOTNAME}*band*thresh.expmap" "${EXPMAP}" - # exposure-corrected image - ln -sv "${ROOTNAME}*band*flux.img" "${IMG_EXPCORR}" -fi - -## main }}} - -exit 0 - diff --git a/scripts/ciao_expcorr_sbp.sh b/scripts/ciao_expcorr_sbp.sh new file mode 100755 index 0000000..fc1e5ed --- /dev/null +++ b/scripts/ciao_expcorr_sbp.sh @@ -0,0 +1,236 @@ +#!/bin/sh +unalias -a +export LC_COLLATE=C + +##################################################################### +## make exposure map and exposure-corrected image +## extract surface brightness profile (revoke 'ciao_sbp.sh') +## +## ChangeLogs: +## v4.1, 2014/10/30, Weitian LI +## updated 'EXPCORR_SCRIPT' & 'EXTRACT_SBP_SCRIPT', +## removed version number in scripts filename. +## 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.sh" +EXTRACT_SBP_SCRIPT="ciao_sbp.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= sbp_reg= nh= z= temp= abund= cellreg= expcorr=\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 + diff --git a/scripts/ciao_expcorr_sbp_v4.sh b/scripts/ciao_expcorr_sbp_v4.sh deleted file mode 100755 index 861c6dd..0000000 --- a/scripts/ciao_expcorr_sbp_v4.sh +++ /dev/null @@ -1,233 +0,0 @@ -#!/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= sbp_reg= nh= z= temp= abund= cellreg= expcorr=\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 - diff --git a/scripts/ciao_genregs.sh b/scripts/ciao_genregs.sh new file mode 100755 index 0000000..643acba --- /dev/null +++ b/scripts/ciao_genregs.sh @@ -0,0 +1,277 @@ +#!/bin/sh +unalias -a +export LC_COLLATE=C + +##################################################################### +## generate 'radius spectra profile' regions +## and 'radius surface profile' regions +## +## ChangeLogs: +## v1, 2013/10/12, LIweitiaNux +## split from 'ciao_expcorr_sbp_v3.sh' +##################################################################### + +SCRIPT_PATH=`readlink -f $0` +SCRIPT_DIR=`dirname ${SCRIPT_PATH}` +XCENTROID_SCRIPT="chandra_xcentroid.sh" +GEN_SPCREG_SCRIPT="chandra_genspcreg.sh" +GEN_SBPREG_SCRIPT="chandra_gensbpreg.sh" + +## about, used in `usage' {{{ +VERSION="v1" +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= reg_in= bkgd= ds9=\n" + printf "\nversion:\n" + printf "${VERSION}, ${UPDATE}\n" + exit ${ERR_USG} + ;; +esac +## usage }}} + +## link needed files {{{ +BKGD_FILE=`ls ../bkg/bkgcorr_bl*.pi 2> /dev/null | head -n 1` +if [ -r "${BKGD_FILE}" ]; then + ln -svf ${BKGD_FILE} . +fi +ASOL_FILE=`ls ../pcad*_asol?.fits 2> /dev/null` +if [ -r "${ASOL_FILE}" ]; then + ln -svf ${ASOL_FILE} . +fi +CELL_REG_FILE=`ls ../evt/celld*.reg 2> /dev/null | grep -v 'orig'` +if [ -r "${CELL_REG_FILE}" ]; then + ln -svf ${CELL_REG_FILE} . +fi +# }}} + +## default parameters {{{ +## clean evt2 file +DFT_EVT=`ls evt2*_clean.fits 2> /dev/null` +## the repro dir +DFT_BASEDIR=".." +# default `asol file' +ASOL="`ls pcadf*_asol1.fits 2> /dev/null | head -n 1`" + +## energy range +# format: `E_START:E_END:E_WIDTH' +DFT_ENERGY=700:7000:100 +E_START=`echo ${DFT_ENERGY} | awk -F':' '{ print $1 }'` +E_END=`echo ${DFT_ENERGY} | awk -F':' '{ print $2 }'` + +## log file +DFT_LOGFILE="genreg_`date '+%Y%m%d'`.log" + +## background spectra +DFT_BKGD=`ls bkgcorr_bl*.pi | head -n 1` +## 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 "${reg_in}" ]; then + REG_IN="${reg_in}" +else + read -p "input region file: " REG_IN + if [ ! -r "${REG_IN}" ]; then + printf "ERROR: cannot access given \`${REG_IN}' evt file\n" + exit ${ERR_REG} + fi +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/\/*$//'` + +#background spectrum +if [ -r "${bkgd}" ] ;then + BKGD=${bkgd} +elif [ -r "${DFT_BKGD}" ] ; then + BKGD="${DFT_BKGD}" + #ln -svf ${DFT_BKGD} . +else + read -p ">background spectrum file: " BKGD + if [ ! -d ${BKGD} ] ; then + printf "ERROR on background spectrum file" + exit ${ERR_BKG} + fi +fi + +## ds9: flag to determine whether to use ds9 check centroid and regions +if [ ! -z "${ds9}" ]; then + case "${ds9}" in + [nN][oO]|[fF]*) + F_DS9="NO" + ;; + *) + F_DS9="YES" + ;; + esac +else + F_DS9="YES" +fi +## parameters }}} + +## 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 }}} + +## filter by energy band +printf "filter out events in energy band: \`${E_START}:${E_END}' ...\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 "======== X-RAY CENTROID =======\n" +CMD="${SCRIPT_DIR}/${XCENTROID_SCRIPT} evt=${EVT} reg=${REG_IN} conv=yes 2>&1 | tee xcentroid.dat" +printf "CMD: $CMD\n" +${SCRIPT_DIR}/${XCENTROID_SCRIPT} evt=${EVT} reg=${REG_IN} conv=yes 2>&1 | tee xcentroid.dat +X=`grep '(X,Y)' xcentroid.dat | tr -d ' XY():' | awk -F',' '{ print $2 }'` +Y=`grep '(X,Y)' xcentroid.dat | tr -d ' XY():' | awk -F',' '{ print $3 }'` +CNTRD_WCS_REG="centroid_wcs.reg" +CNTRD_PHY_REG="centroid_phy.reg" +printf "## X centroid: ($X,$Y)\n" +if [ "${F_DS9}" = "YES" ]; then + printf "check the X centroid ...\n" + ds9 ${EVT_E} -regions ${CNTRD_PHY_REG} -cmap sls -bin factor 4 +fi +X0=$X +Y0=$Y +X=`grep -i 'point' ${CNTRD_PHY_REG} | head -n 1 | tr -d 'a-zA-Z() ' | awk -F',' '{ print $1 }'` +Y=`grep -i 'point' ${CNTRD_PHY_REG} | head -n 1 | tr -d 'a-zA-Z() ' | awk -F',' '{ print $2 }'` +if [ "x${X}" != "x${X0}" ] || [ "x${Y}" != "x${Y0}" ]; then + printf "## X CENTROID CHANGED -> ($X,$Y)\n" + # update ${CNTRD_WCS_REG} + printf "update ${CNTRD_WCS_REG} ...\n" + rm -f ${CNTRD_WCS_REG} + punlearn dmcoords + dmcoords infile="${EVT}" asolfile="${ASOL}" option=sky x=${X} y=${Y} + RA=`pget dmcoords ra` + DEC=`pget dmcoords dec` + echo "point(${RA},${DEC})" > ${CNTRD_WCS_REG} +fi +printf "======== X-RAY CENTROID FINISHED =======\n\n" + +SPC_REG=rspec.reg +SBP_REG=sbprofile.reg + +printf "======== GENERATE SBPROFILE REGIONS =======\n" +CMD="${SCRIPT_DIR}/${GEN_SBPREG_SCRIPT} ${EVT} ${EVT_E} ${X} ${Y} ${BKGD} ${SBP_REG}" +printf "CMD: ${CMD}\n" +${SCRIPT_DIR}/${GEN_SBPREG_SCRIPT} ${EVT} ${EVT_E} ${X} ${Y} ${BKGD} ${SBP_REG} +if [ "${F_DS9}" = "YES" ]; then + printf "check SBP regions ...\n" + ds9 ${EVT_E} -regions ${SBP_REG} -cmap sls -bin factor 4 +fi +printf "======== GENERATE SBPROFILE REGIONS FINISHED =======\n\n" + +printf "======== GENERATE SPECTRUM REGIONS =======\n" +CMD="${SCRIPT_DIR}/${GEN_SPCREG_SCRIPT} ${EVT} ${EVT_E} ${BKGD} ${X} ${Y} ${SPC_REG}" +printf "CMD: ${CMD}\n" +${SCRIPT_DIR}/${GEN_SPCREG_SCRIPT} ${EVT} ${EVT_E} ${BKGD} ${X} ${Y} ${SPC_REG} +if [ "${F_DS9}" = "YES" ]; then + printf "check SPC regions ...\n" + ds9 ${EVT_E} -regions ${SPC_REG} -cmap sls -bin factor 4 +fi +printf "======== GENERATE SPECTRUM REGIONS FINISHED =======\n\n" + +exit 0 + diff --git a/scripts/ciao_genregs_v1.sh b/scripts/ciao_genregs_v1.sh deleted file mode 100755 index 643acba..0000000 --- a/scripts/ciao_genregs_v1.sh +++ /dev/null @@ -1,277 +0,0 @@ -#!/bin/sh -unalias -a -export LC_COLLATE=C - -##################################################################### -## generate 'radius spectra profile' regions -## and 'radius surface profile' regions -## -## ChangeLogs: -## v1, 2013/10/12, LIweitiaNux -## split from 'ciao_expcorr_sbp_v3.sh' -##################################################################### - -SCRIPT_PATH=`readlink -f $0` -SCRIPT_DIR=`dirname ${SCRIPT_PATH}` -XCENTROID_SCRIPT="chandra_xcentroid.sh" -GEN_SPCREG_SCRIPT="chandra_genspcreg.sh" -GEN_SBPREG_SCRIPT="chandra_gensbpreg.sh" - -## about, used in `usage' {{{ -VERSION="v1" -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= reg_in= bkgd= ds9=\n" - printf "\nversion:\n" - printf "${VERSION}, ${UPDATE}\n" - exit ${ERR_USG} - ;; -esac -## usage }}} - -## link needed files {{{ -BKGD_FILE=`ls ../bkg/bkgcorr_bl*.pi 2> /dev/null | head -n 1` -if [ -r "${BKGD_FILE}" ]; then - ln -svf ${BKGD_FILE} . -fi -ASOL_FILE=`ls ../pcad*_asol?.fits 2> /dev/null` -if [ -r "${ASOL_FILE}" ]; then - ln -svf ${ASOL_FILE} . -fi -CELL_REG_FILE=`ls ../evt/celld*.reg 2> /dev/null | grep -v 'orig'` -if [ -r "${CELL_REG_FILE}" ]; then - ln -svf ${CELL_REG_FILE} . -fi -# }}} - -## default parameters {{{ -## clean evt2 file -DFT_EVT=`ls evt2*_clean.fits 2> /dev/null` -## the repro dir -DFT_BASEDIR=".." -# default `asol file' -ASOL="`ls pcadf*_asol1.fits 2> /dev/null | head -n 1`" - -## energy range -# format: `E_START:E_END:E_WIDTH' -DFT_ENERGY=700:7000:100 -E_START=`echo ${DFT_ENERGY} | awk -F':' '{ print $1 }'` -E_END=`echo ${DFT_ENERGY} | awk -F':' '{ print $2 }'` - -## log file -DFT_LOGFILE="genreg_`date '+%Y%m%d'`.log" - -## background spectra -DFT_BKGD=`ls bkgcorr_bl*.pi | head -n 1` -## 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 "${reg_in}" ]; then - REG_IN="${reg_in}" -else - read -p "input region file: " REG_IN - if [ ! -r "${REG_IN}" ]; then - printf "ERROR: cannot access given \`${REG_IN}' evt file\n" - exit ${ERR_REG} - fi -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/\/*$//'` - -#background spectrum -if [ -r "${bkgd}" ] ;then - BKGD=${bkgd} -elif [ -r "${DFT_BKGD}" ] ; then - BKGD="${DFT_BKGD}" - #ln -svf ${DFT_BKGD} . -else - read -p ">background spectrum file: " BKGD - if [ ! -d ${BKGD} ] ; then - printf "ERROR on background spectrum file" - exit ${ERR_BKG} - fi -fi - -## ds9: flag to determine whether to use ds9 check centroid and regions -if [ ! -z "${ds9}" ]; then - case "${ds9}" in - [nN][oO]|[fF]*) - F_DS9="NO" - ;; - *) - F_DS9="YES" - ;; - esac -else - F_DS9="YES" -fi -## parameters }}} - -## 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 }}} - -## filter by energy band -printf "filter out events in energy band: \`${E_START}:${E_END}' ...\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 "======== X-RAY CENTROID =======\n" -CMD="${SCRIPT_DIR}/${XCENTROID_SCRIPT} evt=${EVT} reg=${REG_IN} conv=yes 2>&1 | tee xcentroid.dat" -printf "CMD: $CMD\n" -${SCRIPT_DIR}/${XCENTROID_SCRIPT} evt=${EVT} reg=${REG_IN} conv=yes 2>&1 | tee xcentroid.dat -X=`grep '(X,Y)' xcentroid.dat | tr -d ' XY():' | awk -F',' '{ print $2 }'` -Y=`grep '(X,Y)' xcentroid.dat | tr -d ' XY():' | awk -F',' '{ print $3 }'` -CNTRD_WCS_REG="centroid_wcs.reg" -CNTRD_PHY_REG="centroid_phy.reg" -printf "## X centroid: ($X,$Y)\n" -if [ "${F_DS9}" = "YES" ]; then - printf "check the X centroid ...\n" - ds9 ${EVT_E} -regions ${CNTRD_PHY_REG} -cmap sls -bin factor 4 -fi -X0=$X -Y0=$Y -X=`grep -i 'point' ${CNTRD_PHY_REG} | head -n 1 | tr -d 'a-zA-Z() ' | awk -F',' '{ print $1 }'` -Y=`grep -i 'point' ${CNTRD_PHY_REG} | head -n 1 | tr -d 'a-zA-Z() ' | awk -F',' '{ print $2 }'` -if [ "x${X}" != "x${X0}" ] || [ "x${Y}" != "x${Y0}" ]; then - printf "## X CENTROID CHANGED -> ($X,$Y)\n" - # update ${CNTRD_WCS_REG} - printf "update ${CNTRD_WCS_REG} ...\n" - rm -f ${CNTRD_WCS_REG} - punlearn dmcoords - dmcoords infile="${EVT}" asolfile="${ASOL}" option=sky x=${X} y=${Y} - RA=`pget dmcoords ra` - DEC=`pget dmcoords dec` - echo "point(${RA},${DEC})" > ${CNTRD_WCS_REG} -fi -printf "======== X-RAY CENTROID FINISHED =======\n\n" - -SPC_REG=rspec.reg -SBP_REG=sbprofile.reg - -printf "======== GENERATE SBPROFILE REGIONS =======\n" -CMD="${SCRIPT_DIR}/${GEN_SBPREG_SCRIPT} ${EVT} ${EVT_E} ${X} ${Y} ${BKGD} ${SBP_REG}" -printf "CMD: ${CMD}\n" -${SCRIPT_DIR}/${GEN_SBPREG_SCRIPT} ${EVT} ${EVT_E} ${X} ${Y} ${BKGD} ${SBP_REG} -if [ "${F_DS9}" = "YES" ]; then - printf "check SBP regions ...\n" - ds9 ${EVT_E} -regions ${SBP_REG} -cmap sls -bin factor 4 -fi -printf "======== GENERATE SBPROFILE REGIONS FINISHED =======\n\n" - -printf "======== GENERATE SPECTRUM REGIONS =======\n" -CMD="${SCRIPT_DIR}/${GEN_SPCREG_SCRIPT} ${EVT} ${EVT_E} ${BKGD} ${X} ${Y} ${SPC_REG}" -printf "CMD: ${CMD}\n" -${SCRIPT_DIR}/${GEN_SPCREG_SCRIPT} ${EVT} ${EVT_E} ${BKGD} ${X} ${Y} ${SPC_REG} -if [ "${F_DS9}" = "YES" ]; then - printf "check SPC regions ...\n" - ds9 ${EVT_E} -regions ${SPC_REG} -cmap sls -bin factor 4 -fi -printf "======== GENERATE SPECTRUM REGIONS FINISHED =======\n\n" - -exit 0 - diff --git a/scripts/ciao_procevt.sh b/scripts/ciao_procevt.sh new file mode 100755 index 0000000..d73bfaf --- /dev/null +++ b/scripts/ciao_procevt.sh @@ -0,0 +1,185 @@ +#!/bin/sh +# +########################################################### +## process `evt2 file' generated by `chandra_repro' ## +## to produce a `clean evt2 file' ## +## (remove point source and flares) ## +## for later image and spectral analysis ## +## ## +## NOTES: ## +## based on previous `ciao_process_evt2.sh' ## +## support ACIS-I(chip: 0-3) and ACIS-S(chip: 7) ## +## determine by check `DETNAM' for chip number ## +## if `DETNAM' has `0123', then `ACIS-I' ## +## if `DETNAM' has `7', then `ACIS-S' ## +## ## +## LIweitiaNux ## +## August 16, 2012 ## +########################################################### + +########################################################### +## ChangeLogs: +## v2.2, 2014/10/30, Weitian LI +## small fix to the generation of chips script: +## changed '>>' to '>' +## v2.1, 2012/08/16, LIweitiaNux +## improve invoke `chips', run it in a separate terminal +########################################################### + +## about, used in `usage' {{{ +VERSION="v2" +UPDATE="2012-08-16" +## 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 +## error code }}} + +## usage, help {{{ +case "$1" in + -[hH]*|--[hH]*) + printf "usage:\n" + printf " `basename $0` evt=\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 acisf?????*_repro_evt2.fits 2> /dev/null`" +## 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 }}} + +## 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 "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" +## parameters }}} + +## determine ACIS type {{{ +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" + ROOTNAME="evt2_c`echo ${CCD} | tr ':' '-'`" +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" + ROOTNAME="evt2_c${CCD}" +else + printf "ERROR: unknown detector type: ${DETNAM}\n" + exit ${ERR_DET} +fi +## ACIS type }}} + +## main process {{{ +printf "filter raw evt2 file by ccd_id ...\n" +EVT2_ORIG="${ROOTNAME}_orig.fits" +punlearn dmcopy +dmcopy infile="${EVT}[ccd_id=${CCD}]" outfile=${EVT2_ORIG} clobber=yes + +# detect point sources +printf "using \`celldetect' to find the point sources ...\n" +CELLD="celld_${ROOTNAME}" +[ -e "${CELLD}.reg" ] && mv -fv ${CELLD}.reg ${CELLD}.reg_bak +punlearn celldetect +celldetect infile=${EVT2_ORIG} outfile="${CELLD}.fits" \ + regfile="${CELLD}.reg" clobber=yes + +printf "check the result of \`celldetect' ...\n" +printf "modify if necessary and save as the same name, \`${CELLD}.reg'\n" +cp -fv ${CELLD}.reg ${CELLD}_orig.reg +ds9 ${EVT2_ORIG} -region ${CELLD}.reg + +EVT2_RMSRCS="${ROOTNAME}_rmsrcs.fits" +punlearn dmcopy +dmcopy infile="${EVT2_ORIG}[exclude sky=region(${CELLD}.reg)]" \ + outfile=${EVT2_RMSRCS} clobber=yes + +LC_REG="ex_bkg.reg" +printf "filter flares ...\n" +printf "select a big source region and save as \`${LC_REG}'\n" +touch ${LC_REG} +ds9 ${EVT2_RMSRCS} + +printf "create the lightcurve ...\n" +LC="${LC_REG%.reg}.lc" +[ -e "${LC}" ] && mv -fv ${LC} ${LC}_bak +punlearn dmextract +dmextract infile="${EVT2_RMSRCS}[exclude sky=region(${LC_REG})][bin time=::200]" \ + outfile="${LC}" opt=ltc1 clobber=yes + +# generate a script for `chips' +# for convenience and backup +LC_SCRIPT="_${LC%.lc}.chips" +LC_SCALE=1.2 +GTI="${LC%.lc}.gti" +cat > ${LC_SCRIPT} << _EOF_ +from lightcurves import * +lc_clean("${LC}") +lc_clean("${LC}", scale=${LC_SCALE}, outfile="${GTI}") +print_window("${LC%.lc}_lc.jpg", ["format", "jpg", "clobber", "True"]) +_EOF_ + +printf "generate GTI ...\n" +# [-x] let chips run in a separate terminal +chips -x ${LC_SCRIPT} + +printf "remove the background flare ...\n" +EVT2_CLEAN="${ROOTNAME}_clean.fits" +punlearn dmcopy +dmcopy infile="${EVT2_RMSRCS}[@${GTI}]" outfile=${EVT2_CLEAN} clobber=yes +## main process }}} + +printf "FINISHED\n" + + diff --git a/scripts/ciao_procevt_v2.sh b/scripts/ciao_procevt_v2.sh deleted file mode 100755 index 6cfe45d..0000000 --- a/scripts/ciao_procevt_v2.sh +++ /dev/null @@ -1,182 +0,0 @@ -#!/bin/sh -# -########################################################### -## process `evt2 file' generated by `chandra_repro' ## -## to produce a `clean evt2 file' ## -## (remove point source and flares) ## -## for later image and spectral analysis ## -## ## -## NOTES: ## -## based on previous `ciao_process_evt2.sh' ## -## support ACIS-I(chip: 0-3) and ACIS-S(chip: 7) ## -## determine by check `DETNAM' for chip number ## -## if `DETNAM' has `0123', then `ACIS-I' ## -## if `DETNAM' has `7', then `ACIS-S' ## -## ## -## LIweitiaNux ## -## August 16, 2012 ## -########################################################### - -########################################################### -## ChangeLogs: -## v2.1, 2012/08/16, LIweitiaNux -## improve invoke `chips', run it in a separate terminal -########################################################### - -## about, used in `usage' {{{ -VERSION="v2" -UPDATE="2012-08-16" -## 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 -## error code }}} - -## usage, help {{{ -case "$1" in - -[hH]*|--[hH]*) - printf "usage:\n" - printf " `basename $0` evt=\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 acisf?????*_repro_evt2.fits 2> /dev/null`" -## 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 }}} - -## 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 "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" -## parameters }}} - -## determine ACIS type {{{ -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" - ROOTNAME="evt2_c`echo ${CCD} | tr ':' '-'`" -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" - ROOTNAME="evt2_c${CCD}" -else - printf "ERROR: unknown detector type: ${DETNAM}\n" - exit ${ERR_DET} -fi -## ACIS type }}} - -## main process {{{ -printf "filter raw evt2 file by ccd_id ...\n" -EVT2_ORIG="${ROOTNAME}_orig.fits" -punlearn dmcopy -dmcopy infile="${EVT}[ccd_id=${CCD}]" outfile=${EVT2_ORIG} clobber=yes - -# detect point sources -printf "using \`celldetect' to find the point sources ...\n" -CELLD="celld_${ROOTNAME}" -[ -e "${CELLD}.reg" ] && mv -fv ${CELLD}.reg ${CELLD}.reg_bak -punlearn celldetect -celldetect infile=${EVT2_ORIG} outfile="${CELLD}.fits" \ - regfile="${CELLD}.reg" clobber=yes - -printf "check the result of \`celldetect' ...\n" -printf "modify if necessary and save as the same name, \`${CELLD}.reg'\n" -cp -fv ${CELLD}.reg ${CELLD}_orig.reg -ds9 ${EVT2_ORIG} -region ${CELLD}.reg - -EVT2_RMSRCS="${ROOTNAME}_rmsrcs.fits" -punlearn dmcopy -dmcopy infile="${EVT2_ORIG}[exclude sky=region(${CELLD}.reg)]" \ - outfile=${EVT2_RMSRCS} clobber=yes - -LC_REG="ex_bkg.reg" -printf "filter flares ...\n" -printf "select a big source region and save as \`${LC_REG}'\n" -touch ${LC_REG} -ds9 ${EVT2_RMSRCS} - -printf "create the lightcurve ...\n" -LC="${LC_REG%.reg}.lc" -[ -e "${LC}" ] && mv -fv ${LC} ${LC}_bak -punlearn dmextract -dmextract infile="${EVT2_RMSRCS}[exclude sky=region(${LC_REG})][bin time=::200]" \ - outfile="${LC}" opt=ltc1 clobber=yes - -# generate a script for `chips' -# for convenience and backup -LC_SCRIPT="_${LC%.lc}.chips" -LC_SCALE=1.2 -GTI="${LC%.lc}.gti" -cat >> ${LC_SCRIPT} << _EOF_ -from lightcurves import * -lc_clean("${LC}") -lc_clean("${LC}", scale=${LC_SCALE}, outfile="${GTI}") -print_window("${LC%.lc}_lc.jpg", ["format", "jpg", "clobber", "True"]) -_EOF_ - -printf "generate GTI ...\n" -# [-x] let chips run in a separate terminal -chips -x ${LC_SCRIPT} - -printf "remove the background flare ...\n" -EVT2_CLEAN="${ROOTNAME}_clean.fits" -punlearn dmcopy -dmcopy infile="${EVT2_RMSRCS}[@${GTI}]" outfile=${EVT2_CLEAN} clobber=yes -## main process }}} - -printf "FINISHED\n" - - diff --git a/scripts/ciao_r500avgt.sh b/scripts/ciao_r500avgt.sh new file mode 100755 index 0000000..9a8dffd --- /dev/null +++ b/scripts/ciao_r500avgt.sh @@ -0,0 +1,539 @@ +#!/bin/sh +# +unalias -a +export LC_COLLATE=C +########################################################### +## script to extract spectrum and prepare needed files ## +## for calculating the `average temperature' ## +## within (0.1-0.5 r500) region ## +## ## +## NOTE: ## +## 1) r500 default in unit `pixel', if in `kpc', ## +## then `redshift z' and `calculator' are needed. ## +## 2) `background process' is same as `deproj_spectra' ## +## which supports `spectrum', `local' and `blanksky' ## +## 3) ARF/RMF files either provided or use the ARF/RMF ## +## of the outmost region ## +## ## +## LIweitiaNux ## +## August 22, 2012 ## +########################################################### + +########################################################### +## ChangeLogs +## v1.1, 2012/08/26, LIweitiaNux +## modify `KPC_PER_PIX', able to use the newest version `calc_distance' +## v1.2, 2012/08/26, LIweitiaNux +## fix a bug with `DFT_BKGD' +## v2.0, 2012/09/04, LIweitiaNux +## add parameter `inner' and `outer' to adjust the region range +## modify parameter `r500' to take `kpc' as the default unit +## v2.1, 2012/10/05, LIweitiaNux +## change `DFT_GRP_CMD' to `group 1 128 4 ...' +## v3.0, 2013/02/09, LIweitiaNux +## modify for new process +########################################################### + +## about, used in `usage' {{{ +VERSION="v3.0" +UPDATE="2013-02-09" +## about }}} + +## error code {{{ +ERR_USG=1 +ERR_DIR=11 +ERR_EVT=12 +ERR_BKG=13 +ERR_REG=14 +ERR_JSON=15 +ERR_ASOL=21 +ERR_BPIX=22 +ERR_PBK=23 +ERR_MSK=24 +ERR_BKGTY=31 +ERR_SPEC=32 +ERR_ARF=51 +ERR_RMF=52 +ERR_UNI=61 +## error code }}} + +## usage, help {{{ +case "$1" in + -[hH]*|--[hH]*) + printf "usage:\n" + printf " `basename $0` evt= r500= basedir= info= inner= outer= regin= regout= bkgd= nh= z= arf= rmf= [ grpcmd= log= ]\n" + printf "\nversion:\n" + printf "${VERSION}, ${UPDATE}\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 +COSCALC="`which cosmo_calc calc_distance 2>/dev/null | head -n 1`" +# COSCALC="_path_to_calc_distance_" +# COSCALC="$HOME/bin/mass/calc_distance" +if [ -z "${COSCALC}" ] || [ ! -x ${COSCALC} ]; then + printf "ERROR: \`COSCALC: ${COSCALC}' neither specified nor executable\n" + exit 255 +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`" +# default `bkgd', use `bkgcorr_blanksky*' corrected bkg spectrum +DFT_BKGD="`ls bkgcorr_blanksky_*.pi 2> /dev/null`" +# default basedir +DFT_BASEDIR="../.." +# default `radial region file' +#DFT_REG_IN="_NOT_EXIST_" +DFT_REG_IN="rspec.reg" +# default region range (0.1-0.5 R500) +DFT_INNER="0.1" +DFT_OUTER="0.5" +# default ARF/RMF, the one of the outmost region +DFT_ARF="`ls -1 r?_*.warf 2> /dev/null | tail -n 1`" +DFT_RMF="`ls -1 r?_*.wrmf 2> /dev/null | tail -n 1`" + +# default `group command' for `grppha' +DFT_GRP_CMD="group min 25" +#DFT_GRP_CMD="group 1 128 2 129 256 4 257 512 8 513 1024 16" +#DFT_GRP_CMD="group 1 128 4 129 256 8 257 512 16 513 1024 32" +# default `log file' +DFT_LOGFILE="r500avgt_`date '+%Y%m%d%H'`.log" + +# default JSON pattern +DFT_JSON_PAT="*_INFO.json" +## 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 +} + +## background renormalization (BACKSCAL) {{{ +# renorm background according to particle background +# energy range: 9.5-12.0 keV (channel: 651-822) +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 }'` + punlearn dmkeypar + EXPTIME=`dmkeypar $1 EXPOSURE echo=yes` + BACK=`dmkeypar $1 BACKSCAL echo=yes` + # fix `scientific notation' bug for `bc' + EXPTIME_B=`echo ${EXPTIME} | sed 's/[eE]/\*10\^/' | sed 's/+//'` + BACK_B=`echo "( ${BACK} )" | sed 's/[eE]/\*10\^/' | sed 's/+//'` + PB_FLUX=`echo "scale = 16; ${COUNTS} / ${EXPTIME_B} / ${BACK_B}" | bc -l` + echo ${PB_FLUX} +} + +bkg_renorm() { + # $1: src spectrum, $2: back spectrum + PBFLUX_SRC=`pb_flux $1` + PBFLUX_BKG=`pb_flux $2` + BACK_OLD=`dmkeypar $2 BACKSCAL echo=yes` + BACK_OLD_B=`echo "( ${BACK_OLD} )" | sed 's/[eE]/\*10\^/' | sed 's/+//'` + BACK_NEW=`echo "scale = 16; ${BACK_OLD_B} * ${PBFLUX_BKG} / ${PBFLUX_SRC}" | bc -l` + printf "\`$2': BACKSCAL:\n" + printf " ${BACK_OLD} --> ${BACK_NEW}\n" + punlearn dmhedit + dmhedit infile=$2 filelist=none operation=add \ + key=BACKSCAL value=${BACK_NEW} comment="old value: ${BACK_OLD}" +} +## bkg renorm }}} +## functions end }}} + +## 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 +if [ -d "${basedir}" ]; then + BASEDIR=${basedir} +else + BASEDIR=${DFT_BASEDIR} +fi +if [ ! -z "${json}" ] && [ -r "${BASEDIR}/${json}" ]; then + JSON_FILE="${BASEDIR}/${json}" +elif ls ${BASEDIR}/${DFT_JSON_PAT} > /dev/null 2>&1; then + JSON_FILE=`ls ${BASEDIR}/${DFT_JSON_PAT}` +else + read -p "> JSON_file: " JSON_FILE + if [ ! -r "${JSON_FILE}" ]; then + printf "ERROR: cannot access given \`${JSON_FILE}'\n" + exit ${ERR_JSON} + fi +fi +printf "## use json_file: \`${JSON_FILE}'\n" | ${TOLOG} + +# process `nh' and `redshift' {{{ +NH_JSON=`grep '"nH' ${JSON_FILE} | sed 's/.*"nH.*":\ //' | sed 's/\ *,$//'` +Z_JSON=`grep '"redshift' ${JSON_FILE} | sed 's/.*"redshift.*":\ //' | sed 's/\ *,$//'` +printf "## get nh: \`${NH_JSON}' (from \`${JSON_FILE}')\n" | ${TOLOG} +printf "## get redshift: \`${Z_JSON}' (from \`${JSON_FILE}')\n" | ${TOLOG} +## if `nh' and `redshift' supplied in cmdline, then use them +if [ ! -z "${nh}" ]; then + N_H=${nh} +else + N_H=${NH_JSON} +fi +# redshift +if [ ! -z "${z}" ]; then + REDSHIFT=${z} +else + REDSHIFT=${Z_JSON} +fi +printf "## use nH: ${N_H}\n" | ${TOLOG} +printf "## use redshift: ${REDSHIFT}\n" | ${TOLOG} +# nh & redshift }}} + +# region range {{{ +if [ ! -z "${inner}" ]; then + INNER="${inner}" +else + INNER=${DFT_INNER} +fi +if [ ! -z "${outer}" ]; then + OUTER="${outer}" +else + OUTER=${DFT_OUTER} +fi +printf "## region range: (${INNER} - ${OUTER} R500)\n" | ${TOLOG} +# range }}} + +# process `r500' {{{ +R500_RAW=`grep '"R500.*kpc' ${JSON_FILE} | sed 's/.*"R500.*":\ //' | sed 's/\ *,$//'` +if [ ! -z "${r500}" ]; then + R500_RAW=${r500} +fi +if [ -z "${R500_RAW}" ]; then + printf "## input R500 followed with unit, e.g.: 800kpc, 400pix\n" + read -p "> value of \`R500' (in pixel/kpc): " R500_RAW +fi +R500_VAL=`echo "${R500_RAW}" | tr -d 'a-zA-Z, '` +R500_UNI=`echo "${R500_RAW}" | tr -d '0-9, '` +printf "## get \`R500': ${R500_VAL} in unit \`${R500_UNI}'\n" | ${TOLOG} + +# if in kpc, convert to pix +case "${R500_UNI}" in + [pP]*) + printf "## units in \`pixel', conversion not needed\n" | ${TOLOG} + R500_PIX_B=`echo ${R500_VAL} | sed 's/[eE]/\*10\^/' | sed 's/+//'` + ;; + *) + printf "## units in \`kpc', convert to \`Chandra pixel'\n" | ${TOLOG} + 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" + R500_VAL_B=`echo ${R500_VAL} | sed 's/[eE]/\*10\^/' | sed 's/+//'` + R500_PIX_B=`echo "scale = 4; ${R500_VAL_B} / ( ${KPC_PER_PIX_B} )" | bc -l` + ;; +esac +# calc (inner-outer R500) +R_IN=`echo "scale = 4; ${INNER} * ${R500_PIX_B}" | bc -l` +R_OUT=`echo "scale = 4; ${OUTER} * ${R500_PIX_B}" | bc -l` +printf "## R500 in units pixel: ${R500_PIX_B}\n" | ${TOLOG} +printf "## (${INNER}-${OUTER} R500) range in pixel: ${R_IN} - ${R_OUT}\n" | ${TOLOG} +# r500 }}} + +# 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} + +# input and output region files {{{ +if [ -r "${regin}" ]; then + REG_IN="${regin}" +elif [ -r "${DFT_REG_IN}" ]; then + REG_IN=${DFT_REG_IN} +else + read -p "> previous used radial spec regfile: " REG_IN + if [ ! -r "${REG_IN}" ]; then + printf "ERROR: cannot access given \`${REG_IN}' region file\n" + exit ${ERR_REG} + fi +fi +printf "## use previous regfile: \`${REG_IN}'\n" | ${TOLOG} +if [ ! -z "${regout}" ]; then + REG_OUT="${regout}" +else + REG_OUT="r500avgt_${INNER}-${OUTER}.reg" +fi +[ -e "${REG_OUT}" ] && mv -fv ${REG_OUT} ${REG_OUT}_bak +printf "## set output regfile: \`${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` +XC=`echo ${TMP_REG} | tr -d 'a-zA-Z() ' | awk -F',' '{ print $1 }'` +YC=`echo ${TMP_REG} | tr -d 'a-zA-Z() ' | awk -F',' '{ print $2 }'` +printf "## get center coord: (${XC},${YC})\n" | ${TOLOG} +# region files }}} + +# check given bkgd, determine background {{{ +if [ -r "${bkgd}" ]; then + BKGD=${bkgd} +elif [ -r "${DFT_BKGD}" ]; then + BKGD=${DFT_BKGD} +else + read -p "> background (blanksky_evt | lbkg_reg | bkg_spec): " BKGD + if [ ! -r "${BKGD}" ]; then + printf "ERROR: cannot access given \`${BKGD}'\n" + exit ${ERR_BKG} + fi +fi +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 + printf "## given \`${BKGD}' is a \`text file'\n" + printf "## use it as local bkg region file\n" + printf "## use *LOCAL BKG SPEC*\n" | ${TOLOG} + # just set flags, extract spectrum later + USE_LBKG_REG=YES + USE_BLANKSKY=NO + USE_BKG_SPEC=NO +elif file -bL ${BKGD} | grep -qi 'FITS'; then + printf "## given \`${BKGD}' is a \`FITS file'\n" + # get FITS header keyword + HDUCLAS1=`dmkeypar ${BKGD} HDUCLAS1 echo=yes` + if [ "${HDUCLAS1}" = "EVENTS" ]; then + # event file + printf "## given file is \`event'\n" + # check if `blanksky' or `stowed bkg' + BKG_OBJ=`dmkeypar ${BKGD} OBJECT echo=yes` + if [ "${BKG_OBJ}" = "BACKGROUND DATASET" ] || [ "${BKG_OBJ}" = "ACIS STOWED" ]; then + # valid bkg evt file + printf "## given FITS file is a valid bkgrnd file\n" + printf "## use *BLANKSKY*\n" | ${TOLOG} + USE_BLANKSKY=YES + USE_LBKG_REG=NO + USE_BKG_SPEC=NO + # specify `BLANKSKY' + BLANKSKY=${BKGD} + else + # invalid bkg evt file + printf "ERROR: invalid bkg evt file given\n" + exit ${ERR_BKGTY} + fi + elif [ "${HDUCLAS1}" = "SPECTRUM" ]; then + # spectrum file + printf "## given file is \`spectrum'\n" + printf "## use *BKG SPECTRUM*\n" | ${TOLOG} + USE_BKG_SPEC=YES + USE_BLANKSKY=NO + USE_LBKG_REG=NO + # specify `BKG_SPEC' + BKG_SPEC=${BKGD} + else + # other type + printf "ERROR: other type FITS given\n" + exit ${ERR_BKGTY} + fi +else + printf "ERROR: given \`${BKGD}' type UNKNOWN\n" + exit ${ERR_BKGTY} +fi +# bkgd }}} + +# check `arf' and `rmf' {{{ +if [ -r "${arf}" ]; then + ARF=${arf} +elif [ -r "${DFT_ARF}" ]; then + ARF=${DFT_ARF} +else + read -p "> provide the ARF to use: " ARF + if [ ! -r "${ARF}" ]; then + printf "ERROR: cannot access given \`${ARF}'\n" + exit ${ERR_ARF} + fi +fi +printf "## use ARF: \`${ARF}'\n" | ${TOLOG} +# rmf +if [ -r "${rmf}" ]; then + RMF=${rmf} +elif [ -r "${DFT_RMF}" ]; then + RMF=${DFT_RMF} +else + read -p "> provide the RMF to use: " RMF + if [ ! -r "${RMF}" ]; then + printf "ERROR: cannot access given \`${RMF}'\n" + exit ${ERR_RMF} + fi +fi +printf "## use RMF: \`${RMF}'\n" | ${TOLOG} +# arf & rmf }}} + +# check given `grpcmd' +if [ ! -z "${grpcmd}" ]; then + GRP_CMD="${grpcmd}" +else + GRP_CMD="${DFT_GRP_CMD}" +fi +printf "## use grppha cmd: \`${GRP_CMD}'\n" | ${TOLOG} +## parameters }}} + +################################################## +#### main +## region related {{{ +## generate the needed region file +printf "generate the output region file ...\n" +cat > ${REG_OUT} << _EOF_ +# Region file format: CIAO version 1.0 +pie(${XC},${YC},${R_IN},${R_OUT},0,360) +_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} -region ${REG_OUT} + +## 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'` +if [ "x${INVALID}" != "x" ]; then + printf "*** WARNING: there are pie regions' END_ANGLE > 360\n" | ${TOLOG} + printf "*** will to fix ...\n" + mv -fv ${REG_OUT} ${REG_OUT}_tmp + # using `awk' to fix + awk -F'[,()]' '{ + if ($7 > 360) { + printf "%s(%.2f,%.2f,%.2f,%.2f,%.2f,%.2f)\n", $1,$2,$3,$4,$5,$6,($7-360) + } + else { + print $0 + } + }' ${REG_OUT}_tmp > ${REG_OUT} + rm -f ${REG_OUT}_tmp +fi +## region related }}} + +## generate spectrum {{{ +# object +AVGT_SPEC="${REG_OUT%.reg}.pi" +AVGT_SPEC_GRP="${AVGT_SPEC%.pi}_grp.pi" +printf "extract object spectrum \`${AVGT_SPEC}' ...\n" +punlearn dmextract +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 + +# background +printf "generate the background spectrum ...\n" +AVGT_BKG="${AVGT_SPEC%.pi}_bkg.pi" +if [ "${USE_BLANKSKY}" = "YES" ]; then + # use blanksky as background file + printf "extract spectrum from blanksky ...\n" + punlearn dmextract + dmextract infile="${BLANKSKY}[sky=region(${REG_OUT})][bin PI]" \ + outfile=${AVGT_BKG} wmap="[bin det=8]" clobber=yes +elif [ "${USE_LBKG_REG}" = "YES" ]; then + printf "extract local background ...\n" + punlearn dmextract + dmextract infile="${EVT}[sky=region(${BKGD})][bin PI]" \ + outfile=${AVGT_BKG} wmap="[bin det=8]" clobber=yes +elif [ "${USE_BKG_SPEC}" = "YES" ]; then + printf "copy specified background spectrum ...\n" + cp -fv ${BKG_SPEC} ${AVGT_BKG} +fi + +printf "renormalize the background ...\n" +bkg_renorm ${AVGT_SPEC} ${AVGT_BKG} + +## spectrum }}} + +## generate XSPEC script {{{ +printf "generate a XSPEC script ...\n" +# default output xspec scripts +XSPEC_SCRIPT="xspec_${REG_OUT%.reg}.xcm" +[ -e "${XSPEC_SCRIPT}" ] && mv -fv ${XSPEC_SCRIPT} ${XSPEC_SCRIPT}_bak +cat > ${XSPEC_SCRIPT} << _EOF_ +## XSPEC script +## spectrum analysis to get the average temperatue with (0.1-0.5 R500) +## +## generated by: \``basename $0`' +## date: \``date`' +## + +# xspec settings +statistic chi +abund grsa +query yes + +# data +data ${AVGT_SPEC_GRP} +response ${RMF} +arf ${ARF} +backgrnd ${AVGT_BKG} + +# fitting range +ignore bad +ignore 0.0-0.7,7.0-** + +# plot related +setplot energy + +method leven 10 0.01 +xsect bcmc +cosmo 70 0 0.73 +xset delta 0.01 +systematic 0 + +# model +model wabs*apec + ${N_H} -0.001 0 0 100000 1e+06 + 1.0 0.01 0.008 0.008 64 64 + 0.4 0.001 0 0 5 5 + ${REDSHIFT} -0.01 -0.999 -0.999 10 10 + 0.0 0.01 0 0 1e+24 1e+24 + +## xspec script end +_EOF_ +## xspec script }}} + +exit 0 + diff --git a/scripts/ciao_r500avgt_v3.sh b/scripts/ciao_r500avgt_v3.sh deleted file mode 100755 index 9a8dffd..0000000 --- a/scripts/ciao_r500avgt_v3.sh +++ /dev/null @@ -1,539 +0,0 @@ -#!/bin/sh -# -unalias -a -export LC_COLLATE=C -########################################################### -## script to extract spectrum and prepare needed files ## -## for calculating the `average temperature' ## -## within (0.1-0.5 r500) region ## -## ## -## NOTE: ## -## 1) r500 default in unit `pixel', if in `kpc', ## -## then `redshift z' and `calculator' are needed. ## -## 2) `background process' is same as `deproj_spectra' ## -## which supports `spectrum', `local' and `blanksky' ## -## 3) ARF/RMF files either provided or use the ARF/RMF ## -## of the outmost region ## -## ## -## LIweitiaNux ## -## August 22, 2012 ## -########################################################### - -########################################################### -## ChangeLogs -## v1.1, 2012/08/26, LIweitiaNux -## modify `KPC_PER_PIX', able to use the newest version `calc_distance' -## v1.2, 2012/08/26, LIweitiaNux -## fix a bug with `DFT_BKGD' -## v2.0, 2012/09/04, LIweitiaNux -## add parameter `inner' and `outer' to adjust the region range -## modify parameter `r500' to take `kpc' as the default unit -## v2.1, 2012/10/05, LIweitiaNux -## change `DFT_GRP_CMD' to `group 1 128 4 ...' -## v3.0, 2013/02/09, LIweitiaNux -## modify for new process -########################################################### - -## about, used in `usage' {{{ -VERSION="v3.0" -UPDATE="2013-02-09" -## about }}} - -## error code {{{ -ERR_USG=1 -ERR_DIR=11 -ERR_EVT=12 -ERR_BKG=13 -ERR_REG=14 -ERR_JSON=15 -ERR_ASOL=21 -ERR_BPIX=22 -ERR_PBK=23 -ERR_MSK=24 -ERR_BKGTY=31 -ERR_SPEC=32 -ERR_ARF=51 -ERR_RMF=52 -ERR_UNI=61 -## error code }}} - -## usage, help {{{ -case "$1" in - -[hH]*|--[hH]*) - printf "usage:\n" - printf " `basename $0` evt= r500= basedir= info= inner= outer= regin= regout= bkgd= nh= z= arf= rmf= [ grpcmd= log= ]\n" - printf "\nversion:\n" - printf "${VERSION}, ${UPDATE}\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 -COSCALC="`which cosmo_calc calc_distance 2>/dev/null | head -n 1`" -# COSCALC="_path_to_calc_distance_" -# COSCALC="$HOME/bin/mass/calc_distance" -if [ -z "${COSCALC}" ] || [ ! -x ${COSCALC} ]; then - printf "ERROR: \`COSCALC: ${COSCALC}' neither specified nor executable\n" - exit 255 -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`" -# default `bkgd', use `bkgcorr_blanksky*' corrected bkg spectrum -DFT_BKGD="`ls bkgcorr_blanksky_*.pi 2> /dev/null`" -# default basedir -DFT_BASEDIR="../.." -# default `radial region file' -#DFT_REG_IN="_NOT_EXIST_" -DFT_REG_IN="rspec.reg" -# default region range (0.1-0.5 R500) -DFT_INNER="0.1" -DFT_OUTER="0.5" -# default ARF/RMF, the one of the outmost region -DFT_ARF="`ls -1 r?_*.warf 2> /dev/null | tail -n 1`" -DFT_RMF="`ls -1 r?_*.wrmf 2> /dev/null | tail -n 1`" - -# default `group command' for `grppha' -DFT_GRP_CMD="group min 25" -#DFT_GRP_CMD="group 1 128 2 129 256 4 257 512 8 513 1024 16" -#DFT_GRP_CMD="group 1 128 4 129 256 8 257 512 16 513 1024 32" -# default `log file' -DFT_LOGFILE="r500avgt_`date '+%Y%m%d%H'`.log" - -# default JSON pattern -DFT_JSON_PAT="*_INFO.json" -## 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 -} - -## background renormalization (BACKSCAL) {{{ -# renorm background according to particle background -# energy range: 9.5-12.0 keV (channel: 651-822) -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 }'` - punlearn dmkeypar - EXPTIME=`dmkeypar $1 EXPOSURE echo=yes` - BACK=`dmkeypar $1 BACKSCAL echo=yes` - # fix `scientific notation' bug for `bc' - EXPTIME_B=`echo ${EXPTIME} | sed 's/[eE]/\*10\^/' | sed 's/+//'` - BACK_B=`echo "( ${BACK} )" | sed 's/[eE]/\*10\^/' | sed 's/+//'` - PB_FLUX=`echo "scale = 16; ${COUNTS} / ${EXPTIME_B} / ${BACK_B}" | bc -l` - echo ${PB_FLUX} -} - -bkg_renorm() { - # $1: src spectrum, $2: back spectrum - PBFLUX_SRC=`pb_flux $1` - PBFLUX_BKG=`pb_flux $2` - BACK_OLD=`dmkeypar $2 BACKSCAL echo=yes` - BACK_OLD_B=`echo "( ${BACK_OLD} )" | sed 's/[eE]/\*10\^/' | sed 's/+//'` - BACK_NEW=`echo "scale = 16; ${BACK_OLD_B} * ${PBFLUX_BKG} / ${PBFLUX_SRC}" | bc -l` - printf "\`$2': BACKSCAL:\n" - printf " ${BACK_OLD} --> ${BACK_NEW}\n" - punlearn dmhedit - dmhedit infile=$2 filelist=none operation=add \ - key=BACKSCAL value=${BACK_NEW} comment="old value: ${BACK_OLD}" -} -## bkg renorm }}} -## functions end }}} - -## 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 -if [ -d "${basedir}" ]; then - BASEDIR=${basedir} -else - BASEDIR=${DFT_BASEDIR} -fi -if [ ! -z "${json}" ] && [ -r "${BASEDIR}/${json}" ]; then - JSON_FILE="${BASEDIR}/${json}" -elif ls ${BASEDIR}/${DFT_JSON_PAT} > /dev/null 2>&1; then - JSON_FILE=`ls ${BASEDIR}/${DFT_JSON_PAT}` -else - read -p "> JSON_file: " JSON_FILE - if [ ! -r "${JSON_FILE}" ]; then - printf "ERROR: cannot access given \`${JSON_FILE}'\n" - exit ${ERR_JSON} - fi -fi -printf "## use json_file: \`${JSON_FILE}'\n" | ${TOLOG} - -# process `nh' and `redshift' {{{ -NH_JSON=`grep '"nH' ${JSON_FILE} | sed 's/.*"nH.*":\ //' | sed 's/\ *,$//'` -Z_JSON=`grep '"redshift' ${JSON_FILE} | sed 's/.*"redshift.*":\ //' | sed 's/\ *,$//'` -printf "## get nh: \`${NH_JSON}' (from \`${JSON_FILE}')\n" | ${TOLOG} -printf "## get redshift: \`${Z_JSON}' (from \`${JSON_FILE}')\n" | ${TOLOG} -## if `nh' and `redshift' supplied in cmdline, then use them -if [ ! -z "${nh}" ]; then - N_H=${nh} -else - N_H=${NH_JSON} -fi -# redshift -if [ ! -z "${z}" ]; then - REDSHIFT=${z} -else - REDSHIFT=${Z_JSON} -fi -printf "## use nH: ${N_H}\n" | ${TOLOG} -printf "## use redshift: ${REDSHIFT}\n" | ${TOLOG} -# nh & redshift }}} - -# region range {{{ -if [ ! -z "${inner}" ]; then - INNER="${inner}" -else - INNER=${DFT_INNER} -fi -if [ ! -z "${outer}" ]; then - OUTER="${outer}" -else - OUTER=${DFT_OUTER} -fi -printf "## region range: (${INNER} - ${OUTER} R500)\n" | ${TOLOG} -# range }}} - -# process `r500' {{{ -R500_RAW=`grep '"R500.*kpc' ${JSON_FILE} | sed 's/.*"R500.*":\ //' | sed 's/\ *,$//'` -if [ ! -z "${r500}" ]; then - R500_RAW=${r500} -fi -if [ -z "${R500_RAW}" ]; then - printf "## input R500 followed with unit, e.g.: 800kpc, 400pix\n" - read -p "> value of \`R500' (in pixel/kpc): " R500_RAW -fi -R500_VAL=`echo "${R500_RAW}" | tr -d 'a-zA-Z, '` -R500_UNI=`echo "${R500_RAW}" | tr -d '0-9, '` -printf "## get \`R500': ${R500_VAL} in unit \`${R500_UNI}'\n" | ${TOLOG} - -# if in kpc, convert to pix -case "${R500_UNI}" in - [pP]*) - printf "## units in \`pixel', conversion not needed\n" | ${TOLOG} - R500_PIX_B=`echo ${R500_VAL} | sed 's/[eE]/\*10\^/' | sed 's/+//'` - ;; - *) - printf "## units in \`kpc', convert to \`Chandra pixel'\n" | ${TOLOG} - 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" - R500_VAL_B=`echo ${R500_VAL} | sed 's/[eE]/\*10\^/' | sed 's/+//'` - R500_PIX_B=`echo "scale = 4; ${R500_VAL_B} / ( ${KPC_PER_PIX_B} )" | bc -l` - ;; -esac -# calc (inner-outer R500) -R_IN=`echo "scale = 4; ${INNER} * ${R500_PIX_B}" | bc -l` -R_OUT=`echo "scale = 4; ${OUTER} * ${R500_PIX_B}" | bc -l` -printf "## R500 in units pixel: ${R500_PIX_B}\n" | ${TOLOG} -printf "## (${INNER}-${OUTER} R500) range in pixel: ${R_IN} - ${R_OUT}\n" | ${TOLOG} -# r500 }}} - -# 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} - -# input and output region files {{{ -if [ -r "${regin}" ]; then - REG_IN="${regin}" -elif [ -r "${DFT_REG_IN}" ]; then - REG_IN=${DFT_REG_IN} -else - read -p "> previous used radial spec regfile: " REG_IN - if [ ! -r "${REG_IN}" ]; then - printf "ERROR: cannot access given \`${REG_IN}' region file\n" - exit ${ERR_REG} - fi -fi -printf "## use previous regfile: \`${REG_IN}'\n" | ${TOLOG} -if [ ! -z "${regout}" ]; then - REG_OUT="${regout}" -else - REG_OUT="r500avgt_${INNER}-${OUTER}.reg" -fi -[ -e "${REG_OUT}" ] && mv -fv ${REG_OUT} ${REG_OUT}_bak -printf "## set output regfile: \`${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` -XC=`echo ${TMP_REG} | tr -d 'a-zA-Z() ' | awk -F',' '{ print $1 }'` -YC=`echo ${TMP_REG} | tr -d 'a-zA-Z() ' | awk -F',' '{ print $2 }'` -printf "## get center coord: (${XC},${YC})\n" | ${TOLOG} -# region files }}} - -# check given bkgd, determine background {{{ -if [ -r "${bkgd}" ]; then - BKGD=${bkgd} -elif [ -r "${DFT_BKGD}" ]; then - BKGD=${DFT_BKGD} -else - read -p "> background (blanksky_evt | lbkg_reg | bkg_spec): " BKGD - if [ ! -r "${BKGD}" ]; then - printf "ERROR: cannot access given \`${BKGD}'\n" - exit ${ERR_BKG} - fi -fi -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 - printf "## given \`${BKGD}' is a \`text file'\n" - printf "## use it as local bkg region file\n" - printf "## use *LOCAL BKG SPEC*\n" | ${TOLOG} - # just set flags, extract spectrum later - USE_LBKG_REG=YES - USE_BLANKSKY=NO - USE_BKG_SPEC=NO -elif file -bL ${BKGD} | grep -qi 'FITS'; then - printf "## given \`${BKGD}' is a \`FITS file'\n" - # get FITS header keyword - HDUCLAS1=`dmkeypar ${BKGD} HDUCLAS1 echo=yes` - if [ "${HDUCLAS1}" = "EVENTS" ]; then - # event file - printf "## given file is \`event'\n" - # check if `blanksky' or `stowed bkg' - BKG_OBJ=`dmkeypar ${BKGD} OBJECT echo=yes` - if [ "${BKG_OBJ}" = "BACKGROUND DATASET" ] || [ "${BKG_OBJ}" = "ACIS STOWED" ]; then - # valid bkg evt file - printf "## given FITS file is a valid bkgrnd file\n" - printf "## use *BLANKSKY*\n" | ${TOLOG} - USE_BLANKSKY=YES - USE_LBKG_REG=NO - USE_BKG_SPEC=NO - # specify `BLANKSKY' - BLANKSKY=${BKGD} - else - # invalid bkg evt file - printf "ERROR: invalid bkg evt file given\n" - exit ${ERR_BKGTY} - fi - elif [ "${HDUCLAS1}" = "SPECTRUM" ]; then - # spectrum file - printf "## given file is \`spectrum'\n" - printf "## use *BKG SPECTRUM*\n" | ${TOLOG} - USE_BKG_SPEC=YES - USE_BLANKSKY=NO - USE_LBKG_REG=NO - # specify `BKG_SPEC' - BKG_SPEC=${BKGD} - else - # other type - printf "ERROR: other type FITS given\n" - exit ${ERR_BKGTY} - fi -else - printf "ERROR: given \`${BKGD}' type UNKNOWN\n" - exit ${ERR_BKGTY} -fi -# bkgd }}} - -# check `arf' and `rmf' {{{ -if [ -r "${arf}" ]; then - ARF=${arf} -elif [ -r "${DFT_ARF}" ]; then - ARF=${DFT_ARF} -else - read -p "> provide the ARF to use: " ARF - if [ ! -r "${ARF}" ]; then - printf "ERROR: cannot access given \`${ARF}'\n" - exit ${ERR_ARF} - fi -fi -printf "## use ARF: \`${ARF}'\n" | ${TOLOG} -# rmf -if [ -r "${rmf}" ]; then - RMF=${rmf} -elif [ -r "${DFT_RMF}" ]; then - RMF=${DFT_RMF} -else - read -p "> provide the RMF to use: " RMF - if [ ! -r "${RMF}" ]; then - printf "ERROR: cannot access given \`${RMF}'\n" - exit ${ERR_RMF} - fi -fi -printf "## use RMF: \`${RMF}'\n" | ${TOLOG} -# arf & rmf }}} - -# check given `grpcmd' -if [ ! -z "${grpcmd}" ]; then - GRP_CMD="${grpcmd}" -else - GRP_CMD="${DFT_GRP_CMD}" -fi -printf "## use grppha cmd: \`${GRP_CMD}'\n" | ${TOLOG} -## parameters }}} - -################################################## -#### main -## region related {{{ -## generate the needed region file -printf "generate the output region file ...\n" -cat > ${REG_OUT} << _EOF_ -# Region file format: CIAO version 1.0 -pie(${XC},${YC},${R_IN},${R_OUT},0,360) -_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} -region ${REG_OUT} - -## 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'` -if [ "x${INVALID}" != "x" ]; then - printf "*** WARNING: there are pie regions' END_ANGLE > 360\n" | ${TOLOG} - printf "*** will to fix ...\n" - mv -fv ${REG_OUT} ${REG_OUT}_tmp - # using `awk' to fix - awk -F'[,()]' '{ - if ($7 > 360) { - printf "%s(%.2f,%.2f,%.2f,%.2f,%.2f,%.2f)\n", $1,$2,$3,$4,$5,$6,($7-360) - } - else { - print $0 - } - }' ${REG_OUT}_tmp > ${REG_OUT} - rm -f ${REG_OUT}_tmp -fi -## region related }}} - -## generate spectrum {{{ -# object -AVGT_SPEC="${REG_OUT%.reg}.pi" -AVGT_SPEC_GRP="${AVGT_SPEC%.pi}_grp.pi" -printf "extract object spectrum \`${AVGT_SPEC}' ...\n" -punlearn dmextract -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 - -# background -printf "generate the background spectrum ...\n" -AVGT_BKG="${AVGT_SPEC%.pi}_bkg.pi" -if [ "${USE_BLANKSKY}" = "YES" ]; then - # use blanksky as background file - printf "extract spectrum from blanksky ...\n" - punlearn dmextract - dmextract infile="${BLANKSKY}[sky=region(${REG_OUT})][bin PI]" \ - outfile=${AVGT_BKG} wmap="[bin det=8]" clobber=yes -elif [ "${USE_LBKG_REG}" = "YES" ]; then - printf "extract local background ...\n" - punlearn dmextract - dmextract infile="${EVT}[sky=region(${BKGD})][bin PI]" \ - outfile=${AVGT_BKG} wmap="[bin det=8]" clobber=yes -elif [ "${USE_BKG_SPEC}" = "YES" ]; then - printf "copy specified background spectrum ...\n" - cp -fv ${BKG_SPEC} ${AVGT_BKG} -fi - -printf "renormalize the background ...\n" -bkg_renorm ${AVGT_SPEC} ${AVGT_BKG} - -## spectrum }}} - -## generate XSPEC script {{{ -printf "generate a XSPEC script ...\n" -# default output xspec scripts -XSPEC_SCRIPT="xspec_${REG_OUT%.reg}.xcm" -[ -e "${XSPEC_SCRIPT}" ] && mv -fv ${XSPEC_SCRIPT} ${XSPEC_SCRIPT}_bak -cat > ${XSPEC_SCRIPT} << _EOF_ -## XSPEC script -## spectrum analysis to get the average temperatue with (0.1-0.5 R500) -## -## generated by: \``basename $0`' -## date: \``date`' -## - -# xspec settings -statistic chi -abund grsa -query yes - -# data -data ${AVGT_SPEC_GRP} -response ${RMF} -arf ${ARF} -backgrnd ${AVGT_BKG} - -# fitting range -ignore bad -ignore 0.0-0.7,7.0-** - -# plot related -setplot energy - -method leven 10 0.01 -xsect bcmc -cosmo 70 0 0.73 -xset delta 0.01 -systematic 0 - -# model -model wabs*apec - ${N_H} -0.001 0 0 100000 1e+06 - 1.0 0.01 0.008 0.008 64 64 - 0.4 0.001 0 0 5 5 - ${REDSHIFT} -0.01 -0.999 -0.999 10 10 - 0.0 0.01 0 0 1e+24 1e+24 - -## xspec script end -_EOF_ -## xspec script }}} - -exit 0 - diff --git a/scripts/ciao_sbp.sh b/scripts/ciao_sbp.sh new file mode 100755 index 0000000..d111b7f --- /dev/null +++ b/scripts/ciao_sbp.sh @@ -0,0 +1,529 @@ +#!/bin/sh +# +unalias -a +export LC_COLLATE=C +########################################################### +## 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 ## +## August 16, 2012 ## +########################################################### + +########################################################### +## Changes made by zzh: (2013-02-01) +## removes the region in ccd gap of ACIS_I +## removes the region in the area of point source +## need asol file to prevent offset +########################################################### + +SCRIPT_PATH=`readlink -f $0` +SCRIPT_DIR=`dirname ${SCRIPT_PATH}` +CCDGAP_SCRIPT="chandra_ccdgap_rect.py" + +## about, used in `usage' {{{ +VERSION="v3.1" +UPDATE="2013-02-01" +## about }}} + +## error code {{{ +ERR_USG=1 +ERR_DIR=11 +ERR_EVT=12 +ERR_BKG=13 +ERR_REG=14 +ERR_CELL_REG=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_e= reg= expmap= cellreg= aspec= [bkg= log= ]\n" + printf "\nversion:\n" + printf "${VERSION}, ${UPDATE}\n" + exit ${ERR_USG} + ;; +esac +## usage, help }}} + +## default parameters {{{ +# default energy band +E_RANGE="700:7000" +# default ccd edge cut pixel (20pixel) +DFT_CCD_EDGECUT=25 +# default `event file' which used to match `blanksky' files +#DFT_EVT="_NOT_EXIST_" +DFT_EVT_E="`ls evt2_*_e*.fits 2> /dev/null`" +# default expmap +DFT_EXPMAP="`ls expmap*.fits 2> /dev/null | head -n 1`" +# default `radial region file' to extract surface brightness +#DFT_SBP_REG="_NOT_EXIST_" +DFT_SBP_REG="sbprofile.reg" +# defalut pointsource region file +DFT_CELL_REG="`ls celld*.reg 2> /dev/null`" +# defalut asol file +DFT_ASOL_FILE="`ls pcad*asol*fits 2> /dev/null`" + +# default `log file' +DFT_LOGFILE="sbp_`date '+%Y%m%d'`.log" + +## 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_e}" ]; then + EVT_E=${evt_e} +elif [ -r "${DFT_EVT_E}" ]; then + EVT_E=${DFT_EVT_E} +else + read -p "clean evt2 file: " EVT_E + if [ ! -r "${EVT_E}" ]; then + printf "ERROR: cannot access given \`${EVT_E}' evt file\n" + exit ${ERR_EVT} + fi +fi +printf "## use evt_eng file: \`${EVT_E}'\n" | ${TOLOG} +# check evt file +if [ -r "${expmap}" ]; then + EXPMAP=${expmap} +elif [ -r "${DFT_EXPMAP}" ]; then + EXPMAP=${DFT_EXPMAP} +else + read -p "expocure map file: " EXPMAP + if [ ! -r "${EXPMAP}" ]; then + printf "ERROR: cannot access given \`${EXPMAP}' expmap file\n" + exit ${ERR_EVT} + fi +fi +printf "## use expmap file: \`${EXPMAP}'\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 bkg +if [ ! -z "${bkg}" ]; then + BKG=${bkg} +#else +# read -p "> bkg: " BKG +fi +if [ -r "${BKG}" ]; then + printf "## use bkg: \`${BKG}'\n" | ${TOLOG} +else + BKG="NULL" +fi +# check cell region file +if [ -r "${cellreg}" ]; then + CELL_REG="${cellreg}" +elif [ -r "${DFT_CELL_REG}" ] ; then + CELL_REG="${DFT_CELL_REG}" +else + read -p "> celldetect region file: " CELL_REG + if [ ! -r "${CELL_REG}" ]; then + printf "ERROR: cannot access given \`${CELL_REG}' region file \n" + exit ${ERR_CELL_REG} + fi +fi +printf "## use cell reg file(s): \`${CELL_REG}'\n" | ${TOLOG} +# check asol file +if [ -r "${aspec}" ]; then + ASOL_FILE="${aspec}" +elif [ -r "${DFT_ASOL_FILE}" ] ; then + ASOL_FILE="${DFT_ASOL_FILE}" +else + read -p ">asol file: " ASOL_FILE + if [ ! -r "${ASOL_FILE}" ] ; then + printf " ERROR: cannot access asol file \n" + exit ${ERR_ASOL} + fi +fi +printf "## use asol file(s) : \`${ASOL_FILE}'\n" | ${TOLOG} +## parameters }}} + +## determine ACIS type {{{ +# consistent with `ciao_procevt' +punlearn dmkeypar +DETNAM=`dmkeypar ${EVT_E} DETNAM echo=yes` +if echo ${DETNAM} | grep -q 'ACIS-0123'; then + printf "## \`DETNAM' (${DETNAM}) has chips 0123 -> ACIS-I\n" + ACIS_TYPE="I" + CCD="0:3" +elif echo ${DETNAM} | grep -q 'ACIS-[0-6]*7'; then + printf "## \`DETNAM' (${DETNAM}) has chip 7 -> ACIS-S\n" + ACIS_TYPE="S" + CCD="7" +else + printf "ERROR: unknown detector type: ${DETNAM}\n" + exit ${ERR_DET} +fi +## ACIS type }}} + +## check validity of pie region {{{ +INVALID=`grep -i 'pie' ${SBP_REG} | awk -F'[,()]' '$7 > 360'` +SBP_REG_FIX="_${SBP_REG%.reg}_fix.reg" +if [ "x${INVALID}" != "x" ]; then + printf "*** WARNING: some pie regions' END_ANGLE > 360\n" | ${TOLOG} + printf "*** script will fix ...\n" + cp -fv ${SBP_REG} ${SBP_REG_FIX} + # using `awk' to fix + awk -F'[,()]' '{ + if ($7 > 360) { + printf "%s(%.2f,%.2f,%.2f,%.2f,%.2f,%.2f)\n", $1,$2,$3,$4,$5,$6,($7-360) + } + else { + print $0 + } + }' ${SBP_REG} | sed '/^\#/d' > ${SBP_REG_FIX} +else + cat ${SBP_REG} | sed '/^\#/d' >${SBP_REG_FIX} +fi +## pie validity }}} + +## main process {{{ + +## 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" +punlearn skyfov +skyfov infile="${EVT_E}" outfile="${SKYFOV}" aspec="${ASOL_FILE}" clobber=yes + +## get CCD fov regions {{{ +printf "make regions for CCD ...\n" +TMP_LIST="_tmp_list.txt" +TMP_REC="_tmp_rec.reg" +if [ "${ACIS_TYPE}" = "S" ]; then + # ACIS-S + punlearn dmlist + dmlist infile="${SKYFOV}[ccd_id=${CCD}][cols POS]" opt="data,clean" | awk '{for (i=1;i<=NF;i++) print $i }' |sed -e ':a;N;s/\n/,/;ta' | awk -F"]," '{print "polygon("$2}' | awk -F"NaN" '{print $1}' >${TMP_LIST} + python ${SCRIPT_DIR}/${CCDGAP_SCRIPT} ${TMP_LIST} >${TMP_REC} + XC=` cat ${TMP_REC} | awk -F\( '{print $2}' |awk -F\) '{print $1}' |awk -F\, '{print $1}'` + YC=` cat ${TMP_REC} | awk -F\( '{print $2}' |awk -F\) '{print $1}' |awk -F\, '{print $2}'` + ADD_L=` cat ${TMP_REC} | awk -F\( '{print $2}' |awk -F\) '{print $1}' |awk -F\, '{print $3/2}'` + ANG=` cat ${TMP_REC} | awk -F\( '{print $2}' |awk -F\) '{print $1}' |awk -F\, '{print $5}'` + while [ 1 -eq 1 ]; do + if [ `echo "${ANG} < 0" |bc -l ` -eq 1 ] ; then + ANG=` echo " ${ANG} + 90 " | bc -l ` + elif [ `echo "${ANG} >=90" |bc -l ` -eq 1 ] ; then + ANG=` echo " ${ANG} - 90 " | bc -l` + else + break + fi + done + ANG=`echo "${ANG}/180*3.1415926" |bc -l` + CCD_1_X_RAW=` echo " ${XC} ${ADD_L} ${ANG} "| awk '{print $1-$2*cos($3)-$2*sin($3)}' ` + CCD_2_X_RAW=` echo " ${XC} ${ADD_L} ${ANG} "| awk '{print $1+$2*cos($3)-$2*sin($3)}' ` + CCD_3_X_RAW=` echo " ${XC} ${ADD_L} ${ANG} "| awk '{print $1-$2*cos($3)+$2*sin($3)}' ` + CCD_4_X_RAW=` echo " ${XC} ${ADD_L} ${ANG} "| awk '{print $1+$2*cos($3)+$2*sin($3)}' ` + CCD_1_Y_RAW=` echo " ${YC} ${ADD_L} ${ANG} "| awk '{print $1+$2*cos($3)-$2*sin($3)}' ` + CCD_2_Y_RAW=` echo " ${YC} ${ADD_L} ${ANG} "| awk '{print $1+$2*cos($3)+$2*sin($3)}' ` + CCD_3_Y_RAW=` echo " ${YC} ${ADD_L} ${ANG} "| awk '{print $1-$2*cos($3)-$2*sin($3)}' ` + CCD_4_Y_RAW=` echo " ${YC} ${ADD_L} ${ANG} "| awk '{print $1-$2*cos($3)+$2*sin($3)}' ` + CCD_1_RAW=` echo "${CCD_1_X_RAW},${CCD_1_Y_RAW}"` + CCD_2_RAW=` echo "${CCD_2_X_RAW},${CCD_2_Y_RAW}"` + CCD_3_RAW=` echo "${CCD_3_X_RAW},${CCD_3_Y_RAW}"` + CCD_4_RAW=` echo "${CCD_4_X_RAW},${CCD_4_Y_RAW}"` + REG_CCD_RAW="`echo "polygon(${CCD_1_RAW}, ${CCD_2_RAW}, ${CCD_4_RAW}, ${CCD_3_RAW}) " `" + DX_2T1=$(echo "`echo ${CCD_2_RAW} | awk -F\, '{print $1}'`-`echo ${CCD_1_RAW} |awk -F\, '{print $1}'`" |bc -l) + DY_2T1=$(echo "`echo ${CCD_2_RAW} | awk -F\, '{print $2}'`-`echo ${CCD_1_RAW} |awk -F\, '{print $2}'`" |bc -l) + CCD_WIDTH=`echo "sqrt(${DX_2T1}*${DX_2T1}+${DY_2T1}*${DY_2T1})" | bc -l` + CCD_2T1_MOV_X=`echo "${DFT_CCD_EDGECUT}*${DX_2T1}/${CCD_WIDTH}" | bc -l` + CCD_2T1_MOV_Y=`echo "${DFT_CCD_EDGECUT}*${DY_2T1}/${CCD_WIDTH}" | bc -l` + DX_3T1=$(echo "`echo ${CCD_3_RAW} | awk -F\, '{print $1}'`-`echo ${CCD_1_RAW} |awk -F\, '{print $1}'`" |bc -l) + DY_3T1=$(echo "`echo ${CCD_3_RAW} | awk -F\, '{print $2}'`-`echo ${CCD_1_RAW} |awk -F\, '{print $2}'`" |bc -l) + CCD_3T1_MOV_X=`echo "${DFT_CCD_EDGECUT}*${DX_3T1}/${CCD_WIDTH}" | bc -l` + CCD_3T1_MOV_Y=`echo "${DFT_CCD_EDGECUT}*${DY_3T1}/${CCD_WIDTH}" | bc -l` + CCD_1_X=$(echo "` echo ${CCD_1_RAW} |awk -F\, '{print $1}' ` + `echo ${CCD_2T1_MOV_X}` +`echo ${CCD_3T1_MOV_X}` "| bc -l) + CCD_1_Y=$(echo "` echo ${CCD_1_RAW} |awk -F\, '{print $2}' ` + `echo ${CCD_2T1_MOV_Y}` +`echo ${CCD_3T1_MOV_Y}` "| bc -l) + DX_1T2=$(echo "`echo ${CCD_1_RAW} | awk -F\, '{print $1}'`-`echo ${CCD_2_RAW} |awk -F\, '{print $1}'`" |bc -l) + DY_1T2=$(echo "`echo ${CCD_1_RAW} | awk -F\, '{print $2}'`-`echo ${CCD_2_RAW} |awk -F\, '{print $2}'`" |bc -l) + CCD_1T2_MOV_X=`echo "${DFT_CCD_EDGECUT}*${DX_1T2}/${CCD_WIDTH}" | bc -l` + CCD_1T2_MOV_Y=`echo "${DFT_CCD_EDGECUT}*${DY_1T2}/${CCD_WIDTH}" | bc -l` + DX_4T2=$(echo "`echo ${CCD_4_RAW} | awk -F\, '{print $1}'`-`echo ${CCD_2_RAW} |awk -F\, '{print $1}'`" |bc -l) + DY_4T2=$(echo "`echo ${CCD_4_RAW} | awk -F\, '{print $2}'`-`echo ${CCD_2_RAW} |awk -F\, '{print $2}'`" |bc -l) + CCD_4T2_MOV_X=`echo "${DFT_CCD_EDGECUT}*${DX_4T2}/${CCD_WIDTH}" | bc -l` + CCD_4T2_MOV_Y=`echo "${DFT_CCD_EDGECUT}*${DY_4T2}/${CCD_WIDTH}" | bc -l` + CCD_2_X=$(echo "` echo ${CCD_2_RAW} |awk -F\, '{print $1}' ` + `echo ${CCD_1T2_MOV_X}` +`echo ${CCD_4T2_MOV_X}` "| bc -l) + CCD_2_Y=$(echo "` echo ${CCD_2_RAW} |awk -F\, '{print $2}' ` + `echo ${CCD_1T2_MOV_Y}` +`echo ${CCD_4T2_MOV_Y}` "| bc -l) + DX_1T3=$(echo "`echo ${CCD_1_RAW} | awk -F\, '{print $1}'`-`echo ${CCD_3_RAW} |awk -F\, '{print $1}'`" |bc -l) + DY_1T3=$(echo "`echo ${CCD_1_RAW} | awk -F\, '{print $2}'`-`echo ${CCD_3_RAW} |awk -F\, '{print $2}'`" |bc -l) + CCD_1T3_MOV_X=`echo "${DFT_CCD_EDGECUT}*${DX_1T3}/${CCD_WIDTH}" | bc -l` + CCD_1T3_MOV_Y=`echo "${DFT_CCD_EDGECUT}*${DY_1T3}/${CCD_WIDTH}" | bc -l` + DX_4T3=$(echo "`echo ${CCD_4_RAW} | awk -F\, '{print $1}'`-`echo ${CCD_3_RAW} |awk -F\, '{print $1}'`" |bc -l) + DY_4T3=$(echo "`echo ${CCD_4_RAW} | awk -F\, '{print $2}'`-`echo ${CCD_3_RAW} |awk -F\, '{print $2}'`" |bc -l) + CCD_4T3_MOV_X=`echo "${DFT_CCD_EDGECUT}*${DX_4T3}/${CCD_WIDTH}" | bc -l` + CCD_4T3_MOV_Y=`echo "${DFT_CCD_EDGECUT}*${DY_4T3}/${CCD_WIDTH}" | bc -l` + CCD_3_X=$(echo "` echo ${CCD_3_RAW} |awk -F\, '{print $1}' ` + `echo ${CCD_1T3_MOV_X}` +`echo ${CCD_4T3_MOV_X}` "| bc -l) + CCD_3_Y=$(echo "` echo ${CCD_3_RAW} |awk -F\, '{print $2}' ` + `echo ${CCD_1T3_MOV_Y}` +`echo ${CCD_4T3_MOV_Y}` "| bc -l) + DX_2T4=$(echo "`echo ${CCD_2_RAW} | awk -F\, '{print $1}'`-`echo ${CCD_4_RAW} |awk -F\, '{print $1}'`" |bc -l) + DY_2T4=$(echo "`echo ${CCD_2_RAW} | awk -F\, '{print $2}'`-`echo ${CCD_4_RAW} |awk -F\, '{print $2}'`" |bc -l) + CCD_2T4_MOV_X=`echo "${DFT_CCD_EDGECUT}*${DX_2T4}/${CCD_WIDTH}" | bc -l` + CCD_2T4_MOV_Y=`echo "${DFT_CCD_EDGECUT}*${DY_2T4}/${CCD_WIDTH}" | bc -l` + DX_3T4=$(echo "`echo ${CCD_3_RAW} | awk -F\, '{print $1}'`-`echo ${CCD_4_RAW} |awk -F\, '{print $1}'`" |bc -l) + DY_3T4=$(echo "`echo ${CCD_3_RAW} | awk -F\, '{print $2}'`-`echo ${CCD_4_RAW} |awk -F\, '{print $2}'`" |bc -l) + CCD_3T4_MOV_X=`echo "${DFT_CCD_EDGECUT}*${DX_3T4}/${CCD_WIDTH}" | bc -l` + CCD_3T4_MOV_Y=`echo "${DFT_CCD_EDGECUT}*${DY_3T4}/${CCD_WIDTH}" | bc -l` + CCD_4_X=$(echo "` echo ${CCD_4_RAW} |awk -F\, '{print $1}' ` + `echo ${CCD_2T4_MOV_X}` +`echo ${CCD_3T4_MOV_X}` "| bc -l) + CCD_4_Y=$(echo "` echo ${CCD_4_RAW} |awk -F\, '{print $2}' ` + `echo ${CCD_2T4_MOV_Y}` +`echo ${CCD_3T4_MOV_Y}` "| bc -l) + REG_CCD_CUT=`echo "polygon(${CCD_1_X},${CCD_1_Y},${CCD_2_X},${CCD_2_Y},${CCD_4_X},${CCD_4_Y},${CCD_3_X},${CCD_3_Y})"` + REG_FILE_CCD="_ccd.reg" + [ -e "${REG_FILE_CCD}" ] && mv -f ${REG_FILE_CCD} ${REG_FILE_CCD}_bak + echo "${REG_CCD_CUT}" >>${REG_FILE_CCD} + +elif [ "${ACIS_TYPE}" = "I" ]; then + # ACIS-I + TMP_REG_FILE_CCD="_ccd_tmp.reg" + [ -e "${TMP_REG_FILE_CCD}" ] && mv -f ${TMP_REG_FILE_CCD} ${TMP_REG_FILE_CCD}_bak + for i in `seq 0 3` ; do + punlearn dmlist + dmlist infile="${SKYFOV}[ccd_id=${i}][cols POS]" opt="data,clean" | awk '{for (i=1;i<=NF;i++) print $i }' |sed -e ':a;N;s/\n/,/;ta' | awk -F"]," '{print "polygon("$2}' | awk -F"NaN" '{print $1}' >${TMP_LIST} + python ${SCRIPT_DIR}/${CCDGAP_SCRIPT} ${TMP_LIST} >${TMP_REC} + XC=` cat ${TMP_REC} | awk -F\( '{print $2}' |awk -F\) '{print $1}' |awk -F\, '{print $1}'` + YC=` cat ${TMP_REC} | awk -F\( '{print $2}' |awk -F\) '{print $1}' |awk -F\, '{print $2}'` + ADD_L=` cat ${TMP_REC} | awk -F\( '{print $2}' |awk -F\) '{print $1}' |awk -F\, '{print $3/2}'` + ANG=` cat ${TMP_REC} | awk -F\( '{print $2}' |awk -F\) '{print $1}' |awk -F\, '{print $5}'` + while [ 1 -eq 1 ]; do + if [ `echo "${ANG} < 0" |bc -l ` -eq 1 ] ; then + ANG=` echo " ${ANG} + 90 " | bc -l ` + elif [ `echo "${ANG} >=90" |bc -l ` -eq 1 ] ; then + ANG=` echo " ${ANG} - 90 " | bc -l` + else + break + fi + done + ANG=`echo "${ANG}/180*3.1415926" |bc -l` + CCD_1_X_RAW=` echo " ${XC} ${ADD_L} ${ANG} "| awk '{print $1-$2*cos($3)-$2*sin($3)}' ` + CCD_2_X_RAW=` echo " ${XC} ${ADD_L} ${ANG} "| awk '{print $1+$2*cos($3)-$2*sin($3)}' ` + CCD_3_X_RAW=` echo " ${XC} ${ADD_L} ${ANG} "| awk '{print $1-$2*cos($3)+$2*sin($3)}' ` + CCD_4_X_RAW=` echo " ${XC} ${ADD_L} ${ANG} "| awk '{print $1+$2*cos($3)+$2*sin($3)}' ` + CCD_1_Y_RAW=` echo " ${YC} ${ADD_L} ${ANG} "| awk '{print $1+$2*cos($3)-$2*sin($3)}' ` + CCD_2_Y_RAW=` echo " ${YC} ${ADD_L} ${ANG} "| awk '{print $1+$2*cos($3)+$2*sin($3)}' ` + CCD_3_Y_RAW=` echo " ${YC} ${ADD_L} ${ANG} "| awk '{print $1-$2*cos($3)-$2*sin($3)}' ` + CCD_4_Y_RAW=` echo " ${YC} ${ADD_L} ${ANG} "| awk '{print $1-$2*cos($3)+$2*sin($3)}' ` + CCD_1_RAW=` echo "${CCD_1_X_RAW},${CCD_1_Y_RAW}"` + CCD_2_RAW=` echo "${CCD_2_X_RAW},${CCD_2_Y_RAW}"` + CCD_3_RAW=` echo "${CCD_3_X_RAW},${CCD_3_Y_RAW}"` + CCD_4_RAW=` echo "${CCD_4_X_RAW},${CCD_4_Y_RAW}"` + REG_CCD_RAW=`echo "polygon(${CCD_1_RAW}, ${CCD_2_RAW}, ${CCD_4_RAW}, ${CCD_3_RAW}) " ` + DX_2T1=$(echo "`echo ${CCD_2_RAW} | awk -F\, '{print $1}'`-`echo ${CCD_1_RAW} |awk -F\, '{print $1}'`" |bc -l) + DY_2T1=$(echo "`echo ${CCD_2_RAW} | awk -F\, '{print $2}'`-`echo ${CCD_1_RAW} |awk -F\, '{print $2}'`" |bc -l) + CCD_WIDTH=`echo "sqrt(${DX_2T1}*${DX_2T1}+${DY_2T1}*${DY_2T1})" | bc -l` + CCD_2T1_MOV_X=`echo "${DFT_CCD_EDGECUT}*${DX_2T1}/${CCD_WIDTH}" | bc -l` + CCD_2T1_MOV_Y=`echo "${DFT_CCD_EDGECUT}*${DY_2T1}/${CCD_WIDTH}" | bc -l` + DX_3T1=$(echo "`echo ${CCD_3_RAW} | awk -F\, '{print $1}'`-`echo ${CCD_1_RAW} |awk -F\, '{print $1}'`" |bc -l) + DY_3T1=$(echo "`echo ${CCD_3_RAW} | awk -F\, '{print $2}'`-`echo ${CCD_1_RAW} |awk -F\, '{print $2}'`" |bc -l) + CCD_3T1_MOV_X=`echo "${DFT_CCD_EDGECUT}*${DX_3T1}/${CCD_WIDTH}" | bc -l` + CCD_3T1_MOV_Y=`echo "${DFT_CCD_EDGECUT}*${DY_3T1}/${CCD_WIDTH}" | bc -l` + CCD_1_X=$(echo "` echo ${CCD_1_RAW} |awk -F\, '{print $1}' ` + `echo ${CCD_2T1_MOV_X}` +`echo ${CCD_3T1_MOV_X}` "| bc -l) + CCD_1_Y=$(echo "` echo ${CCD_1_RAW} |awk -F\, '{print $2}' ` + `echo ${CCD_2T1_MOV_Y}` +`echo ${CCD_3T1_MOV_Y}` "| bc -l) + DX_1T2=$(echo "`echo ${CCD_1_RAW} | awk -F\, '{print $1}'`-`echo ${CCD_2_RAW} |awk -F\, '{print $1}'`" |bc -l) + DY_1T2=$(echo "`echo ${CCD_1_RAW} | awk -F\, '{print $2}'`-`echo ${CCD_2_RAW} |awk -F\, '{print $2}'`" |bc -l) + CCD_1T2_MOV_X=`echo "${DFT_CCD_EDGECUT}*${DX_1T2}/${CCD_WIDTH}" | bc -l` + CCD_1T2_MOV_Y=`echo "${DFT_CCD_EDGECUT}*${DY_1T2}/${CCD_WIDTH}" | bc -l` + DX_4T2=$(echo "`echo ${CCD_4_RAW} | awk -F\, '{print $1}'`-`echo ${CCD_2_RAW} |awk -F\, '{print $1}'`" |bc -l) + DY_4T2=$(echo "`echo ${CCD_4_RAW} | awk -F\, '{print $2}'`-`echo ${CCD_2_RAW} |awk -F\, '{print $2}'`" |bc -l) + CCD_4T2_MOV_X=`echo "${DFT_CCD_EDGECUT}*${DX_4T2}/${CCD_WIDTH}" | bc -l` + CCD_4T2_MOV_Y=`echo "${DFT_CCD_EDGECUT}*${DY_4T2}/${CCD_WIDTH}" | bc -l` + CCD_2_X=$(echo "` echo ${CCD_2_RAW} |awk -F\, '{print $1}' ` + `echo ${CCD_1T2_MOV_X}` +`echo ${CCD_4T2_MOV_X}` "| bc -l) + CCD_2_Y=$(echo "` echo ${CCD_2_RAW} |awk -F\, '{print $2}' ` + `echo ${CCD_1T2_MOV_Y}` +`echo ${CCD_4T2_MOV_Y}` "| bc -l) + DX_1T3=$(echo "`echo ${CCD_1_RAW} | awk -F\, '{print $1}'`-`echo ${CCD_3_RAW} |awk -F\, '{print $1}'`" |bc -l) + DY_1T3=$(echo "`echo ${CCD_1_RAW} | awk -F\, '{print $2}'`-`echo ${CCD_3_RAW} |awk -F\, '{print $2}'`" |bc -l) + CCD_1T3_MOV_X=`echo "${DFT_CCD_EDGECUT}*${DX_1T3}/${CCD_WIDTH}" | bc -l` + CCD_1T3_MOV_Y=`echo "${DFT_CCD_EDGECUT}*${DY_1T3}/${CCD_WIDTH}" | bc -l` + DX_4T3=$(echo "`echo ${CCD_4_RAW} | awk -F\, '{print $1}'`-`echo ${CCD_3_RAW} |awk -F\, '{print $1}'`" |bc -l) + DY_4T3=$(echo "`echo ${CCD_4_RAW} | awk -F\, '{print $2}'`-`echo ${CCD_3_RAW} |awk -F\, '{print $2}'`" |bc -l) + CCD_4T3_MOV_X=`echo "${DFT_CCD_EDGECUT}*${DX_4T3}/${CCD_WIDTH}" | bc -l` + CCD_4T3_MOV_Y=`echo "${DFT_CCD_EDGECUT}*${DY_4T3}/${CCD_WIDTH}" | bc -l` + CCD_3_X=$(echo "` echo ${CCD_3_RAW} |awk -F\, '{print $1}' ` + `echo ${CCD_1T3_MOV_X}` +`echo ${CCD_4T3_MOV_X}` "| bc -l) + CCD_3_Y=$(echo "` echo ${CCD_3_RAW} |awk -F\, '{print $2}' ` + `echo ${CCD_1T3_MOV_Y}` +`echo ${CCD_4T3_MOV_Y}` "| bc -l) + DX_2T4=$(echo "`echo ${CCD_2_RAW} | awk -F\, '{print $1}'`-`echo ${CCD_4_RAW} |awk -F\, '{print $1}'`" |bc -l) + DY_2T4=$(echo "`echo ${CCD_2_RAW} | awk -F\, '{print $2}'`-`echo ${CCD_4_RAW} |awk -F\, '{print $2}'`" |bc -l) + CCD_2T4_MOV_X=`echo "${DFT_CCD_EDGECUT}*${DX_2T4}/${CCD_WIDTH}" | bc -l` + CCD_2T4_MOV_Y=`echo "${DFT_CCD_EDGECUT}*${DY_2T4}/${CCD_WIDTH}" | bc -l` + DX_3T4=$(echo "`echo ${CCD_3_RAW} | awk -F\, '{print $1}'`-`echo ${CCD_4_RAW} |awk -F\, '{print $1}'`" |bc -l) + DY_3T4=$(echo "`echo ${CCD_3_RAW} | awk -F\, '{print $2}'`-`echo ${CCD_4_RAW} |awk -F\, '{print $2}'`" |bc -l) + CCD_3T4_MOV_X=`echo "${DFT_CCD_EDGECUT}*${DX_3T4}/${CCD_WIDTH}" | bc -l` + CCD_3T4_MOV_Y=`echo "${DFT_CCD_EDGECUT}*${DY_3T4}/${CCD_WIDTH}" | bc -l` + CCD_4_X=$(echo "` echo ${CCD_4_RAW} |awk -F\, '{print $1}' ` + `echo ${CCD_2T4_MOV_X}` +`echo ${CCD_3T4_MOV_X}` "| bc -l) + CCD_4_Y=$(echo "` echo ${CCD_4_RAW} |awk -F\, '{print $2}' ` + `echo ${CCD_2T4_MOV_Y}` +`echo ${CCD_3T4_MOV_Y}` "| bc -l) + REG_CCD_CUT=`echo "polygon(${CCD_1_X},${CCD_1_Y},${CCD_2_X},${CCD_2_Y},${CCD_4_X},${CCD_4_Y},${CCD_3_X},${CCD_3_Y})"` + echo ${REG_CCD_CUT} >>${TMP_REG_FILE_CCD} +done + REG_FILE_CCD="_ccd.reg" + [ -e "${REG_FILE_CCD}" ] && mv -fv ${REG_FILE_CCD} ${REG_FILE_CCD}_bak + # echo "` cat ${TMP_REG_FILE_CCD} | head -n 1 | tail -n 1` + ` cat ${TMP_REG_FILE_CCD} | head -n 2 | tail -n 1` +`cat ${TMP_REG_FILE_CCD} | head -n 3 | tail -n 1`+`cat ${TMP_REG_FILE_CCD} | head -n 4 | tail -n 1 `" >${REG_FILE_CCD} + cat "${TMP_REG_FILE_CCD}">${REG_FILE_CCD} +else + # + printf "*** ERROR ACIS_TYPE ***\n" + exit 255 +fi +## }}} + +## cut ccd region edge +echo "${REG_CCD_RAW}" >_ccd_raw.reg + +# exit 233 + +## generate new regions within CCD for dmextract +SBP_REG_INCCD="_${SBP_REG%.reg}_inccd.reg" +[ -e "${SBP_REG_INCCD}" ] && mv -fv ${SBP_REG_INCCD} ${SBP_REG_INCCD}_bak + +echo "CMD: cat ${CELL_REG} | grep \( | sed -e ':a;N;s/\n/-/;ta'" +CELL_REG_USE=`cat ${CELL_REG} | grep \( | sed -e ':a;N;s/\n/-/;ta'` +# exit 233 + +if [ "${ACIS_TYPE}" = "S" ]; then + grep -iE '^(pie|annulus)' ${SBP_REG_FIX} | sed "s/$/\ \&\ `cat ${REG_FILE_CCD}`/" | sed "s/$/\ \-\ ${CELL_REG_USE}/" > ${SBP_REG_INCCD} + else + L=`cat ${SBP_REG_FIX} | wc -l ` + + for i in `seq 1 $L` ; do + echo "`cat ${SBP_REG_FIX} |head -n $i | tail -n 1 ` & `cat ${REG_FILE_CCD} | head -n 1 `- ${CELL_REG_USE} | `cat ${SBP_REG_FIX} |head -n $i | tail -n 1` & `cat ${REG_FILE_CCD} | head -n 2| tail -n 1 `- ${CELL_REG_USE} |`cat ${SBP_REG_FIX} |head -n $i | tail -n 1 ` & `cat ${REG_FILE_CCD} | head -n 3 | tail -n 1 `- ${CELL_REG_USE} |`cat ${SBP_REG_FIX} |head -n $i | tail -n 1 ` & `cat ${REG_FILE_CCD} | tail -n 1 `- ${CELL_REG_USE} " >>${SBP_REG_INCCD} + done +fi +# ds9 ${EVT_E} -region ${SBP_REG_INCCD} + +## `surface brightness profile' related data {{{ +## extract sbp +printf "extract surface brightness profile ...\n" +SBP_DAT="${SBP_REG%.reg}.fits" +[ -e "${SBP_DAT}" ] && mv -fv ${SBP_DAT} ${SBP_DAT}_bak +if [ -r "${BKG}" ]; then + EXPO_EVT=`dmkeypar ${EVT_E} EXPOSURE echo=yes` + EXPO_BKG=`dmkeypar ${BKG} EXPOSURE echo=yes` + BKG_NORM=`echo "${EXPO_EVT} ${EXPO_BKG}" | awk '{ printf("%g", $1/$2) }'` + printf " == (BKG subtracted; bkgnorm=${BKG_NORM}, energy:${E_RANGE}) ==\n" + punlearn dmextract + dmextract infile="${EVT_E}[bin sky=@${SBP_REG_INCCD}]" outfile="${SBP_DAT}" \ + exp="${EXPMAP}" bkg="${BKG}[energy=${E_RANGE}][bin sky=@${SBP_REG_INCCD}]" \ + bkgexp=")exp" bkgnorm=${BKG_NORM} opt=generic clobber=yes +else + punlearn dmextract + dmextract infile="${EVT_E}[bin sky=@${SBP_REG_INCCD}]" outfile="${SBP_DAT}" \ + exp="${EXPMAP}" opt=generic clobber=yes +fi + +## add `rmid' column +printf "add \`RMID' & \`R_ERR' column ...\n" +SBP_RMID="${SBP_DAT%.fits}_rmid.fits" +[ -e "${SBP_RMID}" ] && mv -fv ${SBP_RMID} ${SBP_RMID}_bak +punlearn dmtcalc +dmtcalc infile="${SBP_DAT}" outfile="${SBP_RMID}" \ + expression="RMID=(R[0]+R[1])/2,R_ERR=(R[1]-R[0])/2" \ + clobber=yes + +## output needed sbp data to files +printf "output needed sbp data ...\n" +SBP_TXT="${SBP_DAT%.fits}.txt" +SBP_QDP="${SBP_DAT%.fits}.qdp" +[ -e "${SBP_TXT}" ] && mv -fv ${SBP_TXT} ${SBP_TXT}_bak +[ -e "${SBP_QDP}" ] && mv -fv ${SBP_QDP} ${SBP_QDP}_bak +punlearn dmlist +dmlist infile="${SBP_RMID}[cols RMID,R_ERR,SUR_FLUX,SUR_FLUX_ERR]" \ + outfile="${SBP_TXT}" opt="data,clean" + +## QDP for sbp {{{ +printf "generate a handy QDP file for sbp ...\n" +cp -fv ${SBP_TXT} ${SBP_QDP} +# change comment sign +sed -i'' 's/#/!/g' ${SBP_QDP} +# add QDP commands +sed -i'' '1 i\ +READ SERR 1 2' ${SBP_QDP} +sed -i'' '2 i\ +LABEL Y "Surface Flux (photons/cm\\u2\\d/pixel\\u2\\d/s)"' ${SBP_QDP} +sed -i'' '2 i\ +LABEL X "Radius (pixel)"' ${SBP_QDP} +sed -i'' '2 i\ +LABEL T "Surface Brightness Profile"' ${SBP_QDP} +## QDP }}} + +printf "generate sbp fitting needed files ...\n" +SBP_RADIUS="radius_sbp.txt" +SBP_FLUX="flux_sbp.txt" +[ -e "${SBP_RADIUS}" ] && mv -fv ${SBP_RADIUS} ${SBP_RADIUS}_bak +[ -e "${SBP_FLUX}" ] && mv -fv ${SBP_FLUX} ${SBP_FLUX}_bak +punlearn dmlist +dmlist infile="${SBP_RMID}[cols R]" \ + opt="data,clean" | awk '{ print $2 }' > ${SBP_RADIUS} +# change the first line `R[2]' to `0.0' +sed -i'' 's/R.*/0\.0/' ${SBP_RADIUS} +dmlist infile="${SBP_RMID}[cols SUR_FLUX,SUR_FLUX_ERR]" \ + opt="data,clean" > ${SBP_FLUX} +# remove the first comment line +sed -i'' '/#.*/d' ${SBP_FLUX} + +## sbp data }}} + +## main }}} + +exit 0 + diff --git a/scripts/ciao_sbp_v3.1.sh b/scripts/ciao_sbp_v3.1.sh deleted file mode 100755 index d111b7f..0000000 --- a/scripts/ciao_sbp_v3.1.sh +++ /dev/null @@ -1,529 +0,0 @@ -#!/bin/sh -# -unalias -a -export LC_COLLATE=C -########################################################### -## 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 ## -## August 16, 2012 ## -########################################################### - -########################################################### -## Changes made by zzh: (2013-02-01) -## removes the region in ccd gap of ACIS_I -## removes the region in the area of point source -## need asol file to prevent offset -########################################################### - -SCRIPT_PATH=`readlink -f $0` -SCRIPT_DIR=`dirname ${SCRIPT_PATH}` -CCDGAP_SCRIPT="chandra_ccdgap_rect.py" - -## about, used in `usage' {{{ -VERSION="v3.1" -UPDATE="2013-02-01" -## about }}} - -## error code {{{ -ERR_USG=1 -ERR_DIR=11 -ERR_EVT=12 -ERR_BKG=13 -ERR_REG=14 -ERR_CELL_REG=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_e= reg= expmap= cellreg= aspec= [bkg= log= ]\n" - printf "\nversion:\n" - printf "${VERSION}, ${UPDATE}\n" - exit ${ERR_USG} - ;; -esac -## usage, help }}} - -## default parameters {{{ -# default energy band -E_RANGE="700:7000" -# default ccd edge cut pixel (20pixel) -DFT_CCD_EDGECUT=25 -# default `event file' which used to match `blanksky' files -#DFT_EVT="_NOT_EXIST_" -DFT_EVT_E="`ls evt2_*_e*.fits 2> /dev/null`" -# default expmap -DFT_EXPMAP="`ls expmap*.fits 2> /dev/null | head -n 1`" -# default `radial region file' to extract surface brightness -#DFT_SBP_REG="_NOT_EXIST_" -DFT_SBP_REG="sbprofile.reg" -# defalut pointsource region file -DFT_CELL_REG="`ls celld*.reg 2> /dev/null`" -# defalut asol file -DFT_ASOL_FILE="`ls pcad*asol*fits 2> /dev/null`" - -# default `log file' -DFT_LOGFILE="sbp_`date '+%Y%m%d'`.log" - -## 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_e}" ]; then - EVT_E=${evt_e} -elif [ -r "${DFT_EVT_E}" ]; then - EVT_E=${DFT_EVT_E} -else - read -p "clean evt2 file: " EVT_E - if [ ! -r "${EVT_E}" ]; then - printf "ERROR: cannot access given \`${EVT_E}' evt file\n" - exit ${ERR_EVT} - fi -fi -printf "## use evt_eng file: \`${EVT_E}'\n" | ${TOLOG} -# check evt file -if [ -r "${expmap}" ]; then - EXPMAP=${expmap} -elif [ -r "${DFT_EXPMAP}" ]; then - EXPMAP=${DFT_EXPMAP} -else - read -p "expocure map file: " EXPMAP - if [ ! -r "${EXPMAP}" ]; then - printf "ERROR: cannot access given \`${EXPMAP}' expmap file\n" - exit ${ERR_EVT} - fi -fi -printf "## use expmap file: \`${EXPMAP}'\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 bkg -if [ ! -z "${bkg}" ]; then - BKG=${bkg} -#else -# read -p "> bkg: " BKG -fi -if [ -r "${BKG}" ]; then - printf "## use bkg: \`${BKG}'\n" | ${TOLOG} -else - BKG="NULL" -fi -# check cell region file -if [ -r "${cellreg}" ]; then - CELL_REG="${cellreg}" -elif [ -r "${DFT_CELL_REG}" ] ; then - CELL_REG="${DFT_CELL_REG}" -else - read -p "> celldetect region file: " CELL_REG - if [ ! -r "${CELL_REG}" ]; then - printf "ERROR: cannot access given \`${CELL_REG}' region file \n" - exit ${ERR_CELL_REG} - fi -fi -printf "## use cell reg file(s): \`${CELL_REG}'\n" | ${TOLOG} -# check asol file -if [ -r "${aspec}" ]; then - ASOL_FILE="${aspec}" -elif [ -r "${DFT_ASOL_FILE}" ] ; then - ASOL_FILE="${DFT_ASOL_FILE}" -else - read -p ">asol file: " ASOL_FILE - if [ ! -r "${ASOL_FILE}" ] ; then - printf " ERROR: cannot access asol file \n" - exit ${ERR_ASOL} - fi -fi -printf "## use asol file(s) : \`${ASOL_FILE}'\n" | ${TOLOG} -## parameters }}} - -## determine ACIS type {{{ -# consistent with `ciao_procevt' -punlearn dmkeypar -DETNAM=`dmkeypar ${EVT_E} DETNAM echo=yes` -if echo ${DETNAM} | grep -q 'ACIS-0123'; then - printf "## \`DETNAM' (${DETNAM}) has chips 0123 -> ACIS-I\n" - ACIS_TYPE="I" - CCD="0:3" -elif echo ${DETNAM} | grep -q 'ACIS-[0-6]*7'; then - printf "## \`DETNAM' (${DETNAM}) has chip 7 -> ACIS-S\n" - ACIS_TYPE="S" - CCD="7" -else - printf "ERROR: unknown detector type: ${DETNAM}\n" - exit ${ERR_DET} -fi -## ACIS type }}} - -## check validity of pie region {{{ -INVALID=`grep -i 'pie' ${SBP_REG} | awk -F'[,()]' '$7 > 360'` -SBP_REG_FIX="_${SBP_REG%.reg}_fix.reg" -if [ "x${INVALID}" != "x" ]; then - printf "*** WARNING: some pie regions' END_ANGLE > 360\n" | ${TOLOG} - printf "*** script will fix ...\n" - cp -fv ${SBP_REG} ${SBP_REG_FIX} - # using `awk' to fix - awk -F'[,()]' '{ - if ($7 > 360) { - printf "%s(%.2f,%.2f,%.2f,%.2f,%.2f,%.2f)\n", $1,$2,$3,$4,$5,$6,($7-360) - } - else { - print $0 - } - }' ${SBP_REG} | sed '/^\#/d' > ${SBP_REG_FIX} -else - cat ${SBP_REG} | sed '/^\#/d' >${SBP_REG_FIX} -fi -## pie validity }}} - -## main process {{{ - -## 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" -punlearn skyfov -skyfov infile="${EVT_E}" outfile="${SKYFOV}" aspec="${ASOL_FILE}" clobber=yes - -## get CCD fov regions {{{ -printf "make regions for CCD ...\n" -TMP_LIST="_tmp_list.txt" -TMP_REC="_tmp_rec.reg" -if [ "${ACIS_TYPE}" = "S" ]; then - # ACIS-S - punlearn dmlist - dmlist infile="${SKYFOV}[ccd_id=${CCD}][cols POS]" opt="data,clean" | awk '{for (i=1;i<=NF;i++) print $i }' |sed -e ':a;N;s/\n/,/;ta' | awk -F"]," '{print "polygon("$2}' | awk -F"NaN" '{print $1}' >${TMP_LIST} - python ${SCRIPT_DIR}/${CCDGAP_SCRIPT} ${TMP_LIST} >${TMP_REC} - XC=` cat ${TMP_REC} | awk -F\( '{print $2}' |awk -F\) '{print $1}' |awk -F\, '{print $1}'` - YC=` cat ${TMP_REC} | awk -F\( '{print $2}' |awk -F\) '{print $1}' |awk -F\, '{print $2}'` - ADD_L=` cat ${TMP_REC} | awk -F\( '{print $2}' |awk -F\) '{print $1}' |awk -F\, '{print $3/2}'` - ANG=` cat ${TMP_REC} | awk -F\( '{print $2}' |awk -F\) '{print $1}' |awk -F\, '{print $5}'` - while [ 1 -eq 1 ]; do - if [ `echo "${ANG} < 0" |bc -l ` -eq 1 ] ; then - ANG=` echo " ${ANG} + 90 " | bc -l ` - elif [ `echo "${ANG} >=90" |bc -l ` -eq 1 ] ; then - ANG=` echo " ${ANG} - 90 " | bc -l` - else - break - fi - done - ANG=`echo "${ANG}/180*3.1415926" |bc -l` - CCD_1_X_RAW=` echo " ${XC} ${ADD_L} ${ANG} "| awk '{print $1-$2*cos($3)-$2*sin($3)}' ` - CCD_2_X_RAW=` echo " ${XC} ${ADD_L} ${ANG} "| awk '{print $1+$2*cos($3)-$2*sin($3)}' ` - CCD_3_X_RAW=` echo " ${XC} ${ADD_L} ${ANG} "| awk '{print $1-$2*cos($3)+$2*sin($3)}' ` - CCD_4_X_RAW=` echo " ${XC} ${ADD_L} ${ANG} "| awk '{print $1+$2*cos($3)+$2*sin($3)}' ` - CCD_1_Y_RAW=` echo " ${YC} ${ADD_L} ${ANG} "| awk '{print $1+$2*cos($3)-$2*sin($3)}' ` - CCD_2_Y_RAW=` echo " ${YC} ${ADD_L} ${ANG} "| awk '{print $1+$2*cos($3)+$2*sin($3)}' ` - CCD_3_Y_RAW=` echo " ${YC} ${ADD_L} ${ANG} "| awk '{print $1-$2*cos($3)-$2*sin($3)}' ` - CCD_4_Y_RAW=` echo " ${YC} ${ADD_L} ${ANG} "| awk '{print $1-$2*cos($3)+$2*sin($3)}' ` - CCD_1_RAW=` echo "${CCD_1_X_RAW},${CCD_1_Y_RAW}"` - CCD_2_RAW=` echo "${CCD_2_X_RAW},${CCD_2_Y_RAW}"` - CCD_3_RAW=` echo "${CCD_3_X_RAW},${CCD_3_Y_RAW}"` - CCD_4_RAW=` echo "${CCD_4_X_RAW},${CCD_4_Y_RAW}"` - REG_CCD_RAW="`echo "polygon(${CCD_1_RAW}, ${CCD_2_RAW}, ${CCD_4_RAW}, ${CCD_3_RAW}) " `" - DX_2T1=$(echo "`echo ${CCD_2_RAW} | awk -F\, '{print $1}'`-`echo ${CCD_1_RAW} |awk -F\, '{print $1}'`" |bc -l) - DY_2T1=$(echo "`echo ${CCD_2_RAW} | awk -F\, '{print $2}'`-`echo ${CCD_1_RAW} |awk -F\, '{print $2}'`" |bc -l) - CCD_WIDTH=`echo "sqrt(${DX_2T1}*${DX_2T1}+${DY_2T1}*${DY_2T1})" | bc -l` - CCD_2T1_MOV_X=`echo "${DFT_CCD_EDGECUT}*${DX_2T1}/${CCD_WIDTH}" | bc -l` - CCD_2T1_MOV_Y=`echo "${DFT_CCD_EDGECUT}*${DY_2T1}/${CCD_WIDTH}" | bc -l` - DX_3T1=$(echo "`echo ${CCD_3_RAW} | awk -F\, '{print $1}'`-`echo ${CCD_1_RAW} |awk -F\, '{print $1}'`" |bc -l) - DY_3T1=$(echo "`echo ${CCD_3_RAW} | awk -F\, '{print $2}'`-`echo ${CCD_1_RAW} |awk -F\, '{print $2}'`" |bc -l) - CCD_3T1_MOV_X=`echo "${DFT_CCD_EDGECUT}*${DX_3T1}/${CCD_WIDTH}" | bc -l` - CCD_3T1_MOV_Y=`echo "${DFT_CCD_EDGECUT}*${DY_3T1}/${CCD_WIDTH}" | bc -l` - CCD_1_X=$(echo "` echo ${CCD_1_RAW} |awk -F\, '{print $1}' ` + `echo ${CCD_2T1_MOV_X}` +`echo ${CCD_3T1_MOV_X}` "| bc -l) - CCD_1_Y=$(echo "` echo ${CCD_1_RAW} |awk -F\, '{print $2}' ` + `echo ${CCD_2T1_MOV_Y}` +`echo ${CCD_3T1_MOV_Y}` "| bc -l) - DX_1T2=$(echo "`echo ${CCD_1_RAW} | awk -F\, '{print $1}'`-`echo ${CCD_2_RAW} |awk -F\, '{print $1}'`" |bc -l) - DY_1T2=$(echo "`echo ${CCD_1_RAW} | awk -F\, '{print $2}'`-`echo ${CCD_2_RAW} |awk -F\, '{print $2}'`" |bc -l) - CCD_1T2_MOV_X=`echo "${DFT_CCD_EDGECUT}*${DX_1T2}/${CCD_WIDTH}" | bc -l` - CCD_1T2_MOV_Y=`echo "${DFT_CCD_EDGECUT}*${DY_1T2}/${CCD_WIDTH}" | bc -l` - DX_4T2=$(echo "`echo ${CCD_4_RAW} | awk -F\, '{print $1}'`-`echo ${CCD_2_RAW} |awk -F\, '{print $1}'`" |bc -l) - DY_4T2=$(echo "`echo ${CCD_4_RAW} | awk -F\, '{print $2}'`-`echo ${CCD_2_RAW} |awk -F\, '{print $2}'`" |bc -l) - CCD_4T2_MOV_X=`echo "${DFT_CCD_EDGECUT}*${DX_4T2}/${CCD_WIDTH}" | bc -l` - CCD_4T2_MOV_Y=`echo "${DFT_CCD_EDGECUT}*${DY_4T2}/${CCD_WIDTH}" | bc -l` - CCD_2_X=$(echo "` echo ${CCD_2_RAW} |awk -F\, '{print $1}' ` + `echo ${CCD_1T2_MOV_X}` +`echo ${CCD_4T2_MOV_X}` "| bc -l) - CCD_2_Y=$(echo "` echo ${CCD_2_RAW} |awk -F\, '{print $2}' ` + `echo ${CCD_1T2_MOV_Y}` +`echo ${CCD_4T2_MOV_Y}` "| bc -l) - DX_1T3=$(echo "`echo ${CCD_1_RAW} | awk -F\, '{print $1}'`-`echo ${CCD_3_RAW} |awk -F\, '{print $1}'`" |bc -l) - DY_1T3=$(echo "`echo ${CCD_1_RAW} | awk -F\, '{print $2}'`-`echo ${CCD_3_RAW} |awk -F\, '{print $2}'`" |bc -l) - CCD_1T3_MOV_X=`echo "${DFT_CCD_EDGECUT}*${DX_1T3}/${CCD_WIDTH}" | bc -l` - CCD_1T3_MOV_Y=`echo "${DFT_CCD_EDGECUT}*${DY_1T3}/${CCD_WIDTH}" | bc -l` - DX_4T3=$(echo "`echo ${CCD_4_RAW} | awk -F\, '{print $1}'`-`echo ${CCD_3_RAW} |awk -F\, '{print $1}'`" |bc -l) - DY_4T3=$(echo "`echo ${CCD_4_RAW} | awk -F\, '{print $2}'`-`echo ${CCD_3_RAW} |awk -F\, '{print $2}'`" |bc -l) - CCD_4T3_MOV_X=`echo "${DFT_CCD_EDGECUT}*${DX_4T3}/${CCD_WIDTH}" | bc -l` - CCD_4T3_MOV_Y=`echo "${DFT_CCD_EDGECUT}*${DY_4T3}/${CCD_WIDTH}" | bc -l` - CCD_3_X=$(echo "` echo ${CCD_3_RAW} |awk -F\, '{print $1}' ` + `echo ${CCD_1T3_MOV_X}` +`echo ${CCD_4T3_MOV_X}` "| bc -l) - CCD_3_Y=$(echo "` echo ${CCD_3_RAW} |awk -F\, '{print $2}' ` + `echo ${CCD_1T3_MOV_Y}` +`echo ${CCD_4T3_MOV_Y}` "| bc -l) - DX_2T4=$(echo "`echo ${CCD_2_RAW} | awk -F\, '{print $1}'`-`echo ${CCD_4_RAW} |awk -F\, '{print $1}'`" |bc -l) - DY_2T4=$(echo "`echo ${CCD_2_RAW} | awk -F\, '{print $2}'`-`echo ${CCD_4_RAW} |awk -F\, '{print $2}'`" |bc -l) - CCD_2T4_MOV_X=`echo "${DFT_CCD_EDGECUT}*${DX_2T4}/${CCD_WIDTH}" | bc -l` - CCD_2T4_MOV_Y=`echo "${DFT_CCD_EDGECUT}*${DY_2T4}/${CCD_WIDTH}" | bc -l` - DX_3T4=$(echo "`echo ${CCD_3_RAW} | awk -F\, '{print $1}'`-`echo ${CCD_4_RAW} |awk -F\, '{print $1}'`" |bc -l) - DY_3T4=$(echo "`echo ${CCD_3_RAW} | awk -F\, '{print $2}'`-`echo ${CCD_4_RAW} |awk -F\, '{print $2}'`" |bc -l) - CCD_3T4_MOV_X=`echo "${DFT_CCD_EDGECUT}*${DX_3T4}/${CCD_WIDTH}" | bc -l` - CCD_3T4_MOV_Y=`echo "${DFT_CCD_EDGECUT}*${DY_3T4}/${CCD_WIDTH}" | bc -l` - CCD_4_X=$(echo "` echo ${CCD_4_RAW} |awk -F\, '{print $1}' ` + `echo ${CCD_2T4_MOV_X}` +`echo ${CCD_3T4_MOV_X}` "| bc -l) - CCD_4_Y=$(echo "` echo ${CCD_4_RAW} |awk -F\, '{print $2}' ` + `echo ${CCD_2T4_MOV_Y}` +`echo ${CCD_3T4_MOV_Y}` "| bc -l) - REG_CCD_CUT=`echo "polygon(${CCD_1_X},${CCD_1_Y},${CCD_2_X},${CCD_2_Y},${CCD_4_X},${CCD_4_Y},${CCD_3_X},${CCD_3_Y})"` - REG_FILE_CCD="_ccd.reg" - [ -e "${REG_FILE_CCD}" ] && mv -f ${REG_FILE_CCD} ${REG_FILE_CCD}_bak - echo "${REG_CCD_CUT}" >>${REG_FILE_CCD} - -elif [ "${ACIS_TYPE}" = "I" ]; then - # ACIS-I - TMP_REG_FILE_CCD="_ccd_tmp.reg" - [ -e "${TMP_REG_FILE_CCD}" ] && mv -f ${TMP_REG_FILE_CCD} ${TMP_REG_FILE_CCD}_bak - for i in `seq 0 3` ; do - punlearn dmlist - dmlist infile="${SKYFOV}[ccd_id=${i}][cols POS]" opt="data,clean" | awk '{for (i=1;i<=NF;i++) print $i }' |sed -e ':a;N;s/\n/,/;ta' | awk -F"]," '{print "polygon("$2}' | awk -F"NaN" '{print $1}' >${TMP_LIST} - python ${SCRIPT_DIR}/${CCDGAP_SCRIPT} ${TMP_LIST} >${TMP_REC} - XC=` cat ${TMP_REC} | awk -F\( '{print $2}' |awk -F\) '{print $1}' |awk -F\, '{print $1}'` - YC=` cat ${TMP_REC} | awk -F\( '{print $2}' |awk -F\) '{print $1}' |awk -F\, '{print $2}'` - ADD_L=` cat ${TMP_REC} | awk -F\( '{print $2}' |awk -F\) '{print $1}' |awk -F\, '{print $3/2}'` - ANG=` cat ${TMP_REC} | awk -F\( '{print $2}' |awk -F\) '{print $1}' |awk -F\, '{print $5}'` - while [ 1 -eq 1 ]; do - if [ `echo "${ANG} < 0" |bc -l ` -eq 1 ] ; then - ANG=` echo " ${ANG} + 90 " | bc -l ` - elif [ `echo "${ANG} >=90" |bc -l ` -eq 1 ] ; then - ANG=` echo " ${ANG} - 90 " | bc -l` - else - break - fi - done - ANG=`echo "${ANG}/180*3.1415926" |bc -l` - CCD_1_X_RAW=` echo " ${XC} ${ADD_L} ${ANG} "| awk '{print $1-$2*cos($3)-$2*sin($3)}' ` - CCD_2_X_RAW=` echo " ${XC} ${ADD_L} ${ANG} "| awk '{print $1+$2*cos($3)-$2*sin($3)}' ` - CCD_3_X_RAW=` echo " ${XC} ${ADD_L} ${ANG} "| awk '{print $1-$2*cos($3)+$2*sin($3)}' ` - CCD_4_X_RAW=` echo " ${XC} ${ADD_L} ${ANG} "| awk '{print $1+$2*cos($3)+$2*sin($3)}' ` - CCD_1_Y_RAW=` echo " ${YC} ${ADD_L} ${ANG} "| awk '{print $1+$2*cos($3)-$2*sin($3)}' ` - CCD_2_Y_RAW=` echo " ${YC} ${ADD_L} ${ANG} "| awk '{print $1+$2*cos($3)+$2*sin($3)}' ` - CCD_3_Y_RAW=` echo " ${YC} ${ADD_L} ${ANG} "| awk '{print $1-$2*cos($3)-$2*sin($3)}' ` - CCD_4_Y_RAW=` echo " ${YC} ${ADD_L} ${ANG} "| awk '{print $1-$2*cos($3)+$2*sin($3)}' ` - CCD_1_RAW=` echo "${CCD_1_X_RAW},${CCD_1_Y_RAW}"` - CCD_2_RAW=` echo "${CCD_2_X_RAW},${CCD_2_Y_RAW}"` - CCD_3_RAW=` echo "${CCD_3_X_RAW},${CCD_3_Y_RAW}"` - CCD_4_RAW=` echo "${CCD_4_X_RAW},${CCD_4_Y_RAW}"` - REG_CCD_RAW=`echo "polygon(${CCD_1_RAW}, ${CCD_2_RAW}, ${CCD_4_RAW}, ${CCD_3_RAW}) " ` - DX_2T1=$(echo "`echo ${CCD_2_RAW} | awk -F\, '{print $1}'`-`echo ${CCD_1_RAW} |awk -F\, '{print $1}'`" |bc -l) - DY_2T1=$(echo "`echo ${CCD_2_RAW} | awk -F\, '{print $2}'`-`echo ${CCD_1_RAW} |awk -F\, '{print $2}'`" |bc -l) - CCD_WIDTH=`echo "sqrt(${DX_2T1}*${DX_2T1}+${DY_2T1}*${DY_2T1})" | bc -l` - CCD_2T1_MOV_X=`echo "${DFT_CCD_EDGECUT}*${DX_2T1}/${CCD_WIDTH}" | bc -l` - CCD_2T1_MOV_Y=`echo "${DFT_CCD_EDGECUT}*${DY_2T1}/${CCD_WIDTH}" | bc -l` - DX_3T1=$(echo "`echo ${CCD_3_RAW} | awk -F\, '{print $1}'`-`echo ${CCD_1_RAW} |awk -F\, '{print $1}'`" |bc -l) - DY_3T1=$(echo "`echo ${CCD_3_RAW} | awk -F\, '{print $2}'`-`echo ${CCD_1_RAW} |awk -F\, '{print $2}'`" |bc -l) - CCD_3T1_MOV_X=`echo "${DFT_CCD_EDGECUT}*${DX_3T1}/${CCD_WIDTH}" | bc -l` - CCD_3T1_MOV_Y=`echo "${DFT_CCD_EDGECUT}*${DY_3T1}/${CCD_WIDTH}" | bc -l` - CCD_1_X=$(echo "` echo ${CCD_1_RAW} |awk -F\, '{print $1}' ` + `echo ${CCD_2T1_MOV_X}` +`echo ${CCD_3T1_MOV_X}` "| bc -l) - CCD_1_Y=$(echo "` echo ${CCD_1_RAW} |awk -F\, '{print $2}' ` + `echo ${CCD_2T1_MOV_Y}` +`echo ${CCD_3T1_MOV_Y}` "| bc -l) - DX_1T2=$(echo "`echo ${CCD_1_RAW} | awk -F\, '{print $1}'`-`echo ${CCD_2_RAW} |awk -F\, '{print $1}'`" |bc -l) - DY_1T2=$(echo "`echo ${CCD_1_RAW} | awk -F\, '{print $2}'`-`echo ${CCD_2_RAW} |awk -F\, '{print $2}'`" |bc -l) - CCD_1T2_MOV_X=`echo "${DFT_CCD_EDGECUT}*${DX_1T2}/${CCD_WIDTH}" | bc -l` - CCD_1T2_MOV_Y=`echo "${DFT_CCD_EDGECUT}*${DY_1T2}/${CCD_WIDTH}" | bc -l` - DX_4T2=$(echo "`echo ${CCD_4_RAW} | awk -F\, '{print $1}'`-`echo ${CCD_2_RAW} |awk -F\, '{print $1}'`" |bc -l) - DY_4T2=$(echo "`echo ${CCD_4_RAW} | awk -F\, '{print $2}'`-`echo ${CCD_2_RAW} |awk -F\, '{print $2}'`" |bc -l) - CCD_4T2_MOV_X=`echo "${DFT_CCD_EDGECUT}*${DX_4T2}/${CCD_WIDTH}" | bc -l` - CCD_4T2_MOV_Y=`echo "${DFT_CCD_EDGECUT}*${DY_4T2}/${CCD_WIDTH}" | bc -l` - CCD_2_X=$(echo "` echo ${CCD_2_RAW} |awk -F\, '{print $1}' ` + `echo ${CCD_1T2_MOV_X}` +`echo ${CCD_4T2_MOV_X}` "| bc -l) - CCD_2_Y=$(echo "` echo ${CCD_2_RAW} |awk -F\, '{print $2}' ` + `echo ${CCD_1T2_MOV_Y}` +`echo ${CCD_4T2_MOV_Y}` "| bc -l) - DX_1T3=$(echo "`echo ${CCD_1_RAW} | awk -F\, '{print $1}'`-`echo ${CCD_3_RAW} |awk -F\, '{print $1}'`" |bc -l) - DY_1T3=$(echo "`echo ${CCD_1_RAW} | awk -F\, '{print $2}'`-`echo ${CCD_3_RAW} |awk -F\, '{print $2}'`" |bc -l) - CCD_1T3_MOV_X=`echo "${DFT_CCD_EDGECUT}*${DX_1T3}/${CCD_WIDTH}" | bc -l` - CCD_1T3_MOV_Y=`echo "${DFT_CCD_EDGECUT}*${DY_1T3}/${CCD_WIDTH}" | bc -l` - DX_4T3=$(echo "`echo ${CCD_4_RAW} | awk -F\, '{print $1}'`-`echo ${CCD_3_RAW} |awk -F\, '{print $1}'`" |bc -l) - DY_4T3=$(echo "`echo ${CCD_4_RAW} | awk -F\, '{print $2}'`-`echo ${CCD_3_RAW} |awk -F\, '{print $2}'`" |bc -l) - CCD_4T3_MOV_X=`echo "${DFT_CCD_EDGECUT}*${DX_4T3}/${CCD_WIDTH}" | bc -l` - CCD_4T3_MOV_Y=`echo "${DFT_CCD_EDGECUT}*${DY_4T3}/${CCD_WIDTH}" | bc -l` - CCD_3_X=$(echo "` echo ${CCD_3_RAW} |awk -F\, '{print $1}' ` + `echo ${CCD_1T3_MOV_X}` +`echo ${CCD_4T3_MOV_X}` "| bc -l) - CCD_3_Y=$(echo "` echo ${CCD_3_RAW} |awk -F\, '{print $2}' ` + `echo ${CCD_1T3_MOV_Y}` +`echo ${CCD_4T3_MOV_Y}` "| bc -l) - DX_2T4=$(echo "`echo ${CCD_2_RAW} | awk -F\, '{print $1}'`-`echo ${CCD_4_RAW} |awk -F\, '{print $1}'`" |bc -l) - DY_2T4=$(echo "`echo ${CCD_2_RAW} | awk -F\, '{print $2}'`-`echo ${CCD_4_RAW} |awk -F\, '{print $2}'`" |bc -l) - CCD_2T4_MOV_X=`echo "${DFT_CCD_EDGECUT}*${DX_2T4}/${CCD_WIDTH}" | bc -l` - CCD_2T4_MOV_Y=`echo "${DFT_CCD_EDGECUT}*${DY_2T4}/${CCD_WIDTH}" | bc -l` - DX_3T4=$(echo "`echo ${CCD_3_RAW} | awk -F\, '{print $1}'`-`echo ${CCD_4_RAW} |awk -F\, '{print $1}'`" |bc -l) - DY_3T4=$(echo "`echo ${CCD_3_RAW} | awk -F\, '{print $2}'`-`echo ${CCD_4_RAW} |awk -F\, '{print $2}'`" |bc -l) - CCD_3T4_MOV_X=`echo "${DFT_CCD_EDGECUT}*${DX_3T4}/${CCD_WIDTH}" | bc -l` - CCD_3T4_MOV_Y=`echo "${DFT_CCD_EDGECUT}*${DY_3T4}/${CCD_WIDTH}" | bc -l` - CCD_4_X=$(echo "` echo ${CCD_4_RAW} |awk -F\, '{print $1}' ` + `echo ${CCD_2T4_MOV_X}` +`echo ${CCD_3T4_MOV_X}` "| bc -l) - CCD_4_Y=$(echo "` echo ${CCD_4_RAW} |awk -F\, '{print $2}' ` + `echo ${CCD_2T4_MOV_Y}` +`echo ${CCD_3T4_MOV_Y}` "| bc -l) - REG_CCD_CUT=`echo "polygon(${CCD_1_X},${CCD_1_Y},${CCD_2_X},${CCD_2_Y},${CCD_4_X},${CCD_4_Y},${CCD_3_X},${CCD_3_Y})"` - echo ${REG_CCD_CUT} >>${TMP_REG_FILE_CCD} -done - REG_FILE_CCD="_ccd.reg" - [ -e "${REG_FILE_CCD}" ] && mv -fv ${REG_FILE_CCD} ${REG_FILE_CCD}_bak - # echo "` cat ${TMP_REG_FILE_CCD} | head -n 1 | tail -n 1` + ` cat ${TMP_REG_FILE_CCD} | head -n 2 | tail -n 1` +`cat ${TMP_REG_FILE_CCD} | head -n 3 | tail -n 1`+`cat ${TMP_REG_FILE_CCD} | head -n 4 | tail -n 1 `" >${REG_FILE_CCD} - cat "${TMP_REG_FILE_CCD}">${REG_FILE_CCD} -else - # - printf "*** ERROR ACIS_TYPE ***\n" - exit 255 -fi -## }}} - -## cut ccd region edge -echo "${REG_CCD_RAW}" >_ccd_raw.reg - -# exit 233 - -## generate new regions within CCD for dmextract -SBP_REG_INCCD="_${SBP_REG%.reg}_inccd.reg" -[ -e "${SBP_REG_INCCD}" ] && mv -fv ${SBP_REG_INCCD} ${SBP_REG_INCCD}_bak - -echo "CMD: cat ${CELL_REG} | grep \( | sed -e ':a;N;s/\n/-/;ta'" -CELL_REG_USE=`cat ${CELL_REG} | grep \( | sed -e ':a;N;s/\n/-/;ta'` -# exit 233 - -if [ "${ACIS_TYPE}" = "S" ]; then - grep -iE '^(pie|annulus)' ${SBP_REG_FIX} | sed "s/$/\ \&\ `cat ${REG_FILE_CCD}`/" | sed "s/$/\ \-\ ${CELL_REG_USE}/" > ${SBP_REG_INCCD} - else - L=`cat ${SBP_REG_FIX} | wc -l ` - - for i in `seq 1 $L` ; do - echo "`cat ${SBP_REG_FIX} |head -n $i | tail -n 1 ` & `cat ${REG_FILE_CCD} | head -n 1 `- ${CELL_REG_USE} | `cat ${SBP_REG_FIX} |head -n $i | tail -n 1` & `cat ${REG_FILE_CCD} | head -n 2| tail -n 1 `- ${CELL_REG_USE} |`cat ${SBP_REG_FIX} |head -n $i | tail -n 1 ` & `cat ${REG_FILE_CCD} | head -n 3 | tail -n 1 `- ${CELL_REG_USE} |`cat ${SBP_REG_FIX} |head -n $i | tail -n 1 ` & `cat ${REG_FILE_CCD} | tail -n 1 `- ${CELL_REG_USE} " >>${SBP_REG_INCCD} - done -fi -# ds9 ${EVT_E} -region ${SBP_REG_INCCD} - -## `surface brightness profile' related data {{{ -## extract sbp -printf "extract surface brightness profile ...\n" -SBP_DAT="${SBP_REG%.reg}.fits" -[ -e "${SBP_DAT}" ] && mv -fv ${SBP_DAT} ${SBP_DAT}_bak -if [ -r "${BKG}" ]; then - EXPO_EVT=`dmkeypar ${EVT_E} EXPOSURE echo=yes` - EXPO_BKG=`dmkeypar ${BKG} EXPOSURE echo=yes` - BKG_NORM=`echo "${EXPO_EVT} ${EXPO_BKG}" | awk '{ printf("%g", $1/$2) }'` - printf " == (BKG subtracted; bkgnorm=${BKG_NORM}, energy:${E_RANGE}) ==\n" - punlearn dmextract - dmextract infile="${EVT_E}[bin sky=@${SBP_REG_INCCD}]" outfile="${SBP_DAT}" \ - exp="${EXPMAP}" bkg="${BKG}[energy=${E_RANGE}][bin sky=@${SBP_REG_INCCD}]" \ - bkgexp=")exp" bkgnorm=${BKG_NORM} opt=generic clobber=yes -else - punlearn dmextract - dmextract infile="${EVT_E}[bin sky=@${SBP_REG_INCCD}]" outfile="${SBP_DAT}" \ - exp="${EXPMAP}" opt=generic clobber=yes -fi - -## add `rmid' column -printf "add \`RMID' & \`R_ERR' column ...\n" -SBP_RMID="${SBP_DAT%.fits}_rmid.fits" -[ -e "${SBP_RMID}" ] && mv -fv ${SBP_RMID} ${SBP_RMID}_bak -punlearn dmtcalc -dmtcalc infile="${SBP_DAT}" outfile="${SBP_RMID}" \ - expression="RMID=(R[0]+R[1])/2,R_ERR=(R[1]-R[0])/2" \ - clobber=yes - -## output needed sbp data to files -printf "output needed sbp data ...\n" -SBP_TXT="${SBP_DAT%.fits}.txt" -SBP_QDP="${SBP_DAT%.fits}.qdp" -[ -e "${SBP_TXT}" ] && mv -fv ${SBP_TXT} ${SBP_TXT}_bak -[ -e "${SBP_QDP}" ] && mv -fv ${SBP_QDP} ${SBP_QDP}_bak -punlearn dmlist -dmlist infile="${SBP_RMID}[cols RMID,R_ERR,SUR_FLUX,SUR_FLUX_ERR]" \ - outfile="${SBP_TXT}" opt="data,clean" - -## QDP for sbp {{{ -printf "generate a handy QDP file for sbp ...\n" -cp -fv ${SBP_TXT} ${SBP_QDP} -# change comment sign -sed -i'' 's/#/!/g' ${SBP_QDP} -# add QDP commands -sed -i'' '1 i\ -READ SERR 1 2' ${SBP_QDP} -sed -i'' '2 i\ -LABEL Y "Surface Flux (photons/cm\\u2\\d/pixel\\u2\\d/s)"' ${SBP_QDP} -sed -i'' '2 i\ -LABEL X "Radius (pixel)"' ${SBP_QDP} -sed -i'' '2 i\ -LABEL T "Surface Brightness Profile"' ${SBP_QDP} -## QDP }}} - -printf "generate sbp fitting needed files ...\n" -SBP_RADIUS="radius_sbp.txt" -SBP_FLUX="flux_sbp.txt" -[ -e "${SBP_RADIUS}" ] && mv -fv ${SBP_RADIUS} ${SBP_RADIUS}_bak -[ -e "${SBP_FLUX}" ] && mv -fv ${SBP_FLUX} ${SBP_FLUX}_bak -punlearn dmlist -dmlist infile="${SBP_RMID}[cols R]" \ - opt="data,clean" | awk '{ print $2 }' > ${SBP_RADIUS} -# change the first line `R[2]' to `0.0' -sed -i'' 's/R.*/0\.0/' ${SBP_RADIUS} -dmlist infile="${SBP_RMID}[cols SUR_FLUX,SUR_FLUX_ERR]" \ - opt="data,clean" > ${SBP_FLUX} -# remove the first comment line -sed -i'' '/#.*/d' ${SBP_FLUX} - -## sbp data }}} - -## main }}} - -exit 0 - -- cgit v1.2.2