aboutsummaryrefslogtreecommitdiffstats
path: root/bin
diff options
context:
space:
mode:
Diffstat (limited to 'bin')
-rw-r--r--bin/_context.py16
-rwxr-xr-xbin/analyze_path.py52
-rwxr-xr-xbin/calc_centroid.py184
-rwxr-xr-xbin/clean_evt2.py207
-rwxr-xr-xbin/collect_yaml.py74
-rwxr-xr-xbin/cosmo_calc.py133
-rwxr-xr-xbin/event2image.py105
-rwxr-xr-xbin/make_blanksky.py228
-rwxr-xr-xbin/manifest.py26
-rwxr-xr-xbin/renorm_spectrum.py75
-rwxr-xr-xbin/results.py18
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()