diff options
Diffstat (limited to 'scripts')
| -rw-r--r-- | scripts/_context.py | 16 | ||||
| -rwxr-xr-x | scripts/analyze_path.py | 52 | ||||
| -rwxr-xr-x | scripts/calc_centroid.py | 184 | ||||
| -rwxr-xr-x | scripts/clean_evt2.py | 207 | ||||
| -rwxr-xr-x | scripts/collect_yaml.py | 74 | ||||
| -rwxr-xr-x | scripts/cosmo_calc.py | 133 | ||||
| -rwxr-xr-x | scripts/event2image.py | 105 | ||||
| -rwxr-xr-x | scripts/make_blanksky.py | 228 | ||||
| -rwxr-xr-x | scripts/manifest.py | 26 | ||||
| -rwxr-xr-x | scripts/renorm_spectrum.py | 75 | ||||
| -rwxr-xr-x | scripts/results.py | 18 |
11 files changed, 0 insertions, 1118 deletions
diff --git a/scripts/_context.py b/scripts/_context.py deleted file mode 100644 index b23b2c8..0000000 --- a/scripts/_context.py +++ /dev/null @@ -1,16 +0,0 @@ -# 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/scripts/analyze_path.py b/scripts/analyze_path.py deleted file mode 100755 index d2acd42..0000000 --- a/scripts/analyze_path.py +++ /dev/null @@ -1,52 +0,0 @@ -#!/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/scripts/calc_centroid.py b/scripts/calc_centroid.py deleted file mode 100755 index 7df71b1..0000000 --- a/scripts/calc_centroid.py +++ /dev/null @@ -1,184 +0,0 @@ -#!/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/scripts/clean_evt2.py b/scripts/clean_evt2.py deleted file mode 100755 index 5b53922..0000000 --- a/scripts/clean_evt2.py +++ /dev/null @@ -1,207 +0,0 @@ -#!/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/scripts/collect_yaml.py b/scripts/collect_yaml.py deleted file mode 100755 index 3147a0c..0000000 --- a/scripts/collect_yaml.py +++ /dev/null @@ -1,74 +0,0 @@ -#!/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/scripts/cosmo_calc.py b/scripts/cosmo_calc.py deleted file mode 100755 index 1302bdb..0000000 --- a/scripts/cosmo_calc.py +++ /dev/null @@ -1,133 +0,0 @@ -#!/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/scripts/event2image.py b/scripts/event2image.py deleted file mode 100755 index 44189e0..0000000 --- a/scripts/event2image.py +++ /dev/null @@ -1,105 +0,0 @@ -#!/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/scripts/make_blanksky.py b/scripts/make_blanksky.py deleted file mode 100755 index 8229d39..0000000 --- a/scripts/make_blanksky.py +++ /dev/null @@ -1,228 +0,0 @@ -#!/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/scripts/manifest.py b/scripts/manifest.py deleted file mode 100755 index 6b100f0..0000000 --- a/scripts/manifest.py +++ /dev/null @@ -1,26 +0,0 @@ -#!/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/scripts/renorm_spectrum.py b/scripts/renorm_spectrum.py deleted file mode 100755 index 9a0c110..0000000 --- a/scripts/renorm_spectrum.py +++ /dev/null @@ -1,75 +0,0 @@ -#!/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/scripts/results.py b/scripts/results.py deleted file mode 100755 index d325dc3..0000000 --- a/scripts/results.py +++ /dev/null @@ -1,18 +0,0 @@ -#!/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() |