diff options
Diffstat (limited to 'scripts/ciao_deproj_spectra.sh')
| -rwxr-xr-x | scripts/ciao_deproj_spectra.sh | 172 |
1 files changed, 20 insertions, 152 deletions
diff --git a/scripts/ciao_deproj_spectra.sh b/scripts/ciao_deproj_spectra.sh index 1ec5ed6..e134fde 100755 --- a/scripts/ciao_deproj_spectra.sh +++ b/scripts/ciao_deproj_spectra.sh @@ -35,12 +35,14 @@ ## [4] CIAO 4.6 Release Notes ## http://cxc.cfa.harvard.edu/ciao/releasenotes/ciao_4.6_release.html ## -AUTHOR="Weitian LI <liweitianux@gmail.com>" -CREATED="2012/07/24" -UPDATED="2015/06/03" -VERSION="v10.0" +## Author: Weitian LI <liweitianux@gmail.com> +## Created: 2012-07-24 ## -## ChangeLogs: +## Change logs: +## 2017-02-26, Weitian LI +## * Use 'manifest.py', 'results.py', and 'renorm_spectrum.py' +## * Simplify 'specextract' parameters handling +## * Other cleanups ## v10.0, 2015/06/03, Aaron LI ## * Copy needed pfiles to current working directory, and ## set environment variable $PFILES to use these first. @@ -81,19 +83,14 @@ VERSION="v10.0" ## * added `background renormalization' ## -unalias -a -export LC_COLLATE=C - ## usage, help {{{ case "$1" in -[hH]*|--[hH]*) printf "Usage:\n" - printf " `basename $0` evt=<evt2_clean> reg=<radial_reg> bkgd=<blank_evt | lbkg_reg | bkg_spec> basedir=<base_dir> nh=<nH> z=<redshift> [ grouptype=<NUM_CTS|BIN> grouptypeval=<number> binspec=<binspec> log=<log_file> ]\n" + printf " `basename $0` reg=<radial_reg> [ bkgd=<blank_evt | lbkg_reg | bkg_spec> evt=<evt2_clean> nh=<nH> z=<redshift> [ grouptype=<NUM_CTS|BIN> grouptypeval=<number> binspec=<binspec> log=<log_file> ]\n" printf "\nNotes:\n" printf " If grouptype=NUM_CTS, then grouptypeval required.\n" printf " If grouptype=BIN, then binspec required.\n" - printf "\nVersion:\n" - printf " ${VERSION}, ${UPDATED}\n" exit ${ERR_USG} ;; esac @@ -101,14 +98,9 @@ esac ## default parameters {{{ # default `event file' which used to match `blanksky' files -#DFT_EVT="_NOT_EXIST_" -DFT_EVT="`\ls evt2*_clean.fits`" +DFT_EVT=$(manifest.py -b getpath -r evt2_clean) # 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 parameters for 'dmgroup' DFT_GROUPTYPE="NUM_CTS" DFT_GROUPTYPEVAL="20" @@ -121,18 +113,9 @@ DFT_LOGFILE="deproj_spectra_`date '+%Y%m%d'`.log" 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" +ASOL=$(manifest.py -b -s "," getpath -r asol) +BPIX=$(manifest.py -b getpath -r bpix) +MSK=$(manifest.py -b getpath -r msk) ## default parameters }}} ## error code {{{ @@ -164,39 +147,6 @@ getopt_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 {{{ @@ -212,8 +162,6 @@ 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 @@ -222,12 +170,6 @@ 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) @@ -245,7 +187,7 @@ 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 + BKGD=$(manifest.py -b getpath -r bkg_spec) else BKGD=${bkgd} fi @@ -309,33 +251,18 @@ fi # bkgd }}} # check given nH if [ -z "${nh}" ]; then - read -p "> value of nH: " N_H + N_H=$(results.py -b get nh) 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 + REDSHIFT=$(results.py -b get z) 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 dmgroup parameters: grouptype, grouptypeval, binspec if [ -z "${grouptype}" ]; then GROUPTYPE="${DFT_GROUPTYPE}" @@ -371,37 +298,6 @@ 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 }}} ## prepare parameter files (pfiles) {{{ @@ -490,39 +386,11 @@ for i in `seq ${LINES}`; do # NO background response files # NO background spectrum (generate by self) # NO spectrum grouping (group by self using `dmgroup') - # Determine parameters for different versions of specextract {{{ - # 'pbkfile' parameter deprecated in CIAO-4.6 - if `pget specextract pbkfile >/dev/null 2>&1`; then - P_PBKFILE="pbkfile=${PBK}" - else - P_PBKFILE="" - fi - # specextract: revision 2013-06: - # 'correct' parameter renamed to 'correctpsf' - if `pget specextract correct >/dev/null 2>&1`; then - P_CORRECT="correct=no" - else - P_CORRECT="correctpsf=no" - fi - # specextract: revision 2013-12: - # 'weight' parameter controls whether ONLY ARFs are weighted. - # 'weight_rmf' added to control whether RMFs are weighted. - # NOTE: - # (1) only when 'weight=yes' will the 'weight_rmf' parameter be used. - # (2) no longer distingush between unweighted and weighted reponses - # in file extension; only .arf & .rmf are now used. - if `pget specextract weight_rmf >/dev/null 2>&1`; then - P_WEIGHTRMF="weight_rmf=yes" - else - P_WEIGHTRMF="" - fi - # }}} punlearn specextract specextract infile="${EVT}[sky=region(${REG_CIAO})]" \ - outroot="r${i}_${ROOTNAME}" bkgfile="" asp="@${ASOLIS}" \ + outroot="r${i}_${ROOTNAME}" bkgfile="" asp="${ASOL}" \ mskfile="${MSK}" badpixfile="${BPIX}" \ - ${P_PBKFILE} ${P_CORRECT} \ - weight=yes ${P_WEIGHTRMF} \ + correctpsf=no weight=yes weight_rmf=yes \ energy="0.3:11.0:0.01" channel="1:1024:1" \ bkgresp=no combine=no binarfwmap=2 \ grouptype=NONE binspec=NONE \ @@ -547,7 +415,7 @@ for i in `seq ${LINES}`; do ## bkg renormalization {{{ printf "Renormalize background ...\n" - bkg_renorm ${RSPEC_PI} ${RSPEC_BKG_PI} + renorm_spectrum.py -r ${RSPEC_PI} ${RSPEC_BKG_PI} ## bkg renorm }}} ## group spectrum {{{ @@ -663,7 +531,7 @@ ignore **:0.0-0.7,7.0-** method leven 1000 0.01 # change abundance standard -abund ${ABUND} +abund grsa xsect bcmc cosmo 70 0 0.73 xset delta 0.01 @@ -741,7 +609,7 @@ ignore **:0.0-0.7,7.0-** method leven 1000 0.01 # change abundance standard -abund ${ABUND} +abund grsa xsect bcmc cosmo 70 0 0.73 xset delta 0.01 |