Move various Python tools from 'scripts/' to 'bin/'

11 files changed, 1118 insertions, 0 deletions

diff --git a/bin/_context.py b/bin/_context.py
new file mode 100644
index 0000000..b23b2c8
--- /dev/null
+++ b/bin/_context.py
@@ -0,0 +1,16 @@
+# Copyright (c) 2017 Weitian LI <liweitianux@live.com>
+# MIT license
+
+"""
+Portal to 'acispy' module/package
+"""
+
+import os
+import sys
+
+sys.path.insert(
+    0,
+    os.path.abspath(os.path.join(os.path.dirname(__file__), ".."))
+)
+
+import acispy
diff --git a/bin/analyze_path.py b/bin/analyze_path.py
new file mode 100755
index 0000000..d2acd42
--- /dev/null
+++ b/bin/analyze_path.py
@@ -0,0 +1,52 @@
+#!/usr/bin/env python3
+#
+# Copyright (c) 2017 Weitian LI <liweitianux@live.com>
+# MIT license
+#
+# Weitian LI
+# 2017-02-06
+
+"""
+Extract the object name and observation ID from the directory path.
+
+The root directory of the object data has the format:
+    <name>_oi<obsid>
+"""
+
+import os
+import argparse
+
+from _context import acispy
+from acispy.analyze_path import get_name, get_obsid
+
+
+def main():
+    parser = argparse.ArgumentParser(
+        description="Extract object name and ObsID from data directory")
+    parser.add_argument("-b", "--brief", dest="brief",
+                        action="store_true", help="Be brief")
+    parser.add_argument("-n", "--name", dest="name",
+                        action="store_true", help="Only get object name")
+    parser.add_argument("-i", "--obsid", dest="obsid",
+                        action="store_true", help="Only get observation ID")
+    parser.add_argument("path", nargs="?", default=os.getcwd(),
+                        help="Path to the data directory " +
+                        "(default: current working directory)")
+    args = parser.parse_args()
+    path = os.path.abspath(args.path)
+
+    b_get_name = False if args.obsid else True
+    b_get_obsid = False if args.name else True
+
+    if b_get_name:
+        if not args.brief:
+            print("Name:", end=" ")
+        print(get_name(path))
+    if b_get_obsid:
+        if not args.brief:
+            print("ObsID:", end=" ")
+        print(get_obsid(path))
+
+
+if __name__ == "__main__":
+    main()
diff --git a/bin/calc_centroid.py b/bin/calc_centroid.py
new file mode 100755
index 0000000..7df71b1
--- /dev/null
+++ b/bin/calc_centroid.py
@@ -0,0 +1,184 @@
+#!/usr/bin/env python3
+#
+# Copyright (c) 2017 Weitian LI <liweitianux@live.com>
+# MIT license
+
+"""
+Calculate the coordinate of the emission centroid within the image.
+
+The image are smoothed first, and then an iterative procedure with
+two phases is applied to determine the emission centroid.
+"""
+
+import os
+import sys
+import argparse
+import subprocess
+import tempfile
+
+from _context import acispy
+from acispy.manifest import get_manifest
+from acispy.pfiles import setup_pfiles
+from acispy.ds9 import ds9_view
+from acispy.region import Regions
+
+
+def smooth_image(infile, outfile=None,
+                 kernelspec="lib:gaus(2,5,1,10,10)", method="fft",
+                 clobber=False):
+    """
+    Smooth the image by a Gaussian kernel using the ``aconvolve`` tool.
+
+    Parameters
+    ----------
+    infile : str
+        Path to the input image file
+    outfile : str, optional
+        Filename/path of the output smoothed image
+        (default: build in format ``<infile_basename>_aconv.fits``)
+    kernelspec : str, optional
+        Kernel specification for ``aconvolve``
+    method : str, optional
+        Smooth method for ``aconvolve``
+
+    Returns
+    -------
+    outfile : str
+        Filename/path of the smoothed image
+    """
+    clobber = "yes" if clobber else "no"
+    if outfile is None:
+        outfile = os.path.splitext(infile)[0] + "_aconv.fits"
+    subprocess.check_call(["punlearn", "aconvolve"])
+    subprocess.check_call([
+        "aconvolve", "infile=%s" % infile, "outfile=%s" % outfile,
+        "kernelspec=%s" % kernelspec, "method=%s" % method,
+        "clobber=%s" % clobber
+    ])
+    return outfile
+
+
+def get_peak(image):
+    """
+    Get the peak coordinate on the image.
+
+    Returns
+    -------
+    peak : 2-float tuple
+        (Physical) coordinate of the peak.
+    """
+    subprocess.check_call(["punlearn", "dmstat"])
+    subprocess.check_call([
+        "dmstat", "infile=%s" % image,
+        "centroid=no", "media=no", "sigma=no", "clip=no", "verbose=0"
+    ])
+    peak = subprocess.check_output([
+        "pget", "dmstat", "out_max_loc"
+    ]).decode("utf-8").strip()
+    peak = peak.split(",")
+    return (float(peak[0]), float(peak[1]))
+
+
+def get_centroid(image, center, radius=50):
+    """
+    Calculate the centroid on image within the specified circle.
+
+    Parameters
+    ----------
+    image : str
+        Path to the image file.
+    center : 2-float tuple
+        Central (physical) coordinate of the circle.
+    radius : float
+        Radius (pixel) of the circle.
+
+    Returns
+    -------
+    centroid : 2-float tuple
+        (Physical) coordinate of the centroid.
+    """
+    x, y = center
+    with tempfile.NamedTemporaryFile(mode="w+") as fp:
+        fp.file.write("circle(%f,%f,%f)\n" % (x, y, radius))
+        fp.file.flush()
+        subprocess.check_call(["punlearn", "dmstat"])
+        subprocess.check_call([
+            "dmstat", "infile=%s[sky=region(%s)]" % (image, fp.name),
+            "centroid=yes", "media=no", "sigma=no", "clip=no", "verbose=0"
+        ])
+    centroid = subprocess.check_output([
+        "pget", "dmstat", "out_cntrd_phys"
+    ]).decode("utf-8").strip()
+    centroid = centroid.split(",")
+    return (float(centroid[0]), float(centroid[1]))
+
+
+def main():
+    parser = argparse.ArgumentParser(
+        description="Calculate the emission centroid within the image")
+    parser.add_argument("-i", "--infile", dest="infile", required=True,
+                        help="input image file (e.g., 0.7-2.0 keV)")
+    parser.add_argument("-o", "--outfile", dest="outfile",
+                        default="centroid.reg",
+                        help="output centroid region file " +
+                        "(default: centroid.reg")
+    parser.add_argument("-R", "--radius1", dest="radius1",
+                        type=float, default=100,
+                        help="circle radius [pixel] for first phase " +
+                        "centroid calculation (default: 100 pixel)")
+    parser.add_argument("-r", "--radius2", dest="radius2",
+                        type=float, default=50,
+                        help="circle radius [pixel] for second phase " +
+                        "calculation to tune centroid (default: 50 pixel)")
+    parser.add_argument("-n", "--niter", dest="niter",
+                        type=int, default=10,
+                        help="iterations for each phase (default: 10)")
+    parser.add_argument("-s", "--start", dest="start",
+                        help="a region file containing a circle/point " +
+                        "that specifies the starting point " +
+                        "(default: using the peak of the image)")
+    parser.add_argument("-V", "--view", dest="view", action="store_true",
+                        help="open DS9 to view output centroid")
+    parser.add_argument("-C", "--clobber", dest="clobber", action="store_true",
+                        help="overwrite existing files")
+    args = parser.parse_args()
+
+    setup_pfiles(["aconvolve", "dmstat"])
+
+    print("Smooth input image using 'aconvolve' ...", file=sys.stderr)
+    img_smoothed = smooth_image(args.infile, clobber=args.clobber)
+
+    if args.start:
+        print("Get starting point from region file: %s" % args.start,
+              file=sys.stderr)
+        region = Regions(args.start).regions[0]
+        center = (region.xc, region.yc)
+    else:
+        print("Use peak as the starting point ...", file=sys.stderr)
+        center = get_peak(img_smoothed)
+    print("Starting point: (%f, %f)" % center, file=sys.stderr)
+
+    centroid = center
+    for phase, radius in enumerate([args.radius1, args.radius2]):
+        print("Calculate centroid phase %d (circle radius: %.1f)" %
+              (phase+1, radius), file=sys.stderr)
+        for i in range(args.niter):
+            print("%d..." % (i+1), end="", flush=True, file=sys.stderr)
+            centroid = get_centroid(img_smoothed, center=centroid,
+                                    radius=radius)
+        print("Done!", file=sys.stderr)
+
+    with open(args.outfile, "w") as f:
+        f.write("point(%f,%f)\n" % centroid)
+    print("Saved centroid to file:", args.outfile, file=sys.stderr)
+    if args.view:
+        ds9_view(img_smoothed, regfile=args.outfile)
+
+    # Add calculated centroid region to manifest
+    manifest = get_manifest()
+    key = "reg_centroid"
+    manifest.setpath(key, args.outfile)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/bin/clean_evt2.py b/bin/clean_evt2.py
new file mode 100755
index 0000000..5b53922
--- /dev/null
+++ b/bin/clean_evt2.py
@@ -0,0 +1,207 @@
+#!/usr/bin/env python3
+#
+# Copyright (c) 2017 Weitian LI <liweitianux@live.com>
+# MIT license
+
+"""
+Further process the reprocessed evt2 file to make a cleaned one
+for later scientific analysis.
+
+The following steps are carried out:
+1. Filter out the chips of interest;
+2. Detect (point) sources, visually check and manually update;
+3. Remove the updated sources;
+4. Extract light curve of the background regions (away from the object);
+5. Check and clip the light curve to create a GTI;
+6. Remove flares by filtering the GTI to create the finally cleaned evt2.
+"""
+
+import os
+import argparse
+import subprocess
+import shutil
+import logging
+
+from _context import acispy
+from acispy.manifest import get_manifest
+from acispy.pfiles import setup_pfiles
+from acispy.acis import ACIS
+from acispy.ds9 import ds9_view
+
+
+logging.basicConfig(level=logging.INFO)
+logger = logging.getLogger(__name__)
+
+
+def filter_chips(infile, outfile, chips, clobber=False):
+    """
+    Filter out the chips of interest, e.g., ``ccd_id=7`` for ACIS-S
+    observation, and ``ccd_id=0:3`` for ACIS-I observation.
+    """
+    chips = chips.replace("-", ":")
+    clobber = "yes" if clobber else "no"
+    logger.info("Filter out chips of interest: %s" % chips)
+    logger.info("Outfile: %s" % outfile)
+    subprocess.check_call(["punlearn", "dmcopy"])
+    subprocess.check_call([
+        "dmcopy", "infile=%s[ccd_id=%s]" % (infile, chips),
+        "outfile=%s" % outfile, "clobber=%s" % clobber
+    ])
+    logger.info("Done!\n")
+
+
+def detect_sources(infile, outfile, clobber=False):
+    """
+    Detect (point) sources using ``celldetect``; examine the detected
+    sources in DS9, and update the source regions and save.
+    """
+    src_fits = os.path.splitext(outfile)[0] + ".fits"
+    clobber = "yes" if clobber else "no"
+    logger.info("Detect sources using 'celldetect' ...")
+    logger.info("Outfile: %s" % outfile)
+    subprocess.check_call(["punlearn", "celldetect"])
+    subprocess.check_call([
+        "celldetect", "infile=%s" % infile,
+        "outfile=%s" % src_fits, "regfile=%s" % outfile,
+        "clobber=%s" % clobber
+    ])
+    os.remove(src_fits)
+    shutil.copy(outfile, outfile+".orig")
+    logger.warning("Check and update detected source regions; " +
+                   "and save/overwrite file: %s" % outfile)
+    ds9_view(infile, regfile=outfile)
+    logger.info("Done!\n")
+
+
+def remove_sources(infile, outfile, srcfile, clobber=False):
+    """
+    Remove detected sources
+    """
+    clobber = "yes" if clobber else "no"
+    logger.info("Remove detected sources ...")
+    logger.info("Outfile: %s" % outfile)
+    subprocess.check_call(["punlearn", "dmcopy"])
+    subprocess.check_call([
+        "dmcopy", "infile=%s[exclude sky=region(%s)]" % (infile, srcfile),
+        "outfile=%s" % outfile, "clobber=%s" % clobber
+    ])
+    logger.info("Done!\n")
+
+
+def extract_lightcurve(infile, outfile, bintime=200, clobber=False):
+    """
+    Extract the light curve from regions away from the object.
+    """
+    clobber = "yes" if clobber else "no"
+    logger.info("Extract light curve for GTI generation ...")
+    logger.info("Outfile: %s" % outfile)
+    regfile = os.path.splitext(outfile)[0] + ".reg"
+    # Credit: https://stackoverflow.com/a/12654798/4856091
+    open(regfile, "a").close()
+    logger.warning("Select a large region containing most of the object, " +
+                   "but also leaving enough area outside as background; " +
+                   "and save as file: %s" % regfile)
+    ds9_view(infile)
+    fsky = "exclude sky=region(%s)" % regfile
+    fbintime = "bin time=::%s" % bintime
+    subprocess.check_call(["punlearn", "dmextract"])
+    subprocess.check_call([
+        "dmextract", "infile=%s[%s][%s]" % (infile, fsky, fbintime),
+        "outfile=%s" % outfile, "opt=ltc1", "clobber=%s" % clobber
+    ])
+    logger.info("Done!\n")
+
+
+def make_gti(infile, outfile, scale=1.2, clobber=False):
+    """
+    Examine the light curve for flares, and clip it to make the GTI.
+    """
+    logger.info("Examine the light curve and create GTI ...")
+    logger.info("Outfile: %s" % outfile)
+
+    chipsfile = os.path.splitext(outfile)[0] + ".chips"
+    if (not clobber) and (os.path.exists(outfile) or
+                          os.path.exists(chipsfile)):
+        raise OSError("'%s' or '%s' already exists" % (outfile, chipsfile))
+
+    outimg = os.path.splitext(outfile)[0] + "_lc.jpg"
+    lines = [
+        "from lightcurves import lc_clean",
+        "lc_clean('%s')" % infile,
+        "lc_clean('%s', scale=%s, outfile='%s')" % (infile, scale, outfile),
+        "print_window('%s', ['format', 'jpg', 'clobber', 'True'])" % outimg
+    ]
+    open(chipsfile, "w").write("\n".join(lines) + "\n")
+    subprocess.check_call(["chips", "-x", chipsfile])
+
+    if not os.path.exists(outfile):
+        # workaround the problem that ``chips`` sometimes just failed
+        logger.warning("*** Failed to create GTI: %s ***" % outfile)
+        logger.warning("You need to create the GTI manually.")
+        input("When finished GTI creation, press Enter to continue ...")
+    logger.info("Done!\n")
+
+
+def filter_gti(infile, outfile, gti, clobber=False):
+    """
+    Removing flares by filtering on GTI to create the finally cleaned
+    evt2 file.
+    """
+    clobber = "yes" if clobber else "no"
+    logger.info("Filter on GTI to create cleaned evt2 file ...")
+    logger.info("Outfile: %s" % outfile)
+    subprocess.check_call(["punlearn", "dmcopy"])
+    subprocess.check_call([
+        "dmcopy", "infile=%s[@%s]" % (infile, gti),
+        "outfile=%s" % outfile, "clobber=%s" % clobber
+    ])
+    logger.info("Done!\n")
+
+
+def main():
+    parser = argparse.ArgumentParser(
+        description="Make a cleaned evt2 for scientific analysis")
+    parser.add_argument("-i", "--infile", dest="infile",
+                        help="input evt2 produced by 'chandra_repro' " +
+                        "(default: request from manifest)")
+    parser.add_argument("-C", "--clobber", dest="clobber", action="store_true",
+                        help="overwrite existing file")
+    args = parser.parse_args()
+
+    setup_pfiles(["dmkeypar", "dmcopy", "celldetect", "dmextract"])
+
+    manifest = get_manifest()
+    if args.infile:
+        infile = args.infile
+    else:
+        infile = manifest.getpath("evt2", relative=True)
+    chips = ACIS.get_chips_str(infile, sep="-")
+    logger.info("infile: %s" % infile)
+    logger.info("chips: %s" % chips)
+
+    evt2_chips = "evt2_c{chips}_orig.fits".format(chips=chips)
+    evt2_rmsrc = "evt2_c{chips}_rmsrc.fits".format(chips=chips)
+    evt2_clean = "evt2_c{chips}_clean.fits".format(chips=chips)
+    srcfile = "sources_celld.reg"
+    lcfile = "ex_bkg.lc"
+    gtifile = os.path.splitext(lcfile)[0] + ".gti"
+
+    filter_chips(infile, evt2_chips, chips, clobber=args.clobber)
+    detect_sources(evt2_chips, srcfile, clobber=args.clobber)
+    remove_sources(evt2_chips, evt2_rmsrc, srcfile, clobber=args.clobber)
+    extract_lightcurve(evt2_rmsrc, lcfile, clobber=args.clobber)
+    make_gti(lcfile, gtifile, clobber=args.clobber)
+    filter_gti(evt2_rmsrc, evt2_clean, gtifile, clobber=args.clobber)
+
+    # Add cleaned evt2 to manifest
+    key = "evt2_clean"
+    manifest.setpath(key, evt2_clean)
+    logger.info("Added '%s' to manifest: %s" % (key, manifest.get(key)))
+
+    # Remove useless intermediate files
+    os.remove(evt2_chips)
+    os.remove(evt2_rmsrc)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/bin/collect_yaml.py b/bin/collect_yaml.py
new file mode 100755
index 0000000..3147a0c
--- /dev/null
+++ b/bin/collect_yaml.py
@@ -0,0 +1,74 @@
+#!/usr/bin/env python3
+#
+# Copyright (c) 2017 Weitian LI <liweitianux@live.com>
+# MIT license
+#
+# Weitian LI
+# 2017-02-12
+
+"""
+Collect YAML manifest files, and convert collected results to CSV
+format for later use.
+"""
+
+import sys
+import argparse
+import csv
+
+from _context import acispy
+from acispy.manifest import Manifest
+
+
+def main():
+    parser = argparse.ArgumentParser(description="Collect YAML manifest files")
+    parser.add_argument("-k", "--keys", dest="keys", required=True,
+                        help="YAML keys to be collected (in order); " +
+                        "can be comma-separated string, or a file " +
+                        "containing the keys one-per-line")
+    parser.add_argument("-b", "--brief", dest="brief",
+                        action="store_true",
+                        help="be brief and do not print header")
+    parser.add_argument("-v", "--verbose", dest="verbose",
+                        action="store_true",
+                        help="show verbose information")
+    parser.add_argument("-o", "--outfile", dest="outfile", default=sys.stdout,
+                        help="output CSV file to save collected data")
+    parser.add_argument("-i", "--infile", dest="infile",
+                        nargs="+", required=True,
+                        help="list of input YAML manifest files")
+    args = parser.parse_args()
+
+    try:
+        keys = [k.strip() for k in open(args.keys).readlines()]
+    except FileNotFoundError:
+        keys = [k.strip() for k in args.keys.split(",")]
+
+    if args.verbose:
+        print("keys:", keys, file=sys.stderr)
+        print("infile:", args.infile, file=sys.stderr)
+        print("outfile:", args.outfile, file=sys.stderr)
+
+    results = []
+    for fp in args.infile:
+        manifest = Manifest(fp)
+        res = manifest.gets(keys, splitlist=True)
+        if args.verbose:
+            print("FILE:{0}: {1}".format(fp, list(res.values())),
+                  file=sys.stderr)
+        results.append(res)
+
+    try:
+        of = open(args.outfile, "w")
+    except TypeError:
+        of = args.outfile
+    writer = csv.writer(of)
+    if not args.brief:
+        writer.writerow(results[0].keys())
+    for res in results:
+        writer.writerow(res.values())
+    if of is not sys.stdout:
+        of.close()
+
+
+if __name__ == "__main__":
+    main()
diff --git a/bin/cosmo_calc.py b/bin/cosmo_calc.py
new file mode 100755
index 0000000..1302bdb
--- /dev/null
+++ b/bin/cosmo_calc.py
@@ -0,0 +1,133 @@
+#!/usr/bin/env python3
+#
+# Copyright (c) 2017 Weitian LI <liweitianux@live.com>
+# MIT license
+
+"""
+Cosmology calculator with support of Chandra ACIS-specific quantities.
+"""
+
+import sys
+import argparse
+from collections import OrderedDict
+
+from _context import acispy
+from acispy.cosmo import Calculator
+
+
+# Supported quantities
+QUANTITIES = OrderedDict([
+    ("luminosity_distance", {
+        "unit": "Mpc",
+        "label": "Luminosity distance",
+        "kwargs": ["z", "unit"],
+    }),
+    ("angular_diameter_distance", {
+        "unit": "Mpc",
+        "label": "Angular diameter distance",
+        "kwargs": ["z", "unit"],
+    }),
+    ("kpc_per_arcsec", {
+        "unit": None,
+        "label": "kpc/arcsec (DA)",
+        "kwargs": ["z"],
+    }),
+    ("kpc_per_pix", {
+        "unit": None,
+        "label": "kpc/pix (DA)",
+        "kwargs": ["z"],
+    }),
+    ("cm_per_pix", {
+        "unit": None,
+        "label": "cm/pix (DA)",
+        "kwargs": ["z"],
+    }),
+    ("norm_apec", {
+        "unit": "cm^-5",
+        "label": "norm (APEC)",
+        "kwargs": ["z"],
+    }),
+])
+
+
+def get_quantities(args):
+    # convert ``argparse.Namespace`` to a dictionary
+    args = vars(args)
+    q_all = list(QUANTITIES.keys())
+    q_active = [q for q in q_all if args[q]]
+    if len(q_active) == 0:
+        q_active = q_all
+    return q_active
+
+
+def calc_quantity(q, calculator, args):
+    args = vars(args)
+    kwargs = {arg: args[arg] for arg in QUANTITIES[q]["kwargs"]
+              if args[arg] is not None}
+    value = getattr(calculator, q)(**kwargs)
+    label = QUANTITIES[q]["label"]
+    unit = args["unit"] if args["unit"] is not None else QUANTITIES[q]["unit"]
+    if args["brief"]:
+        print(value)
+    else:
+        print("%s: %s  # [%s]" % (label, value, unit))
+
+
+def main():
+    parser = argparse.ArgumentParser(
+        description="Cosmology calculator with Chandra-specific quantities")
+    parser.add_argument("-H", "--hubble", dest="H0",
+                        type=float, default=71.0,
+                        help="Present-day Hubble parameter " +
+                        "(default: 71 km/s/Mpc)")
+    parser.add_argument("-M", "--omega-m", dest="Om0",
+                        type=float, default=0.27,
+                        help="Present-day matter density (default: 0.27")
+    parser.add_argument("-b", "--brief", dest="brief", action="store_true",
+                        help="be brief")
+    parser.add_argument("-U", "--unit", dest="unit",
+                        help="unit for output quantity if supported")
+    parser.add_argument("-L", "--luminosity-distance",
+                        dest="luminosity_distance",
+                        action="store_true",
+                        help="calculate the luminosity distance (DL)")
+    parser.add_argument("-A", "--angular-diameter-distance",
+                        dest="angular_diameter_distance",
+                        action="store_true",
+                        help="calculate the angular diameter distance (DA)")
+    parser.add_argument("--kpc-per-arcsec", dest="kpc_per_arcsec",
+                        action="store_true",
+                        help="calculate the transversal length [kpc] " +
+                        "w.r.t. 1 arcsec at DA(z)")
+    parser.add_argument("--kpc-per-pix", dest="kpc_per_pix",
+                        action="store_true",
+                        help="calculate the transversal length [kpc] " +
+                        "w.r.t. 1 ACIS pixel (0.492 arcsec) at DA(z)")
+    parser.add_argument("--cm-per-pix", dest="cm_per_pix",
+                        action="store_true",
+                        help="calculate the transversal length [cm] " +
+                        "w.r.t. 1 ACIS pixel (0.492 arcsec) at DA(z)")
+    parser.add_argument("--norm-apec", dest="norm_apec",
+                        action="store_true",
+                        help="calculate the normalization factor " +
+                        "of the XSPEC APEC model assuming EM=1")
+    parser.add_argument("z", type=float, help="redshift")
+    args = parser.parse_args()
+
+    cosmocalc = Calculator(H0=args.H0, Om0=args.Om0)
+
+    q_active = get_quantities(args)
+    if len(q_active) > 1:
+        if args.unit is not None:
+            args.unit = None
+            print("WARNING: ignored argument --unit", file=sys.stderr)
+        if args.brief:
+            args.brief = False
+            print("WARNING: ignored argument --brief", file=sys.stderr)
+
+    for q in q_active:
+        calc_quantity(q=q, calculator=cosmocalc, args=args)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/bin/event2image.py b/bin/event2image.py
new file mode 100755
index 0000000..44189e0
--- /dev/null
+++ b/bin/event2image.py
@@ -0,0 +1,105 @@
+#!/usr/bin/env python3
+#
+# Copyright (c) 2017 Weitian LI <liweitianux@live.com>
+# MIT license
+
+"""
+Make image by binning the event file, and update the manifest.
+
+TODO: use logging module instead of print()
+"""
+
+import sys
+import argparse
+import subprocess
+
+from _context import acispy
+from acispy.manifest import get_manifest
+from acispy.pfiles import setup_pfiles
+from acispy.acis import ACIS
+
+
+def make_image(infile, outfile, chips, erange, fov, clobber=False):
+    """
+    Make image by binning the event file.
+
+    Parameters
+    ----------
+    infile : str
+        Path to the input event file
+    outfile : str
+        Filename and path of the output image file
+    chips : str
+        Chips of interest, e.g., ``7`` or ``0-3``
+    erange : str
+        Energy range of interest, e.g., ``700-7000``
+    fov : str
+        Path to the FoV file
+    """
+    chips = chips.replace("-", ":")
+    erange = erange.replace("-", ":")
+    clobber = "yes" if clobber else "no"
+    fregion = "sky=region(%s[ccd_id=%s])" % (fov, chips)
+    fenergy = "energy=%s" % erange
+    fbin = "bin sky=::1"
+    subprocess.check_call(["punlearn", "dmcopy"])
+    subprocess.check_call([
+        "dmcopy", "infile=%s[%s][%s][%s]" % (infile, fregion, fenergy, fbin),
+        "outfile=%s" % outfile, "clobber=%s" % clobber
+    ])
+
+
+def main():
+    parser = argparse.ArgumentParser(
+        description="Make image by binning the event file")
+    parser.add_argument("-L", "--elow", dest="elow", type=int, default=700,
+                        help="lower energy limit [eV] of the output image " +
+                        "(default: 700 [eV])")
+    parser.add_argument("-H", "--ehigh", dest="ehigh", type=int,
+                        help="upper energy limit [eV] of the output image " +
+                        "(default: 7000 [eV])")
+    parser.add_argument("-i", "--infile", dest="infile",
+                        help="event file from which to create the image " +
+                        "(default: evt2_clean from manifest)")
+    parser.add_argument("-o", "--outfile", dest="outfile",
+                        help="output image filename (default: " +
+                        "build in format 'img_c<chip>_e<elow>-<ehigh>.fits')")
+    parser.add_argument("-v", "--verbose", dest="verbose", action="store_true",
+                        help="show verbose information")
+    parser.add_argument("-C", "--clobber", dest="clobber", action="store_true",
+                        help="overwrite existing file")
+    args = parser.parse_args()
+
+    setup_pfiles(["dmkeypar", "dmcopy"])
+
+    manifest = get_manifest()
+    fov = manifest.getpath("fov", relative=True)
+    if args.infile:
+        infile = args.infile
+    else:
+        infile = manifest.getpath("evt2_clean", relative=True)
+    chips = ACIS.get_chips_str(infile, sep="-")
+    erange = "{elow}-{ehigh}".format(elow=args.elow, ehigh=args.ehigh)
+    if args.elow >= args.ehigh:
+        raise ValueError("invalid energy range: %s" % erange)
+    if args.outfile:
+        outfile = args.outfile
+    else:
+        outfile = "img_c{chips}_e{erange}.fits".format(
+            chips=chips, erange=erange)
+    if args.verbose:
+        print("infile:", infile, file=sys.stderr)
+        print("outfile:", outfile, file=sys.stderr)
+        print("fov:", fov, file=sys.stderr)
+        print("chips:", chips, file=sys.stderr)
+        print("erange:", erange, file=sys.stderr)
+
+    make_image(infile, outfile, chips, erange, fov, args.clobber)
+
+    # Add created image to manifest
+    key = "img_e{erange}".format(erange=erange)
+    manifest.setpath(key, outfile)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/bin/make_blanksky.py b/bin/make_blanksky.py
new file mode 100755
index 0000000..8229d39
--- /dev/null
+++ b/bin/make_blanksky.py
@@ -0,0 +1,228 @@
+#!/usr/bin/env python3
+#
+# Copyright (c) 2017 Weitian LI <liweitianux@live.com>
+# MIT license
+
+"""
+Make the "blank-sky" background for spectral analysis.
+
+This tool first finds the corresponding blank-sky files for this observation,
+then apply further filtering and gain corrections, and finally reproject
+the coordinates to matching the observational data.
+
+
+Reference
+---------
+* Analyzing the ACIS Background with "Blank-sky" Files
+  http://cxc.harvard.edu/ciao/threads/acisbackground/
+"""
+
+import os
+import argparse
+import subprocess
+import shutil
+import stat
+import logging
+
+from _context import acispy
+from acispy.manifest import get_manifest
+from acispy.pfiles import setup_pfiles
+from acispy.acis import ACIS
+from acispy.header import write_keyword, copy_keyword
+
+
+logging.basicConfig(level=logging.INFO)
+logger = logging.getLogger(__name__)
+
+
+def find_blanksky(infile, copy=True, clobber=False):
+    """
+    Find the blanksky files for the input file using ``acis_bkgrnd_lookup``
+    tool, and copy the files to current directory if ``copy=True``.
+    """
+    subprocess.check_call(["punlearn", "acis_bkgrnd_lookup"])
+    subprocess.check_call(["acis_bkgrnd_lookup", "infile=%s" % infile])
+    bkgfiles = subprocess.check_output([
+        "pget", "acis_bkgrnd_lookup", "outfile"
+    ]).decode("utf-8").strip()
+    bkgfiles = bkgfiles.split(",")
+    logger.info("Found blanksky files: {0}".format(bkgfiles))
+    if copy:
+        logger.info("Copy blanksky files into CWD ...")
+        for f in bkgfiles:
+            if os.path.exists(os.path.basename(f)) and (not clobber):
+                raise OSError("File already exists: %s" % os.path.basename(f))
+            shutil.copy(f, ".")
+        bkgfiles = [os.path.basename(f) for f in bkgfiles]
+    return bkgfiles
+
+
+def merge_blanksky(infiles, outfile, clobber=False):
+    """
+    Merge multiple blanksky files (e.g., 4 blanksky files for ACIS-I)
+    into a single file.
+    """
+    if len(infiles) == 1:
+        logger.info("Only one blanksky file, no need to merge.")
+        if os.path.exists(outfile) and (not clobber):
+            raise OSError("File already exists: %s" % outfile)
+        shutil.move(infiles[0], outfile)
+    else:
+        logger.info("Merge multiple blanksky files ...")
+        clobber = "yes" if clobber else "no"
+        subprocess.check_call(["punlearn", "dmmerge"])
+        subprocess.check_call([
+            "dmmerge", "infile=%s" % ",".join(infiles),
+            "outfile=%s" % outfile, "clobber=%s" % clobber
+        ])
+        write_keyword(outfile, keyword="DETNAM", value="ACIS-0123")
+        for f in infiles:
+            os.remove(f)
+    # Add write permission to the background file
+    st = os.stat(outfile)
+    os.chmod(outfile, st.st_mode | stat.S_IWUSR)
+
+
+def clean_vfaint(infile, outfile, clobber=False):
+    """
+    Clean the background file by only keeping events with ``status=0``
+    for ``VFAINT`` mode observations.
+    """
+    subprocess.check_call(["punlearn", "dmkeypar"])
+    datamode = subprocess.check_output([
+        "dmkeypar", "infile=%s" % infile,
+        "keyword=DATAMODE", "echo=yes"
+    ]).decode("utf-8").strip()
+    logger.info("DATAMODE: %s" % datamode)
+    if datamode.upper() == "VFAINT":
+        logger.info("Clean background file using 'status=0' ...")
+        clobber = "yes" if clobber else "no"
+        subprocess.check_call(["punlearn", "dmcopy"])
+        subprocess.check_call([
+            "dmcopy", "infile=%s[status=0]" % infile,
+            "outfile=%s" % outfile, "clobber=%s" % clobber
+        ])
+        os.remove(infile)
+    else:
+        logger.info("No need to clean the background file.")
+        if os.path.exists(outfile) and (not clobber):
+            raise OSError("File already exists: %s" % outfile)
+        shutil.move(infile, outfile)
+
+
+def apply_gainfile(infile, outfile, evt2, clobber=False):
+    """
+    Check whether the ``GAINFILE`` of the background file matches that
+    of the observational evt2 file.
+
+    If they are different, then reprocess the background file with the
+    same ``GAINFILE`` of the evt2 file.
+    """
+    subprocess.check_call(["punlearn", "dmkeypar"])
+    gainfile_bkg = subprocess.check_output([
+        "dmkeypar", "infile=%s" % infile,
+        "keyword=GAINFILE", "echo=yes"
+    ]).decode("utf-8").strip()
+    gainfile_evt2 = subprocess.check_output([
+        "dmkeypar", "infile=%s" % evt2,
+        "keyword=GAINFILE", "echo=yes"
+    ]).decode("utf-8").strip()
+    if gainfile_bkg == gainfile_evt2:
+        logger.info("GAINFILE's are the same for background and evt2.")
+        if os.path.exists(outfile) and (not clobber):
+            raise OSError("File already exists: %s" % outfile)
+        shutil.move(infile, outfile)
+    else:
+        # Reprocess the background
+        gainfile = os.path.join(os.environ["CALDB"],
+                                "data/chandra/acis/det_gain",
+                                gainfile_evt2)
+        logger.info("Reprocess background using gainfile: %s ..." % gainfile)
+        clobber = "yes" if clobber else "no"
+        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"
+        ]
+        subprocess.check_call(["punlearn", "acis_process_events"])
+        subprocess.check_call([
+            "acis_process_events", "infile=%s" % infile,
+            "outfile=%s" % outfile, "acaofffile=NONE", "stop=none",
+            "doevtgrade=no", "apply_cti=yes", "apply_tgain=no",
+            "calculate_pi=yes", "pix_adj=NONE", "gainfile=%s" % gainfile,
+            "eventdef={%s}" % ",".join(eventdef),
+            "clobber=%s" % clobber
+        ])
+        os.remove(infile)
+
+
+def reproject_blanksky(infile, outfile, evt2, asol, clobber=False):
+    """
+    Reproject the background to match the coordinates of the observational
+    evt2 data.
+    """
+    clobber = "yes" if clobber else "no"
+    logger.info("Reprojecting the background to match the evt2 file ...")
+    subprocess.check_call(["punlearn", "reproject_events"])
+    subprocess.check_call([
+        "reproject_events", "infile=%s" % infile, "outfile=%s" % outfile,
+        "aspect=%s" % asol, "match=%s" % evt2, "random=0",
+        "clobber=%s" % clobber
+    ])
+    os.remove(infile)
+
+
+def main():
+    parser = argparse.ArgumentParser(
+        description="Make blanksky background")
+    parser.add_argument("-i", "--infile", dest="infile",
+                        help="input evt2 file " +
+                        "(default: 'evt2_clean' from manifest)")
+    parser.add_argument("-o", "--outfile", dest="outfile",
+                        help="output blanksky background file " +
+                        "(default: 'blanksky_c{chips}.fits')")
+    parser.add_argument("-C", "--clobber", dest="clobber", action="store_true",
+                        help="overwrite existing file")
+    args = parser.parse_args()
+
+    setup_pfiles(["acis_bkgrnd_lookup", "dmcopy", "dmmerge", "dmkeypar",
+                  "dmhedit", "reproject_events", "acis_process_events"])
+
+    manifest = get_manifest()
+    if args.infile:
+        infile = args.infile
+    else:
+        infile = manifest.getpath("evt2_clean", relative=True)
+    chips = ACIS.get_chips_str(infile, sep="-")
+    if args.outfile:
+        outfile = args.outfile
+    else:
+        outfile = "blanksky_c{chips}.fits".format(chips=chips)
+    asol = manifest.getpath("asol", relative=True, sep=",")
+    logger.info("infile: %s" % infile)
+    logger.info("outfile: %s" % outfile)
+    logger.info("chips: %s" % chips)
+    logger.info("asol: %s" % asol)
+
+    bkgfiles = find_blanksky(infile, copy=True, clobber=args.clobber)
+    bkg_orig = os.path.splitext(outfile)[0] + "_orig.fits"
+    merge_blanksky(bkgfiles, bkg_orig, clobber=args.clobber)
+    bkg_clean = os.path.splitext(outfile)[0] + "_clean.fits"
+    clean_vfaint(bkg_orig, bkg_clean, clobber=args.clobber)
+    bkg_gained = os.path.splitext(outfile)[0] + "_gained.fits"
+    apply_gainfile(bkg_clean, bkg_gained, evt2=infile, clobber=args.clobber)
+    # Add the PNT header keywords
+    copy_keyword(infile, bkg_gained,
+                 keyword=["RA_PNT", "DEC_PNT", "ROLL_PNT"])
+    reproject_blanksky(bkg_gained, outfile, evt2=infile,
+                       asol=asol, clobber=args.clobber)
+
+    # Add blanksky background to manifest
+    key = "bkg_blank"
+    manifest.setpath(key, outfile)
+    logger.info("Added '%s' to manifest: %s" % (key, manifest.get(key)))
+
+
+if __name__ == "__main__":
+    main()
diff --git a/bin/manifest.py b/bin/manifest.py
new file mode 100755
index 0000000..6b100f0
--- /dev/null
+++ b/bin/manifest.py
@@ -0,0 +1,26 @@
+#!/usr/bin/env python3
+#
+# Copyright (c) 2017 Weitian LI <liweitianux@live.com>
+# MIT license
+#
+# Weitian LI
+# 2017-02-11
+
+"""
+Manage the observation manifest in YAML format.
+
+NOTE
+----
+Use `ruamel.yaml`_ instead of `PyYAML`_ to preserve the comments
+and other structures in the YAML file.
+
+.. _`ruamel.yaml`: https://bitbucket.org/ruamel/yaml
+.. _`PyYAML`: http://pyyaml.org/
+"""
+
+from _context import acispy
+from acispy import manifest
+
+
+if __name__ == "__main__":
+    manifest.main()
diff --git a/bin/renorm_spectrum.py b/bin/renorm_spectrum.py
new file mode 100755
index 0000000..9a0c110
--- /dev/null
+++ b/bin/renorm_spectrum.py
@@ -0,0 +1,75 @@
+#!/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 _context import acispy
+from acispy.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",
+                        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()
diff --git a/bin/results.py b/bin/results.py
new file mode 100755
index 0000000..d325dc3
--- /dev/null
+++ b/bin/results.py
@@ -0,0 +1,18 @@
+#!/usr/bin/env python3
+#
+# Copyright (c) 2017 Weitian LI <liweitianux@live.com>
+# MIT license
+#
+# Weitian LI
+# 2017-02-11
+
+"""
+Manage the analysis results in YAML format.
+"""
+
+from _context import acispy
+from acispy import results
+
+
+if __name__ == "__main__":
+    results.main()