astro/xmm: Add several previously developed scripts

7 files changed, 1019 insertions, 0 deletions

diff --git a/astro/xmm/epic_repro.sh b/astro/xmm/epic_repro.sh
new file mode 100755
index 0000000..9dfce84
--- /dev/null
+++ b/astro/xmm/epic_repro.sh
@@ -0,0 +1,83 @@
+#!/bin/sh
+#
+##################################################################
+## reprocess data                                               ##
+## apply the most recent calibrations                           ##
+##                                                              ##
+## LIweitiaNux, February 17, 2012                               ##
+##################################################################
+
+## init SAS if needed
+if ! echo $PATH | tr ':' '\n' | grep 'xmmsas' > /dev/null; then
+    sasinit > /dev/null     # 'sasinit' is a func defined in ~/.bashrc
+fi
+
+## set needed variables
+WD=`pwd -P`
+ODF=${WD}/ODF
+if ! [ -d ${ODF} ]; then
+    printf "ERROR: NOT found 'ODF' directory in the current place\n"
+    printf "    You may in wrong directory\n"
+    exit 1
+fi
+export SAS_ODF=${ODF}
+export SAS_CCF=${ODF}/ccf.cif
+PPS=${WD}/PPS       # dir to save processed files
+PROC=${WD}/PROC     # processing dir
+[ -d ${PPS} ] && rm -f ${PPS}/* || mkdir ${PPS}
+[ -d ${PROC} ] && rm -f ${PROC}/* || mkdir ${PROC}
+printf "ODF=${ODF}\n"
+printf "PPS=${PPS}\n"
+printf "PROC=${PROC}\n"
+
+## prepare the data
+## Ref: http://heasarc.gsfc.nasa.gov/docs/xmm/abc/node7.html
+cd ${ODF}
+[ -e ./ccf.cif ] && rm -f ./ccf.cif
+cifbuild        # generate new 'ccf.cif'
+if ls *SUM.SAS; then
+    rm -f *SUM.SAS
+fi
+odfingest       # generate new '*SUM.SAS'
+export SAS_ODF=${ODF}/`ls *SUM.SAS`
+
+## rerunning the pipeline
+## Ref: http://heasarc.gsfc.nasa.gov/docs/xmm/abc/node8.html
+cd ${PPS}
+printf "enter directory ${PPS}\n"
+printf "rerunning the pipeline ...\n"
+printf "run 'emproc' to processe MOS* data ...\n"
+emproc      # process 'MOS*' data
+printf "run 'epproc' to processe PN data ...\n"
+epproc      # process 'PN' data
+
+OBS_ID=`ls *MOS1*ImagingEvts* | sed 's/_EMOS.*//' | sed 's/[0-9]*_//'`
+MOS1_EVT=`ls *EMOS1*ImagingEvts*`
+MOS2_EVT=`ls *EMOS2*ImagingEvts*`
+PN_EVT=`ls *PN*ImagingEvts*`
+ATTHK=`ls *AttHk*`
+
+## store variables for later use
+VARS_FILE=${WD}/variables
+printf "save some variables to file '${VARS_FILE}' for later use ...\n"
+[ -e ${VARS_FILE} ] && mv -f ${VARS_FILE} ${VARS_FILE}.bak
+touch ${VARS_FILE}
+printf "#!/bin/sh\n\n" >> ${VARS_FILE}
+printf "# `date`\n\n" >> ${VARS_FILE}
+printf "export SAS_CCFPATH=${SAS_CCFPATH}\n" >> ${VARS_FILE}
+printf "export SAS_CCF=${SAS_CCF}\n" >> ${VARS_FILE}
+printf "export SAS_ODF=${SAS_ODF}\n\n" >> ${VARS_FILE}
+printf "WORK_DIR=${WD}\n" >> ${VARS_FILE}
+printf "ODF=${ODF}\n" >> ${VARS_FILE}
+printf "PPS=${PPS}\n" >> ${VARS_FILE}
+printf "PROC=${PROC}\n\n" >> ${VARS_FILE}
+printf "OBS_ID=${OBS_ID}\n" >> ${VARS_FILE}
+printf "MOS1_EVT=${MOS1_EVT}\n" >> ${VARS_FILE}
+printf "MOS2_EVT=${MOS2_EVT}\n" >> ${VARS_FILE}
+printf "PN_EVT=${PN_EVT}\n" >> ${VARS_FILE}
+printf "ATTHK=${ATTHK}\n\n" >> ${VARS_FILE}
+printf "DONE\n"
+
+printf "!!!ALL FINISHED!!!
+exit 0
+
diff --git a/astro/xmm/make_regfits.sh b/astro/xmm/make_regfits.sh
new file mode 100755
index 0000000..9a38e50
--- /dev/null
+++ b/astro/xmm/make_regfits.sh
@@ -0,0 +1,49 @@
+#/bin/sh
+#
+# Make a new region FITS file based on the original bkg-region, which
+# is created by SAS tool 'region', and custom ds9 region list.
+# The newly created region FITS file can be used to create a new mask.
+#
+# Aaron LI
+# Created: 2015-11-06
+# Updated: 2015-11-09
+#
+# ChangeLog:
+#   * Update sed match pattern to allow '-' minus sign
+#
+
+if [ $# -ne 3 ]; then
+    echo "Usage:"
+    echo "    `basename $0` <orig_reg.fits> <ds9.reg> <new_reg.fits>"
+    echo ""
+    echo "Note: only CIRCLE region are supported!"
+    exit 1
+fi
+
+ORIG_REGFITS="$1"
+DS9REG="$2"
+NEW_REGFITS="$3"
+[ ! -e "${ORIG_REGFITS}" ] && \
+    echo "ERROR: ${ORIG_REGFITS} does not exist" && exit 11
+[ ! -e "${DS9REG}" ] && echo "ERROR: ${DS9REG} does not exist" && exit 12
+[ -e "${NEW_REGFITS}" ] && mv -fv ${NEW_REGFITS} ${NEW_REGFITS}_bak
+
+
+TMP_CDFILE="cdfile.tmp$$"
+TMP_DATAFILE="datafile.tmp$$"
+
+flcol "${ORIG_REGFITS}[REGION]" > ${TMP_CDFILE}
+sed -i'' '/^___Column_Names__/d' ${TMP_CDFILE}
+
+grep -i '^circle' ${DS9REG} | tr 'a-z' 'A-Z' | \
+    sed -e 's/\(CIRCLE\)(\([0-9.-]\+\),\([0-9.-]\+\),\([0-9.-]\+\))/"!\1"  \2 0 0 0  \3 0 0 0  \4 0 0 0  0 0 0 0  1/' | \
+    tr '"' "'" > ${TMP_DATAFILE}
+
+fcreate cdfile=${TMP_CDFILE} datafile=${TMP_DATAFILE} \
+    outfile=${NEW_REGFITS} tbltype="binary" extname="REGION"
+
+cphead "${ORIG_REGFITS}[Primary]" "${NEW_REGFITS}[Primary]"
+cphead "${ORIG_REGFITS}[REGION]" "${NEW_REGFITS}[REGION]"
+
+rm -f ${TMP_CDFILE} ${TMP_DATAFILE}
+
diff --git a/astro/xmm/make_spectrum.sh b/astro/xmm/make_spectrum.sh
new file mode 100755
index 0000000..762f971
--- /dev/null
+++ b/astro/xmm/make_spectrum.sh
@@ -0,0 +1,95 @@
+#!/bin/sh
+#
+# Generate object and background spectra, and generate also RMF and ARF.
+#
+#
+# Weitian LI
+# Created: 2015-11-09
+# Updated: 2015-11-09
+#
+
+
+case "$1" in
+    -[hH]*|--[hH]*)
+        printf "usage:\n"
+        printf "    `basename $0` det=<detector> regtxt=<regtxt> rootname=<rootname> grpcnts=<group_mincounts>\n"
+        exit 1
+        ;;
+esac
+
+## 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 }}}
+
+getopt_keyval "$@"
+
+# Target
+DET="${det}"
+REGTXT="${regtxt}"
+ROOTNAME="${rootname}"
+GRPCNTS="${grpcnts:-50}"
+echo "GRPCNTS: ${GRPCNTS}"
+
+DET_TYPE=`echo "${DET}" | tr -d '12'`
+PREFIX=`\ls ${DET}*-clean.fits | sed 's/^\(mos\|pn\)\([12]*S[0-9]\{3\}\).*$/\2/'`
+[ ! -e "${REGTXT}" ] && echo "ERROR: ${REGTXT} not exist!" && exit 11
+
+# clean-up previously generated files
+#rm -fv ${DET_TYPE}${PREFIX}-corn.fits
+rm -fv ${DET_TYPE}${PREFIX}-obj-im.fits
+rm -fv ${DET_TYPE}${PREFIX}-exp-im.fits ${DET_TYPE}${PREFIX}-obj-im-sp-det.fits
+rm -fv ${DET_TYPE}${PREFIX}-mask-im.fits
+rm -fv ${DET_TYPE}${PREFIX}-obj.pi ${DET_TYPE}${PREFIX}-back.pi
+rm -fv ${DET_TYPE}${PREFIX}.rmf ${DET_TYPE}${PREFIX}.arf
+
+if [ "${DET}" = "pn" ]; then
+    rm -fv ${DET_TYPE}${PREFIX}-obj-oot.pi ${DET_TYPE}${PREFIX}-obj-im-oot.fits
+    #rm -fv ${DET_TYPE}${PREFIX}-corn-oot.pi
+    SPEC_CMD="pn-spectra prefix=${PREFIX} region=${REGTXT} caldb=${SAS_ESAS_CALDB} mask=1 elow=0 ehigh=0 quad1=1 quad2=1 quad3=1 quad4=1"
+    BACK_CMD="pn_back prefix=${PREFIX} caldb=${SAS_ESAS_CALDB} elow=0 ehigh=0 quad1=1 quad2=1 quad3=1 quad4=1"
+else
+    SPEC_CMD="mos-spectra prefix=${PREFIX} region=${REGTXT} caldb=${SAS_ESAS_CALDB} mask=1 elow=0 ehigh=0 ccd1=1 ccd2=1 ccd3=1 ccd4=1 ccd5=1 ccd6=1 ccd7=1"
+    BACK_CMD="mos_back prefix=${PREFIX} caldb=${SAS_ESAS_CALDB} elow=0 ehigh=0 ccd1=1 ccd2=1 ccd3=1 ccd4=1 ccd5=1 ccd6=1 ccd7=1"
+fi
+
+eval ${SPEC_CMD}
+eval ${BACK_CMD}
+
+# Rename products
+mv -v ${DET_TYPE}${PREFIX}-obj.pi             ${DET}_${ROOTNAME}.pi
+mv -v ${DET_TYPE}${PREFIX}-back.pi            ${DET}_${ROOTNAME}_back.pi
+mv -v ${DET_TYPE}${PREFIX}.rmf                ${DET}_${ROOTNAME}.rmf
+mv -v ${DET_TYPE}${PREFIX}.arf                ${DET}_${ROOTNAME}.arf
+mv -v ${DET_TYPE}${PREFIX}-obj-im-sp-det.fits ${DET}_${ROOTNAME}-sp.fits
+mv -v ${DET_TYPE}${PREFIX}-spec.qdp           ${DET}_${ROOTNAME}-spec.qdp
+mv -v ${DET_TYPE}${PREFIX}-aug.qdp            ${DET}_${ROOTNAME}-aug.qdp
+if [ "${DET}" = "pn" ]; then
+    mv -v ${DET_TYPE}${PREFIX}-obj-oot.pi ${DET}_${ROOTNAME}-oot.pi
+    mv -v ${DET_TYPE}${PREFIX}-obj-os.pi ${DET}_${ROOTNAME}-os.pi
+    mv -v ${DET_TYPE}${PREFIX}-aug-spec.qdp ${DET}_${ROOTNAME}-aug-spec.qdp
+    mv -v ${DET_TYPE}${PREFIX}-aug-rev-hard.qdp ${DET}_${ROOTNAME}-aug-rev-hard.qdp
+    mv -v ${DET_TYPE}${PREFIX}-aug-rev-rate.qdp ${DET}_${ROOTNAME}-aug-rev-rate.qdp
+fi
+
+# Group spectrum
+specgroup spectrumset="${DET}_${ROOTNAME}.pi" groupedset="${DET}_${ROOTNAME}_grp.pi" rmfset="${DET}_${ROOTNAME}.rmf" arfset="${DET}_${ROOTNAME}.arf" mincounts=${GRPCNTS}
+
+# Clean up again
+rm -fv ${DET_TYPE}${PREFIX}-?ff.pi  ${DET_TYPE}${PREFIX}-?fc.pi
+rm -fv ${DET_TYPE}${PREFIX}-?obj.pi ${DET_TYPE}${PREFIX}-?oc.pi
+if [ "${DET}" = "pn" ]; then
+    rm -fv ${DET_TYPE}${PREFIX}-?ff-oot.pi  ${DET_TYPE}${PREFIX}-?fc-oot.pi
+    rm -fv ${DET_TYPE}${PREFIX}-?obj-oot.pi ${DET_TYPE}${PREFIX}-?oc-oot.pi
+fi
+
diff --git a/astro/xmm/mos_source_detect.sh b/astro/xmm/mos_source_detect.sh
new file mode 100755
index 0000000..af04e2f
--- /dev/null
+++ b/astro/xmm/mos_source_detect.sh
@@ -0,0 +1,212 @@
+#!/bin/sh
+#
+##################################################################
+## filter, remove source, gti, etc.                             ##
+## Ref: http://heasarc.gsfc.nasa.gov/docs/xmm/abc/node8.html    ##
+##                                                              ##
+## LIweitiaNux, February 17, 2012                               ##
+##################################################################
+
+usage() {
+    printf "Usage:\n"
+    printf "     `basename $0` broad|soft|hard [E_START:E_END]\n"
+}
+
+## process parameters
+if [ $# -eq 1 ]; then
+    E_BAND=$1
+    if [ "${E_BAND}" = "broad" ]; then
+        E_RANGE=300:10000
+    elif [ "${E_BAND}" = "soft" ]; then
+        E_RANGE=300:2000
+    elif [ "${E_BAND}" = "hard" ]; then
+        E_RANGE=2000:10000
+    elif [ "${E_BAND}" = "xid" ]; then
+        E_RANGE=500:4500
+    else
+        printf "ERROR: unknown energy name\n"
+        usage
+        exit 1
+    fi
+elif [ $# -eq 2 ]; then
+    E_BAND=$1
+    E_RANGE=$2
+else
+    printf "ERROR: no parameters given\n"
+    usage
+    exit 2
+fi
+# energy range
+PIMIN=`echo ${E_RANGE} | cut -d':' -f1`
+PIMAX=`echo ${E_RANGE} | cut -d':' -f2`
+if [ "x${PIMIN}" = "x" ] || [ "x${PIMAX}" = "x" ]; then
+    printf "ERROR: given parameters NOT match\n"
+    usage
+    exit 3
+fi
+# get 'ecf'
+read -p "    ECF of MOS1 for this energy band: " ECF_M1
+read -p "    ECF of MOS2 for this energy band: " ECF_M2
+if [ "x${ECF_M1}" = "x" ] || [ "x${ECF_M2}" = "x" ]; then
+    printf "ERROR: ECF given error\n"
+    exit 10
+fi
+
+WD=`pwd -P`
+VARS_FILE=${WD}/variables
+if ! [ -e ${VARS_FILE} ]; then
+    printf "ERROR: NOT found '${VARS_FILE}' in the current directory\n"
+    printf "    You may in wrong directory\n"
+    exit 1
+fi
+## load previous saved variables
+. ${VARS_FILE}
+printf "variables loaded from file '${VARS_FILE}'\n"
+
+## link needed files
+printf "link needed files ...\n"
+cd ${PROC}
+[ -e ${ATTHK} ] || ln -sv ../PPS/${ATTHK} .
+[ -e ${MOS1_EVT} ] || ln -sv ../PPS/${MOS1_EVT} .
+[ -e ${MOS2_EVT} ] || ln -sv ../PPS/${MOS2_EVT} .
+
+## filter the original event file
+printf "filter by energy on the original event file ...\n"
+M1_EVT="${OBS_ID}_mos1.evt"
+M2_EVT="${OBS_ID}_mos2.evt"
+evselect table=${MOS1_EVT} withfilteredset=yes \
+    filteredset=${M1_EVT} filtertype=expression \
+    expression="(PATTERN<=12) && (PI in [200:12000]) && (FLAG==0) && #XMMEA_EM" \
+    keepfilteroutput=yes updateexposure=yes filterexposure=yes
+evselect table=${MOS2_EVT} withfilteredset=yes \
+    filteredset=${M2_EVT} filtertype=expression \
+    expression="(PATTERN<=12) && (PI in [200:12000]) && (FLAG==0) && #XMMEA_EM" \
+    keepfilteroutput=yes updateexposure=yes filterexposure=yes
+
+## make images
+printf "bin event file to make image ...\n"
+M1_IMG="${OBS_ID}_mos1_${E_BAND}.img"
+M2_IMG="${OBS_ID}_mos2_${E_BAND}.img"
+evselect table=${M1_EVT} withimageset=yes imageset=${M1_IMG} \
+    imagebinning=binSize xcolumn=X ximagebinsize=22 \
+    ycolumn=Y yimagebinsize=22 filtertype=expression \
+    expression="(FLAG==0) && (PI in [${E_RANGE}])"
+evselect table=${M2_EVT} withimageset=yes imageset=${M2_IMG} \
+    imagebinning=binSize xcolumn=X ximagebinsize=22 \
+    ycolumn=Y yimagebinsize=22 filtertype=expression \
+    expression="(FLAG==0) && (PI in [${E_RANGE}])"
+
+## make exposure map
+printf "make exposure map ...\n"
+M1_EXP="${OBS_ID}_mos1_${E_BAND}_exp.img"
+M2_EXP="${OBS_ID}_mos2_${E_BAND}_exp.img"
+eexpmap attitudeset=${ATTHK} eventset=${M1_EVT} imageset=${M1_IMG} \
+    expimageset=${M1_EXP} pimin=${PIMIN} pimax=${PIMAX}
+eexpmap attitudeset=${ATTHK} eventset=${M1_EVT} imageset=${M2_IMG} \
+    expimageset=${M2_EXP} pimin=${PIMIN} pimax=${PIMAX}
+
+## make detection mask for detector
+printf "make detection mask for detector ...\n"
+M1_MASK="${OBS_ID}_mos1_${E_BAND}_mask.img"
+M2_MASK="${OBS_ID}_mos2_${E_BAND}_mask.img"
+emask expimageset=${M1_EXP} detmaskset=${M1_MASK}
+emask expimageset=${M2_EXP} detmaskset=${M2_MASK}
+
+## 'local mode' sliding box detection
+## may need to increase parameter 'imagebuffersize', say 'imagebuffersize=2000'
+printf "sliding box detection, in *local* mode ...\n"
+M1_BOXLIST_L="${OBS_ID}_mos1_${E_BAND}_boxlist_local.fits"
+M2_BOXLIST_L="${OBS_ID}_mos2_${E_BAND}_boxlist_local.fits"
+eboxdetect usemap=no likemin=8 withdetmask=yes detmasksets=${M1_MASK} \
+    imagesets=${M1_IMG} expimagesets=${M1_EXP} pimin=${PIMIN} pimax=${PIMAX} \
+    boxlistset=${M1_BOXLIST_L}
+eboxdetect usemap=no likemin=8 withdetmask=yes detmasksets=${M2_MASK} \
+    imagesets=${M2_IMG} expimagesets=${M2_EXP} pimin=${PIMIN} pimax=${PIMAX} \
+    boxlistset=${M2_BOXLIST_L}
+
+## generate spline background map from the non-source regions
+printf "generate background map from the non-source regions ...\n"
+M1_BKG="${OBS_ID}_mos1_${E_BAND}_bkg.img"
+M2_BKG="${OBS_ID}_mos2_${E_BAND}_bkg.img"
+esplinemap bkgimageset=${M1_BKG} imageset=${M1_IMG} scut=0.005 \
+    boxlistset=${M1_BOXLIST_L} nsplinenodes=16 withdetmask=yes \
+    detmaskset=${M1_MASK} withexpimage=yes expimageset=${M1_EXP}
+esplinemap bkgimageset=${M2_BKG} imageset=${M2_IMG} scut=0.005 \
+    boxlistset=${M2_BOXLIST_L} nsplinenodes=16 withdetmask=yes \
+    detmaskset=${M2_MASK} withexpimage=yes expimageset=${M2_EXP}
+
+## sliding box source detection again, in *map* mode
+## improve the detection sensitivity
+## also NOTE the parameter 'imagebuffersize'
+printf "sliding box detection again, in *map* mode ...\n"
+M1_BOXLIST_M="${OBS_ID}_mos1_${E_BAND}_boxlist_map.fits"
+M2_BOXLIST_M="${OBS_ID}_mos2_${E_BAND}_boxlist_map.fits"
+eboxdetect usemap=yes likemin=8 withdetmask=yes detmasksets=${M1_MASK} \
+    imagesets=${M1_IMG} expimagesets=${M1_EXP} pimin=${PIMIN} \
+    pimax=${PIMAX} boxlistset=${M1_BOXLIST_M} \
+    bkgimagesets=${M1_BKG}
+eboxdetect usemap=yes likemin=8 withdetmask=yes detmasksets=${M2_MASK} \
+    imagesets=${M2_IMG} expimagesets=${M2_EXP} pimin=${PIMIN} \
+    pimax=${PIMAX} boxlistset=${M2_BOXLIST_M} \
+    bkgimagesets=${M2_BKG}
+
+## final source list
+## The energy conversion values (ECFs) can be supplied to convert the source
+## count rates into fluxes. The ECFs for each detector and energy band depend
+## on the pattern selection and filter used during the observation. For more
+## information, please consult the calibration paper ``SSC-LUX-TN-0059'',
+## available at the XMM-Newton Science Operations Center or see Table 3.2 in
+## the '2XMM Catalogue User Guide'.
+printf "final source list ...\n"
+M1_EMLLIST="${OBS_ID}_mos1_${E_BAND}_emllist.fits"
+M2_EMLLIST="${OBS_ID}_mos2_${E_BAND}_emllist.fits"
+emldetect imagesets=${M1_IMG} expimagesets=${M1_EXP} bkgimagesets=${M1_BKG} \
+    pimin=${PIMIN} pimax=${PIMAX} boxlistset=${M1_BOXLIST_M} \
+    ecf=${ECF_M1} mlmin=10.0 mllistset=${M1_EMLLIST}
+emldetect imagesets=${M2_IMG} expimagesets=${M2_EXP} bkgimagesets=${M2_BKG} \
+    pimin=${PIMIN} pimax=${PIMAX} boxlistset=${M2_BOXLIST_M} \
+    ecf=${ECF_M2} mlmin=10.0 mllistset=${M2_EMLLIST}
+
+## optional, make a sensitivity map
+printf "optional, make a sensitivity map ...\n"
+M1_SENS="${OBS_ID}_mos1_${E_BAND}_sens.img"
+M2_SENS="${OBS_ID}_mos2_${E_BAND}_sens.img"
+esensmap expimagesets=${M1_EXP} bkgimagesets=${M1_BKG} \
+    detmasksets=${M1_MASK} sensimageset=${M1_SENS} mlmin=10.0
+esensmap expimagesets=${M2_EXP} bkgimagesets=${M2_BKG} \
+    detmasksets=${M2_MASK} sensimageset=${M2_SENS} mlmin=10.0
+
+## write source list to 'ds9 reg' file and view
+M1_SRC="${OBS_ID}_mos1_${E_BAND}_src.reg"
+M2_SRC="${OBS_ID}_mos2_${E_BAND}_src.reg"
+srcdisplay boxlistset=${M1_EMLLIST} imageset=${M1_IMG} \
+    withregionfile=yes regionfile=${M1_SRC} sourceradius=0.003
+srcdisplay boxlistset=${M2_EMLLIST} imageset=${M2_IMG} \
+    withregionfile=yes regionfile=${M2_SRC} sourceradius=0.003
+
+## remove duplicated lines in 'regionfile'
+printf "remove duplicated lines in reg files ...\n"
+mv ${M1_SRC} ${M1_SRC}_tmp
+mv ${M2_SRC} ${M2_SRC}_tmp
+cat ${M1_SRC}_tmp | uniq > ${M1_SRC}
+cat ${M2_SRC}_tmp | uniq > ${M2_SRC}
+rm ${M1_SRC}_tmp ${M2_SRC}_tmp
+
+## save variables
+printf "save variables ...\n"
+CHECK="`grep 'EXP' ${VARS_FILE}`"
+if [ "x${CHECK}" = "x" ]; then
+    printf "# `date`\n" >> ${VARS_FILE}
+    printf "M1_EVT=${M1_EVT}\n" >> ${VARS_FILE}
+    printf "M1_IMG=${M1_IMG}\n" >> ${VARS_FILE}
+    printf "M1_BKG=${M1_BKG}\n" >> ${VARS_FILE}
+    printf "M1_EXP=${M1_EXP}\n" >> ${VARS_FILE}
+    printf "M1_SENS=${M1_SENS}\n" >> ${VARS_FILE}
+    printf "M2_EVT=${M2_EVT}\n" >> ${VARS_FILE}
+    printf "M2_IMG=${M2_IMG}\n" >> ${VARS_FILE}
+    printf "M2_BKG=${M2_BKG}\n" >> ${VARS_FILE}
+    printf "M2_EXP=${M2_EXP}\n" >> ${VARS_FILE}
+    printf "M2_SENS=${M2_SENS}\n" >> ${VARS_FILE}
+    printf "\n" >> ${VARS_FILE}
+fi
+
diff --git a/astro/xmm/sky2det.sh b/astro/xmm/sky2det.sh
new file mode 100755
index 0000000..6a7ca22
--- /dev/null
+++ b/astro/xmm/sky2det.sh
@@ -0,0 +1,53 @@
+#!/bin/sh
+#
+# Convert from SKY(X,Y) coordinate to (DETX,DETY) coordinate for XMM.
+# Using the SAS tool "ecoordconv"
+#
+#
+# Weitian LI
+# Created: 2015-11-09
+# Updated: 2015-11-09
+#
+
+if [ $# -ne 3 ]; then
+    echo "usage:"
+    echo "    `basename $0` <image> <sky.reg> <det.reg>"
+    exit 1
+fi
+
+IMG="$1"
+SKY_REG="$2"
+DET_REG="$3"
+[ -e "${DET_REG}" ] && mv -fv ${DET_REG} ${DET_REG}_bak
+
+sky2det() {
+    _img=$1
+    _x=$2
+    _y=$3
+    _detxy=`ecoordconv imageset=${_img} x=${_x} y=${_y} coordtype=POS -w 0 -V 0 | \
+        grep -E '^\s*DETX:\s*DETY:' | awk '{ print $3, $4 }'`
+    echo "${_detxy}"
+}
+
+
+cat "${SKY_REG}" | while read line; do
+    if echo "${line}" | grep -iq '^circle'; then
+        # Across a circle region
+        echo "CIRCLE: ${line}"
+        CIRCLE_X=`echo "${line}" | awk -F'[(,)]' '{ print $2 }'`
+        CIRCLE_Y=`echo "${line}" | awk -F'[(,)]' '{ print $3 }'`
+        CIRCLE_R=`echo "${line}" | awk -F'[(,)]' '{ print $4 }'`
+        DETXY=`sky2det ${IMG} ${CIRCLE_X} ${CIRCLE_Y}`
+        CIRCLE_DETX=`echo "${DETXY}" | awk '{ print $1 }'`
+        CIRCLE_DETY=`echo "${DETXY}" | awk '{ print $2 }'`
+        CIRCLE_REG="circle(${CIRCLE_DETX},${CIRCLE_DETY},${CIRCLE_R})"
+        echo "${CIRCLE_REG}" >> ${DET_REG}
+    elif echo "${line}" | grep -iq '^physical'; then
+        echo "detector" >> ${DET_REG}
+    else
+        # Just move to output region file
+        echo "${line}" >> ${DET_REG}
+    fi
+done
+
+
diff --git a/astro/xmm/sky2det_manual.sh b/astro/xmm/sky2det_manual.sh
new file mode 100644
index 0000000..589a0b2
--- /dev/null
+++ b/astro/xmm/sky2det_manual.sh
@@ -0,0 +1,80 @@
+#!/bin/sh
+#
+# Convert from SKY(X,Y) coordinate to (DETX,DETY) coordinate for XMM.
+# The required conversion coefficients are extracted from input FITS header.
+#
+#
+# Weitian LI
+# Created: 2015-11-09
+# Updated: 2015-11-09
+#
+
+case "$1" in
+    -[hH]*|--[hH]*)
+        printf "usage:\n"
+        printf "    `basename $0` evt=<evt> x=<sky_x> y=<sky_y>\n"
+        exit 1
+        ;;
+esac
+
+## 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
+}
+
+sky2wcs() {
+    sky_val=$1
+    crpx=$2
+    crvl=$3
+    cdlt=$4
+    wcs_val=`python -c "import math; print(${crvl} + math.tan(${cdlt} * (${sky_val} - (${crpx}))))"`
+    echo ${wcs_val}
+}
+## functions }}}
+
+# process cmdline args using `getopt_keyval'
+getopt_keyval "$@"
+
+[ ! -e "${evt}" ] && echo "ERROR: ${evt} not exist" && exit 11
+
+# Get the WCS conversion coefficients
+fkeypar "${evt}" REFXCTYP
+REFXCTYP=`pget fkeypar value`
+fkeypar "${evt}" REFXCRPX
+REFXCRPX=`pget fkeypar value`
+fkeypar "${evt}" REFXCRVL
+REFXCRVL=`pget fkeypar value`
+fkeypar "${evt}" REFXCDLT
+REFXCDLT=`pget fkeypar value`
+fkeypar "${evt}" REFYCTYP
+REFYCTYP=`pget fkeypar value`
+fkeypar "${evt}" REFYCRPX
+REFYCRPX=`pget fkeypar value`
+fkeypar "${evt}" REFYCRVL
+REFYCRVL=`pget fkeypar value`
+fkeypar "${evt}" REFYCDLT
+REFYCDLT=`pget fkeypar value`
+echo "(X,Y) => (RA,DEC)"
+echo "    RA  = ${REFXCRVL} + TAN[ (${REFXCDLT}) * (X - (${REFXCRPX})) ]"
+echo "    DEC = ${REFYCRVL} + TAN[ (${REFYCDLT}) * (Y - (${REFYCRPX})) ]"
+
+RA=`sky2wcs ${x} ${REFXCRPX} ${REFXCRVL} ${REFXCDLT}`
+DEC=`sky2wcs ${y} ${REFYCRPX} ${REFYCRVL} ${REFYCDLT}`
+echo "sky(${x},${y}) => wcs(${RA},${DEC})"
+
+DETXY=`esky2det datastyle=user ra=${RA} dec=${DEC} outunit=det calinfostyle=set calinfoset="${evt}" | \
+    grep -A 1 -E '^#\s+detX\s+detY' | tail -n 1`
+detx=`echo ${DETXY} | awk '{ print $1 }'`
+dety=`echo ${DETXY} | awk '{ print $2 }'`
+echo "(RA,DEC) => DETXY(${detx},${dety})"
+
diff --git a/astro/xmm/xspec_bkgcorr_xmm.tcl b/astro/xmm/xspec_bkgcorr_xmm.tcl
new file mode 100644
index 0000000..240d4a2
--- /dev/null
+++ b/astro/xmm/xspec_bkgcorr_xmm.tcl
@@ -0,0 +1,447 @@
+#####################################################################
+## XSPEC TCL script
+##
+## To apply background spectrum correction to MOS1 and MOS2 background.
+## It fakes the XB (Galatic + CXB) and NXB (SP + instrument lines)
+## spectra according to the current fitting results, and then invoke
+## "mathpha" to make the final corrected background spectrum.
+##
+## NOTE:
+## The exposure time for the faked spectra is multiplied by a scale
+## factor in order to obtain satisfactory statistics for the higher
+## energy regime.  The same scale factor is also applied to the original
+## FWC background spectrum, considering that many channels have "counts"
+## smaller than one, especially in the higher energy regime.  Otherwise,
+## many of these values will be rounded to *zero* by the "mathpha"
+## operating in the "COUNT" mode.
+## See also the following description to variable "SCALE_FACTOR",
+## procedures "spectrum_rate2counts" and "scale_spectrum".
+##
+##
+## NOTES:
+## Work with XSPEC v12.x
+## Based on "xspec_bkgcorr_v2.tcl" for Chandra background spectrum correction
+##
+## Weitian LI <liweitianux@gmail.com>
+## Created: 2015-11-07
+## Updated: 2015-11-12
+##
+## ChangeLogs:
+## 2015-11-10:
+##   * Remove "diag_rsp"
+##   * Add proc update_header:
+##     Update specified keywords instead of just copy the original header
+##   * Use lists to record the information for all data groups,
+##     therefore this tool can handle any number of data groups.
+## 2015-11-11:
+##   * Add back "diag_rsp"
+## 2015-11-12:
+##   * Update to use "ftkeypar" and "fthedit"
+##   * Add procedure "scale_spectrum" to scale the FWC background
+##   * Add procedure "make_filename"
+##   * Add procedure "spectrum_rate2counts" to convert RATE spectrum to
+##     COUNTS, therefore this script can properly handle PN spectrum
+##   * Update script description
+#####################################################################
+
+
+## Global variables {{{
+set datagrp_total [ tcloutr datagrp ]
+
+# Lists to record the information of each data group
+set rootname       {}
+set rmf            {}
+set arf            {}
+set diag_rsp       {}
+set src_spec       {}
+set bkg_spec       {}
+set src_exposure   {}
+set bkg_exposure   {}
+set src_backscale  {}
+set bkg_backscale  {}
+set modelpars      {}
+set back_modelname {}
+set back_modelcomp {}
+set back_modelpars {}
+# This "modelcomp" is not a list
+set modelcomp      {}
+
+# FITS header keywords needed to be copyed to the corrected spectrum
+set KEYWORDS       {TELESCOP INSTRUME FILTER}
+
+# Scale factor for faking spectra
+# By lengthen the exposure time for the fake process, the total counts
+# of the faked spectra is much larger, therefore the statistical issuses
+# due to small counts (especially in the higher energy regime) is mitigated.
+set SCALE_FACTOR   100
+## Global variables }}}
+
+## Procedures
+# Make new filename by appending the suffix to the stem of filename.
+# e.g., "stem.ext" + "sfix" => "stem_sfix.ext"
+proc make_filename {filename suffix} {
+    regexp -- {(.+)\.([^\.]+)} $filename match stem_ ext_
+    set newname "${stem_}_${suffix}.${ext_}"
+    return $newname
+}
+
+# Determine root name and get other filenames (e.g., rmf, arf, back)
+proc get_filenames {datagrp} {
+    global rootname
+    global src_spec
+    global bkg_spec
+    global rmf
+    global arf
+    global diag_rsp
+    set src_spec_ [ exec basename [ tcloutr filename $datagrp ] ]
+    regexp -- {(.+)_(bin|grp|group)([0-9]*)\.(pi|fits)} $src_spec_ match rootname_ grptype_ grpval_ fext_
+    set bkg_spec_ [ tcloutr backgrnd $datagrp ]
+    set response_ [ tcloutr response $datagrp ]
+    set rmf_      [ lindex $response_ 0 ]
+    set diag_rsp_ [ lindex $response_ 1 ]
+    set arf_      [ lindex [ tcloutr arf      $datagrp ] 0 ]
+    lappend rootname $rootname_
+    lappend src_spec $src_spec_
+    lappend bkg_spec $bkg_spec_
+    lappend rmf      $rmf_
+    lappend arf      $arf_
+    lappend diag_rsp $diag_rsp_
+}
+
+# Get exposure and backscale values
+proc get_exposure_backscale {datagrp} {
+    global src_exposure
+    global bkg_exposure
+    global src_backscale
+    global bkg_backscale
+    set src_exposure_ [ tcloutr expos $datagrp s ]
+    set bkg_exposure_ [ tcloutr expos $datagrp b ]
+    scan [ tcloutr backscal $datagrp s ] "%f" src_backscale_
+    scan [ tcloutr backscal $datagrp b ] "%f" bkg_backscale_
+    lappend src_exposure  $src_exposure_
+    lappend bkg_exposure  $bkg_exposure_
+    lappend src_backscale $src_backscale_
+    lappend bkg_backscale $bkg_backscale_
+}
+
+# Get the model components
+proc get_modelcomp {} {
+    global datagrp_total
+    global modelcomp
+    # NOTE: "model" contains all the model components corresponding to
+    # echo data group, and contains also the background model components.
+    # Therefore, the "model" contents repeats the same model components
+    # for each data group.
+    # e.g.,
+    # if data group 1 has model components: "wabs*apec", and there are
+    # altogether 2 data groups, then the "model" content is:
+    #     "wabs(apec) wabs(apec)"
+    # which is repeated!
+    set model [ tcloutr model ]
+    set indices [ lsearch -all $model {*:} ]
+    if { [ llength $indices ] == 0 } {
+        set modelcomp $model
+    } else {
+        set modelcomp [ lrange $model 0 [ expr {[ lindex $indices 0 ] - 1} ] ]
+    }
+    # Deal with the repeated model components if exist multiple data groups
+    set comp_length [ expr { [ llength $modelcomp ] / $datagrp_total } ]
+    set modelcomp [ lrange $modelcomp 0 [ expr {$comp_length - 1} ] ]
+}
+
+# Get the model parameters
+proc get_modelpars {datagrp} {
+    global datagrp_total
+    global modelpars
+    set modelpars_ ""
+    # NOTE: "modpar" is the total number of parameters,
+    # not just for one data group
+    set npar_ [ expr { [ tcloutr modpar ] / $datagrp_total } ]
+    set par_begin_ [ expr {$npar_ * ($datagrp - 1) + 1} ]
+    set par_end_   [ expr {$npar_ + $par_begin_ - 1} ]
+    for {set i $par_begin_} {$i <= $par_end_} {incr i} {
+        scan [ tcloutr param $i ] "%f" pval_
+        set modelpars_ "$modelpars_ $pval_ &"
+    }
+    lappend modelpars $modelpars_
+}
+
+# Get the model name, components and parameters of the background model
+# corresponding to the given data group
+proc get_backmodel {datagrp} {
+    global back_modelname
+    global back_modelcomp
+    global back_modelpars
+    set model_ [ tcloutr model ]
+    set indices_ [ lsearch -all $model_ {*:} ]
+    set name_idx_ [ lindex $indices_ [ expr {$datagrp - 1} ] ]
+    set back_modelname_ [ regsub {:\s*} [ lindex $model_ $name_idx_ ] "" ]
+    set ncomp_ [ tcloutr modcomp $back_modelname_ ]
+    set back_modelcomp_ [ lrange $model_ [ expr {$name_idx_ + 1} ] [ expr {$name_idx_ + $ncomp_ } ] ]
+    set npar_ [ tcloutr modpar $back_modelname_ ]
+    set back_modelpars_ ""
+    for {set i 1} {$i <= $npar_} {incr i} {
+        scan [ tcloutr param ${back_modelname_}:$i ] "%f" pval_
+        set back_modelpars_ "$back_modelpars_ $pval_ &"
+    }
+    lappend back_modelname $back_modelname_
+    lappend back_modelcomp $back_modelcomp_
+    lappend back_modelpars $back_modelpars_
+}
+
+# Save current XSPEC fitting results
+proc save_xspec_fit {} {
+    set now [ exec date +%Y%m%d%H%M ]
+    set xspec_outfile "xspec_${now}.xcm"
+    if { [ file exists ${xspec_outfile} ] } {
+        exec mv -fv ${xspec_outfile} ${xspec_outfile}_bak
+    }
+    save all "${xspec_outfile}"
+}
+
+# Load model
+proc load_model {comp pars} {
+    model clear
+    model $comp & $pars /*
+}
+
+# Set the norm of the ICM APEC component to be zero,
+# which is the last parameter of current loaded model.
+proc set_icm_norm {} {
+    set npar_ [ tcloutr modpar ]
+    newpar $npar_ 0.0
+    puts "ICM_apec_norm:$npar_: [ tcloutr param $npar_ ]"
+}
+
+# Fake spectrum
+proc fake_spectrum {outfile rmf arf exptime} {
+    if {[ file exists ${outfile} ]} {
+        exec mv -fv ${outfile} ${outfile}_bak
+    }
+    data none
+    fakeit none & $rmf & $arf & y &  & $outfile & $exptime & /*
+}
+
+# Convert a spectrum of "RATE" to a spectrum of "COUNTS" by multiplying
+# the exposure time.  The "STAT_ERR" column is also scaled if exists.
+# NOTE:
+# The MOS1/2 FWC background spectrum generated by "mos_back" has column
+# "COUNTS", however, the PN FWC background spectrum generated by "pn_back"
+# has column "RATE".
+# Therefore, the PN background spectra should be converted to "COUNTS" format,
+# since the other faked CXB & NXB spectra are also in "COUNTS" format.
+# The exposure time is acquired from the "EXPOSURE" keyword in the header
+# of input spectrum.
+#
+# Return:
+#   * 0  : input spectrum is already in "COUNTS" format, no conversion needed
+#   * 1  : input spectrum is in "RATE" format, and is converted to "COUNTS"
+#          format and write to outfile
+#   * -1 : invalid input spectrum
+proc spectrum_rate2counts {outfile infile} {
+    if {[ file exists ${outfile} ]} {
+        exec mv -fv ${outfile} ${outfile}_bak
+    }
+    # Get number of columns/fields
+    exec ftkeypar "${infile}+1" TFIELDS
+    set fields_ [ exec pget ftkeypar value ]
+    set colnames {}
+    for {set i 1} {$i <= $fields_} {incr i} {
+        exec ftkeypar "${infile}+1" "TTYPE${i}"
+        set val_ [ exec pget ftkeypar value ]
+        set colname_ [ string toupper [ string trim $val_ {' } ] ]
+        lappend colnames $colname_
+    }
+    if { [ lsearch $colnames {COUNTS} ] != -1 } {
+        # Input spectrum is already in "COUNTS" format
+        return 0
+    } elseif { [ lsearch $colnames {RATE} ] != -1 } {
+        # Get exposure time
+        exec ftkeypar "${infile}+1" EXPOSURE
+        set exposure_ [ exec pget ftkeypar value ]
+        # Convert column "RATE"
+        set tmpfile "${outfile}_tmp[pid]"
+        exec cp -fv ${infile} ${tmpfile}
+        # COUNTS = RATE * EXPOSURE
+        set fcmd_ "ftcalc \"${tmpfile}\" \"${outfile}\" column=RATE expression=\"RATE * $exposure_\" history=yes"
+        puts $fcmd_
+        exec sh -c $fcmd_
+        exec rm -fv ${tmpfile}
+        # Update column name from "RATE" to "COUNTS", and column units
+        set field_idx_ [ expr { [ lsearch $colnames {RATE} ] + 1 } ]
+        set fcmd_ "fthedit \"${outfile}+1\" keyword=\"TTYPE${field_idx_}\" operation=add value=\"COUNTS\""
+        puts $fcmd_
+        exec sh -c $fcmd_
+        set fcmd_ "fthedit \"${outfile}+1\" keyword=\"TUNIT${field_idx_}\" operation=add value=\"count\""
+        puts $fcmd_
+        exec sh -c $fcmd_
+        # Convert column "STAT_ERR"
+        if { [ lsearch $colnames {STAT_ERR} ] != -1 } {
+            exec mv -fv ${outfile} ${tmpfile}
+            # STAT_ERR = STAT_ERR * EXPOSURE
+            set fcmd_ "ftcalc \"${tmpfile}\" \"${outfile}\" column=STAT_ERR expression=\"STAT_ERR * $exposure_\" history=yes"
+            puts $fcmd_
+            exec sh -c $fcmd_
+            exec rm -fv ${tmpfile}
+        }
+        return 1
+    } else {
+        # Invalid input spectrum
+        return -1
+    }
+}
+
+# Scale the spectrum "COUNTS" and "STAT_ERR" (if exists) columns by a factor.
+# NOTE:
+# This procedure is mainly used to scale the FWC background, within which
+# many channels have "counts" less than one.  When using "mathpha" to combine
+# this FWC background and other faked spectra in "COUNT" mode, the channels
+# with "counts" less than one (especially the higher energy regime) may be
+# truncated/rounded to integer, therefore the higher energy regime of the
+# combined spectrum may be greatly underestimated, and will cause the
+# following fitting very bad in the higher energy regime.
+# See also the description about the "SCALE_FACTOR" above.
+proc scale_spectrum {outfile infile factor} {
+    if {[ file exists ${outfile} ]} {
+        exec mv -fv ${outfile} ${outfile}_bak
+    }
+    exec cp -fv ${infile} ${outfile}
+    # Get number of columns/fields
+    exec ftkeypar "${outfile}+1" TFIELDS
+    set fields_ [ exec pget ftkeypar value ]
+    for {set i 1} {$i <= $fields_} {incr i} {
+        exec ftkeypar "${outfile}+1" "TTYPE${i}"
+        set val_ [ exec pget ftkeypar value ]
+        set colname_ [ string toupper [ string trim $val_ {' } ] ]
+        # Scale column "COUNTS"
+        if {$colname_ == "COUNTS"} {
+            set tmpfile "${outfile}_tmp[pid]"
+            exec mv -fv ${outfile} ${tmpfile}
+            # COUNTS = COUNTS * factor
+            set fcmd_ "ftcalc \"${tmpfile}\" \"${outfile}\" column=COUNTS expression=\"COUNTS * $factor\" history=yes"
+            puts $fcmd_
+            exec sh -c $fcmd_
+            exec rm -fv ${tmpfile}
+        }
+        # Scale column "STAT_ERR"
+        if {$colname_ == "STAT_ERR"} {
+            set tmpfile "${outfile}_tmp[pid]"
+            exec mv -fv ${outfile} ${tmpfile}
+            # STAT_ERR = STAT_ERR * factor
+            set fcmd_ "ftcalc \"${tmpfile}\" \"${outfile}\" column=STAT_ERR expression=\"STAT_ERR * $factor\" history=yes"
+            puts $fcmd_
+            exec sh -c $fcmd_
+            exec rm -fv ${tmpfile}
+        }
+    }
+}
+
+# Combine faked spectra to original FWC background spectrum with "mathpha"
+proc combine_spectra {outfile fwc_bkg cxb nxb exposure backscale} {
+    set combine_expr_ "$fwc_bkg + $cxb + $nxb"
+    set comment1_ "Corrected background spectrum; FWC based"
+    set combine_cmd_ "mathpha expr='${combine_expr_}' outfil=${outfile} exposure=${exposure} backscal=${backscale} units=C areascal=% properr=yes ncomments=1 comment1='${comment1_}'"
+    if {[ file exists ${outfile} ]} {
+        exec mv -fv ${outfile} ${outfile}_bak
+    }
+    puts $combine_cmd_
+    exec sh -c $combine_cmd_
+}
+
+# Copy header keywords from original background spectrum to the
+# newly combined background spectrum
+proc update_header {newfile origfile} {
+    global KEYWORDS
+    for {set i 0} {$i < [ llength $KEYWORDS ]} {incr i} {
+        set key_ [ lindex $KEYWORDS $i ]
+        exec ftkeypar "${origfile}+1" $key_
+        set val_ [ exec pget ftkeypar value ]
+        #set fcmd_ "fparkey value=\"$val_\" fitsfile=\"${newfile}+1\" keyword=\"$key_\" add=yes"
+        set fcmd_ "fthedit \"${newfile}+1\" keyword=\"$key_\" operation=add value=\"$val_\""
+        puts $fcmd_
+        exec sh -c $fcmd_
+    }
+}
+
+
+# Save current fit results
+save_xspec_fit
+
+# Get the current model components, which is the same for all data groups
+get_modelcomp
+
+for {set dg 1} {$dg <= $datagrp_total} {incr dg} {
+    get_filenames          $dg
+    get_exposure_backscale $dg
+    get_modelpars          $dg
+    get_backmodel          $dg
+}
+
+puts "modelcomp: $modelcomp"
+puts "rootname: $rootname"
+puts "rmf: $rmf"
+puts "arf: $arf"
+puts "diag_rsp: $diag_rsp"
+puts "src_spec: $src_spec"
+puts "bkg_spec: $bkg_spec"
+puts "src_exposure: $src_exposure"
+puts "bkg_exposure: $bkg_exposure"
+puts "src_backscale: $src_backscale"
+puts "bkg_backscale: $bkg_backscale"
+puts "modelpars: $modelpars"
+puts "back_modelname: $back_modelname"
+puts "back_modelcomp: $back_modelcomp"
+puts "back_modelpars: $back_modelpars"
+puts "SCALE_FACTOR: $SCALE_FACTOR"
+
+# DEBUG
+set DEBUG 0
+if {$DEBUG != 1} {
+
+# Clear all loaded model data and models
+data none
+model clear
+
+for {set idg 0} {$idg < $datagrp_total} {incr idg} {
+    set rootname_       [ lindex $rootname       $idg ]
+    set bkg_spec_       [ lindex $bkg_spec       $idg ]
+    set rmf_            [ lindex $rmf            $idg ]
+    set arf_            [ lindex $arf            $idg ]
+    set diag_rsp_       [ lindex $diag_rsp       $idg ]
+    set bkg_exposure_   [ lindex $bkg_exposure   $idg ]
+    set bkg_backscale_  [ lindex $bkg_backscale  $idg ]
+    set modelpars_      [ lindex $modelpars      $idg ]
+    set back_modelcomp_ [ lindex $back_modelcomp $idg ]
+    set back_modelpars_ [ lindex $back_modelpars $idg ]
+    set fake_exposure   [ expr {$bkg_exposure_  * $SCALE_FACTOR} ]
+    set fake_backscale  [ expr {$bkg_backscale_ * $SCALE_FACTOR} ]
+    set fake_cxb        "fake_cxb_${rootname_}.pi"
+    set fake_nxb        "fake_nxb_${rootname_}.pi"
+    set bkgspec_cnts    [ make_filename $bkg_spec_ "cnts" ]
+    set bkgcorr_outfile "bkgcorr_${rootname_}.pi"
+    # Fake CXB and NXB
+    load_model    $modelcomp $modelpars_
+    set_icm_norm
+    fake_spectrum $fake_cxb $rmf_ $arf_ $fake_exposure
+    load_model    $back_modelcomp_ $back_modelpars_
+    fake_spectrum $fake_nxb $diag_rsp_ "" $fake_exposure
+    # Convert FWC background spectrum from "RATE" to "COUNTS" for PN
+    set ret [ spectrum_rate2counts $bkgspec_cnts $bkg_spec_ ]
+    if { $ret == 0 } {
+        set bkgspec_cnts $bkg_spec_
+    } elseif { $ret == 1 } {
+        puts "Converted RATE spectrum to COUNTS: $bkgspec_cnts"
+    } else {
+        return -code error "ERROR: invalid spectrum for 'spectrum_rate2counts'"
+    }
+    set bkgspec_scaled [ make_filename $bkgspec_cnts "sf${SCALE_FACTOR}" ]
+    # Scale original FWC background before combination
+    scale_spectrum  $bkgspec_scaled $bkgspec_cnts $SCALE_FACTOR
+    # Combine faked spectra with original FWC background
+    combine_spectra $bkgcorr_outfile $bkgspec_scaled $fake_cxb $fake_nxb $bkg_exposure_ $fake_backscale
+    update_header   $bkgcorr_outfile $bkg_spec_
+}
+
+}
+
+# vim: set ts=8 sw=4 tw=0 fenc=utf-8 ft=tcl: #