Add 'renorm_spectrum.py' to replace 'chandra_bkg_rescale.sh'

Also remove 'chandra_pb_flux.sh'

3 files changed, 74 insertions, 149 deletions

diff --git a/scripts/chandra_bkg_rescale.sh b/scripts/chandra_bkg_rescale.sh
deleted file mode 100755
index 8d4abc7..0000000
--- a/scripts/chandra_bkg_rescale.sh
+++ /dev/null
@@ -1,80 +0,0 @@
-#!/bin/sh
-#
-# background rescale, by adjusting `BACKSCAL'
-# according to the photon flux values in `9.5-12.0 keV'
-#
-# Weitian LI <liweitianux@gmail.com>
-# 2012/08/14
-#
-# Changelogs:
-# 2015/06/03, Aaron LI
-#   * Copy needed pfiles to tmp directory,
-#     set environment variable $PFILES to use these first.
-#     and remove them after usage.
-#
-
-## background rescale (BACKSCAL) {{{
-# rescale 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_rescale() {
-    # $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 rescale }}}
-
-if [ $# -ne 2 ] || [ "x$1" = "x-h" ]; then
-    printf "usage:\n"
-    printf "    `basename $0` <src_spec> <bkg_spec>\n"
-    printf "\nNOTE:\n"
-    printf "<bkg_spec> is the spectrum to be adjusted\n"
-    exit 1
-fi
-
-## prepare parameter files (pfiles) {{{
-CIAO_TOOLS="dmstat dmkeypar dmhedit"
-
-PFILES_TMPDIR="/tmp/pfiles-$$"
-[ -d "${PFILES_TMPDIR}" ] && rm -rf ${PFILES_TMPDIR} || mkdir ${PFILES_TMPDIR}
-
-# Copy necessary pfiles for localized usage
-for tool in ${CIAO_TOOLS}; do
-    pfile=`paccess ${tool}`
-    [ -n "${pfile}" ] && punlearn ${tool} && cp -Lvf ${pfile} ${PFILES_TMPDIR}/
-done
-
-# Modify environment variable 'PFILES' to use local pfiles first
-export PFILES="${PFILES_TMPDIR}:${PFILES}"
-## pfiles }}}
-
-# perform `bkg_rescale'
-bkg_rescale "$1" "$2"
-
-# clean pfiles
-rm -rf ${PFILES_TMPDIR}
-
-exit 0
-
diff --git a/scripts/chandra_pb_flux.sh b/scripts/chandra_pb_flux.sh
deleted file mode 100755
index 59d4418..0000000
--- a/scripts/chandra_pb_flux.sh
+++ /dev/null
@@ -1,69 +0,0 @@
-#!/bin/sh
-#
-# chandra particle background
-# 9.5-12.0 keV (channel: 651-822)
-# PI = [ energy(eV) / 14.6 eV + 1 ]
-#
-# Weitian <liweitianux@gmail.com>
-# 2012/07/30
-#
-# ChangeLog:
-# v2.0, 2015/06/03, Aaron LI
-#   * Copy needed pfiles to tmp directory,
-#     set environment variable $PFILES to use these first.
-#     and remove them after usage.
-# v1.1: August 2, 2012
-#   * fix bugs with scientific notation in `bc'
-#
-
-if [ $# -eq 0 ] || [ "x$1" = "x-h" ]; then
-    echo "usage:"
-    echo "    `basename $0` <spec> ..."
-    exit 1
-fi
-
-## energy range: 9.5 -- 12.0 keV
-EN_LOW="9.5"
-EN_HI="12.0"
-CH_LOW=651
-CH_HI=822
-
-echo "Energy: $EN_LOW -- $EN_HI (keV)"
-echo "Channel: $CH_LOW -- $CH_HI"
-
-## prepare parameter files (pfiles) {{{
-CIAO_TOOLS="dmstat dmkeypar"
-
-PFILES_TMPDIR="/tmp/pfiles-$$"
-[ -d "${PFILES_TMPDIR}" ] && rm -rf ${PFILES_TMPDIR} || mkdir ${PFILES_TMPDIR}
-
-# Copy necessary pfiles for localized usage
-for tool in ${CIAO_TOOLS}; do
-    pfile=`paccess ${tool}`
-    [ -n "${pfile}" ] && punlearn ${tool} && cp -Lvf ${pfile} ${PFILES_TMPDIR}/
-done
-
-# Modify environment variable 'PFILES' to use local pfiles first
-export PFILES="${PFILES_TMPDIR}:${PFILES}"
-## pfiles }}}
-
-while ! [ -z $1 ]; do
-    f=$1
-    shift
-    echo "FILE: $f"
-    punlearn dmstat
-    COUNTS=`dmstat "$f[channel=${CH_LOW}:${CH_HI}][cols counts]" | grep 'sum:' | awk '{ print $2 }'`
-    punlearn dmkeypar
-    EXPTIME=`dmkeypar $f EXPOSURE echo=yes`
-    BACK=`dmkeypar $f 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 "    counts / exptime / backscal: ${COUNTS} / ${EXPTIME} / ${BACK}"
-    echo "    ${PB_FLUX}"
-done
-
-# clean pfiles
-rm -rf ${PFILES_TMPDIR}
-
diff --git a/scripts/renorm_spectrum.py b/scripts/renorm_spectrum.py
new file mode 100755
index 0000000..a09c90d
--- /dev/null
+++ b/scripts/renorm_spectrum.py
@@ -0,0 +1,74 @@
+#!/usr/bin/env python3
+#
+# Copyright (c) 2017 Weitian LI <liweitianux@live.com>
+# MIT license
+
+"""
+Renormalize the (background) spectrum by equaling its particle background
+flux (e.g., 9.5-12.0 keV) with respect to its corresponding source spectrum.
+
+The ``BACKSCAL`` keyword of the background spectrum is modified to realize
+the above renormalization.
+"""
+
+import sys
+import argparse
+import subprocess
+
+from spectrum import Spectrum
+
+
+def renorm_spectrum(specfile, specfile_ref, elow=9500, ehigh=12000):
+    """
+    Modify the ``BACKSCAL`` of ``specfile`` in order to equal its
+    particle background flux to that of ``specfile_ref``.
+
+    Parameters
+    ----------
+    specfile : str
+        (background) spectrum to be renormalized/modified.
+    specfile_ref : str
+        (source/reference) spectrum
+    elow, ehigh : float, optional
+        Lower and upper energy limit of the particle background.
+    """
+    spec = Spectrum(specfile)
+    spec_ref = Spectrum(specfile_ref)
+    flux = spec.calc_pb_flux(elow=elow, ehigh=ehigh)
+    flux_ref = spec_ref.calc_pb_flux(elow=elow, ehigh=ehigh)
+    bs_old = spec.BACKSCAL
+    bs_new = bs_old * flux / flux_ref
+    subprocess.check_call(["punlearn", "dmhedit"])
+    subprocess.check_call([
+        "dmhedit", "infile=%s" % specfile,
+        "filelist=none", "operation=add",
+        "key=BACKSCAL", "value=%s" % bs_new,
+        "comment='Old BACKSCAL: %s'" % bs_old
+    ])
+    print("%s:BACKSCAL: %f -> %f" % (specfile, bs_old, bs_new),
+          file=sys.stderr)
+
+
+def main():
+    parser = argparse.ArgumentParser(
+        description="Renormalize background spectrum w.r.t. source spectrum")
+    parser.add_argument("-L", "--energy-low", dest="elow",
+                        type=int, default=9500,
+                        help="Lower energy limit of the particle " +
+                        "background [eV] (default: 9500 eV)")
+    parser.add_argument("-H", "--energy-high", dest="ehigh",
+                        type=int, default=12000,
+                        help="Upper energy limit of the particle " +
+                        "background [eV] (default: 12000 eV)")
+    parser.add_argument("-r", "--spec-ref", dest="spec_ref", required=True,
+                        help="Reference (source) spectrum")
+    parser.add_argument("spec", dest="spec",
+                        help="(background) spectrum to be renormalized")
+    args = parser.parse_args()
+
+    renorm_spectrum(specfile=args.spec, specfile_ref=args.spec_ref,
+                    elow=args.elow, ehigh=args.ehigh)
+
+
+if __name__ == "__main__":
+    main()