Add make_expmap.py: Make spectral-weighted exposure map

1 files changed, 247 insertions, 0 deletions

diff --git a/bin/make_expmap.py b/bin/make_expmap.py
new file mode 100755
index 0000000..ff1fbb0
--- /dev/null
+++ b/bin/make_expmap.py
@@ -0,0 +1,247 @@
+#!/usr/bin/env python3
+#
+# Copyright (c) 2017 Weitian LI <liweitianux@live.com>
+# MIT license
+
+"""
+Create the spectral-weighted exposure map, which will be used to produce
+the exposure-corrected images and to extract the surface brightness profiles.
+
+NOTE
+----
+We do not use the CIAO tool ``fluximage`` to create the exposure map,
+because it requires an event file (other than an image) as input file.
+However, we want to create the exposure map directly based on the
+(previously created) image file.
+(But an event file is required to create the aspect histograms for each
+chip.)
+
+NOTE/CAVEAT
+-----------
+To compute an image which gives the integrated flux over the full energy
+range, it may be best to first compute flux-corrected images in several
+narrow energy bands (where the ARF is nearly flat) and then sum those
+fluxed images together.  Weighted exposure maps work well for an energy
+band where the ARF variation isn't very large, but for a full-band
+0.5-10 keV image, it may not be a good idea to compute the flux by
+dividing the counts image by a single number. This is especially true
+for cases where the source spectrum varies significantly within the image;
+in that case, there is no general way to compute a single set of weights
+which will be sensible for every part of the image.
+
+References
+----------
+* Chandra Memo: An Introduction to Exposure Map
+  http://cxc.harvard.edu/ciao/download/doc/expmap_intro.ps
+* CIAO: Single Chip ACIS Exposure Map and Exposure-corrected Image
+  http://cxc.harvard.edu/ciao/threads/expmap_acis_single/
+* CIAO: Multiple Chip ACIS Exposure Map and Exposure-corrected Image
+  http://cxc.harvard.edu/ciao/threads/expmap_acis_multi/
+* CIAO: Calculating Spectral Weights for mkinstmap
+  http://cxc.harvard.edu/ciao/threads/spectral_weights/
+"""
+
+import os
+import argparse
+import subprocess
+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
+
+
+logging.basicConfig(level=logging.INFO)
+logger = logging.getLogger(__name__)
+
+
+def get_xygrid(image):
+    """
+    Get the ``xygrid`` of the input image for later ``mkexpmap``'s use.
+    """
+    logger.info("Get the 'xygrid' of image: %s" % image)
+    subprocess.check_call(["punlearn", "get_sky_limits"])
+    subprocess.check_call([
+        "get_sky_limits", "image=%s" % image, "verbose=0"
+    ])
+    xygrid = subprocess.check_output([
+        "pget", "get_sky_limits", "xygrid"
+    ]).decode("utf-8").strip()
+    logger.info("xygrid: %s" % xygrid)
+    return xygrid
+
+
+def make_aspect_histogram(outfile, asol, evtfile, chip, clobber=False):
+    """
+    Create the aspect histogram for each chip, detailing the aspect
+    history of the observation.
+    """
+    clobber = "yes" if clobber else "no"
+    logger.info("Make aspect histogram for chip: %s" % chip)
+    subprocess.check_call(["punlearn", "asphist"])
+    subprocess.check_call([
+        "asphist", "infile=%s" % asol, "outfile=%s" % outfile,
+        "evtfile=%s[ccd_id=%s]" % (evtfile, chip),
+        "clobber=%s" % clobber
+    ])
+
+
+def make_instrument_map(outfile, spectrumfile, chip, obsfile,
+                        maskfile, badpixfile, clobber=False):
+    """
+    Make the spectral-weighted instrument map (effective area vs.
+    detector position).
+
+    NOTE
+    ----
+    The ``obsfile`` should be a FITS file containing keywords which specify
+    the mission, detector, SIM offsets, observation date, etc.  This should
+    be an event file; can also be a binned image file (tested OK).
+    """
+    logger.info("Set bad pixel file: %s" % badpixfile)
+    subprocess.check_call(["punlearn", "ardlib"])
+    subprocess.check_call(["punlearn", "acis_set_ardlib"])
+    subprocess.check_call([
+        "acis_set_ardlib", "badpixfile=%s" % badpixfile, "verbose=0"
+    ])
+
+    clobber = "yes" if clobber else "no"
+    pixelgrid = "1:1024:#1024,1:1024:#1024"
+    logger.info("Make spectral-weighted instrument map for chip: %s" % chip)
+    subprocess.check_call(["punlearn", "mkinstmap"])
+    subprocess.check_call([
+        "mkinstmap", "outfile=%s" % outfile,
+        "spectrumfile=%s" % spectrumfile,
+        "pixelgrid=%s" % pixelgrid,
+        "detsubsys=ACIS-%s" % chip,
+        "obsfile=%s" % obsfile,
+        "maskfile=%s" % maskfile,
+        "clobber=%s" % clobber, "mode=h"
+    ])
+
+
+def make_exposure_map(outfile, asphistfile, instmapfile, xygrid,
+                      clobber=False):
+    """
+    Create the spectral-weighted exposure map by projecting the instrument
+    map onto the sky with the aspect information.
+    """
+    clobber = "yes" if clobber else "no"
+    logger.info("Make exposure map by projecting instrument map ...")
+    subprocess.check_call(["punlearn", "mkexpmap"])
+    subprocess.check_call([
+        "mkexpmap", "outfile=%s" % outfile,
+        "asphistfile=%s" % asphistfile,
+        "instmapfile=%s" % instmapfile,
+        "xygrid=%s" % xygrid,
+        "normalize=no", "useavgaspect=no",
+        "clobber=%s" % clobber, "mode=h"
+    ])
+
+
+def combine_expmaps(outfile, expmaps, clobber=False):
+    """
+    Combine multiple exposure maps of each chip into a single one
+    for ACIS-I (i.e., chips 0123).
+    """
+    if len(expmaps) > 1:
+        logger.info("Combine exposure maps: {0}".format(", ".join(expmaps)))
+        operation = "imgout=%s" % "+".join([
+            "img%d" % (i+1) for i in range(len(expmaps))
+        ])
+        clobber = "yes" if clobber else "no"
+        subprocess.check_call(["punlearn", "dmimgcalc"])
+        subprocess.check_call([
+            "dmimgcalc", "infile=%s" % ",".join(expmaps),
+            "infile2=none", "outfile=%s" % outfile,
+            "operation=%s" % operation,
+            "clobber=%s" % clobber
+        ])
+        for f in expmaps:
+            os.remove(f)
+    else:
+        logger.info("No need to combine exposure maps.")
+        if os.path.exists(outfile) and (not clobber):
+            raise OSError("File already exists: %s" % outfile)
+        os.rename(expmaps[0], outfile)
+
+
+def main():
+    parser = argparse.ArgumentParser(
+        description="Make spectral-weighted exposure map")
+    parser.add_argument("-C", "--clobber", dest="clobber", action="store_true",
+                        help="overwrite existing file")
+    parser.add_argument("-w", "--weights", dest="weights",
+                        help="spectral weights file (default: " +
+                        "'spec_weights' from manifest)")
+    parser.add_argument("-e", "--evtfile", dest="evtfile",
+                        help="event file for aspect histogram creation " +
+                        "(default: 'evt2_clean' from manifest)")
+    parser.add_argument("-i", "--infile", dest="infile", required=True,
+                        help="input image file")
+    parser.add_argument("-o", "--outfile", dest="outfile", required=True,
+                        help="filename of output exposure map")
+    args = parser.parse_args()
+
+    setup_pfiles(["get_sky_limits", "asphist", "mkinstmap", "ardlib",
+                  "acis_set_ardlib", "mkexpmap", "dmimgcalc"])
+
+    manifest = get_manifest()
+    if args.weights:
+        weights = args.weights
+    else:
+        weights = manifest.getpath("spec_weights", relative=True)
+    if args.evtfile:
+        evtfile = args.evtfile
+    else:
+        evtfile = manifest.getpath("evt2_clean", relative=True)
+    asol = manifest.getpath("asol", relative=True, sep=",")
+    bpix = manifest.getpath("bpix", relative=True)
+    msk = manifest.getpath("msk", relative=True)
+    chips = ACIS.get_chips_str(args.infile)
+    logger.info("infile: %s" % args.infile)
+    logger.info("output expmap: %s" % args.outfile)
+    logger.info("weights: %s" % weights)
+    logger.info("evtfile: %s" % evtfile)
+    logger.info("chips: %s" % chips)
+    logger.info("bpix: %s" % bpix)
+    logger.info("asol: %s" % asol)
+    logger.info("msk: %s" % msk)
+
+    xygrid = get_xygrid(args.infile)
+    expmaps = []
+
+    for c in chips:
+        logger.info("Processing chip %s ..." % c)
+        asphist = "asphist_c{chip}.fits".format(chip=c)
+        instmap = "instmap_c{chip}.fits".format(chip=c)
+        expmap = "expmap_c{chip}.fits".format(chip=c)
+        make_aspect_histogram(outfile=asphist, asol=asol, evtfile=evtfile,
+                              chip=c, clobber=args.clobber)
+        make_instrument_map(outfile=instmap, spectrumfile=weights,
+                            chip=c, obsfile=args.infile, maskfile=msk,
+                            badpixfile=bpix, clobber=args.clobber)
+        make_exposure_map(outfile=expmap, asphistfile=asphist,
+                          instmapfile=instmap, xygrid=xygrid,
+                          clobber=args.clobber)
+        expmaps.append(expmap)
+        # Remove intermediate files
+        os.remove(asphist)
+        os.remove(instmap)
+
+    combine_expmaps(outfile=args.outfile, expmaps=expmaps,
+                    clobber=args.clobber)
+    detnam = "ACIS-{0}".format(chips)
+    logger.info("Update keyword 'DETNAM' to %s" % detnam)
+    write_keyword(args.outfile, keyword="DETNAM", value=detnam)
+
+    # Add created exposure map and exposure-corrected image to manifest
+    key = "expmap"
+    manifest.setpath(key, args.outfile)
+    logger.info("Added '%s' to manifest: %s" % (key, manifest.get(key)))
+
+
+if __name__ == "__main__":
+    main()