aboutsummaryrefslogtreecommitdiffstats
path: root/scripts
diff options
context:
space:
mode:
authorAaron LI <aaronly.me@outlook.com>2017-02-20 20:13:02 +0800
committerAaron LI <aaronly.me@outlook.com>2017-02-20 20:13:02 +0800
commited49ec62e2ffcfe8af72dd9e860117276c54f0cc (patch)
tree8561858d9845f493f51b5e4188d2f57a216a5b13 /scripts
parente211dd31667551d2123c6788b29d7cea04da3c50 (diff)
downloadchandra-acis-analysis-ed49ec62e2ffcfe8af72dd9e860117276c54f0cc.tar.bz2
Add scripts/make_blanksky.py (to replace ciao_blanksky.sh)
Diffstat (limited to 'scripts')
-rwxr-xr-xscripts/make_blanksky.py228
1 files changed, 228 insertions, 0 deletions
diff --git a/scripts/make_blanksky.py b/scripts/make_blanksky.py
new file mode 100755
index 0000000..8229d39
--- /dev/null
+++ b/scripts/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()