path: root/astro/suzaku/make_photon_list.py

#!/usr/bin/env python3
#
# Copyright (c) 2017 Weitian LI <weitian@aaronly.me>
# MIT License
#

"""
Make/simulate the X-ray photon list from the object's image and
spectral models.

The simulated X-ray photon list will be used to simulate the
Suzaku event observation by ``xissim`` tool.

This script is intended to replace and extend the abilities of the
``mkphlist`` tool.

NOTE
----
The environment variable ``HEADAS`` should be set in order to help
locate the ``PyXspec`` module and XSPEC shared libraries.

References
----------
* mkphlist: https://heasarc.gsfc.nasa.gov/lheasoft/ftools/headas/mkphlist.txt
* xissim: https://heasarc.gsfc.nasa.gov/lheasoft/ftools/headas/xissim.txt
* PyXspec: https://heasarc.gsfc.nasa.gov/xanadu/xspec/python/html/index.html

Example Configuration File
-----------------------------------------------------------------------
# image to determine the photon counts distribution
image: imgbox800_e500-7000_sm.fits
# region (annuli below) center; in "image" coordinate
center: [400, 399]
nh: 0.03  # 1e22 [cm^-2]
redshift: 0.0137
# simulated photon energy range [keV]
erange: [0.3, 10.0]
# number of energy bins (logarithmic)
ebins: 1000
# total photon counts that will be generated
counts: 300000
# exposure [ks]
exposure: 50
# a set of annular regions, with several pie regions inside each
# annulus; each pie region can have a different spectral model.
regions:
    # annulus 1, with 3 pies
    - radius: [0, 100]
      angle: [0, 120, 200]
      temperature: [1.0, 1.5, 2.0]
      abundance: [0.5, 1.0, 1.5]
      weight: [1, 2, 1.5]
    # annulus 2, with 3 pies
    - radius: [100, 200]
      angle: [0, 90, 250]
      temperature: [0.5, 1.0, 1.5]
      abundance: [1.5, 2.0, 1.0]
      weight: [0.5, 1, 1.5]
    # annulus 3, with 4 pies
    - radius: [200, 400]
      angle: [50, 150, 220, 300]
      temperature: [0.8, 1.2, 1.5, 1.3]
      abundance: [1.1, 2.0, 1.5, 1.2]
      weight: [0.2, 1.5, 0.7, 2]
clobber: True
outfiles:
    photons_table: photons.fits
    counts_map: counts_map.fits
    temperature_map: temperature_map.fits
    abundance_map: abundance_map.fits
-----------------------------------------------------------------------
"""

import os
import sys

try:
    headas = os.environ["HEADAS"]
    healib = os.path.join(headas, "lib")
except KeyError:
    raise ValueError("env variable 'HEADAS' not set")

if ("LD_LIBRARY_PATH" not in os.environ) or (
        os.environ["LD_LIBRARY_PATH"].find(healib) < 0):
    os.environ["LD_LIBRARY_PATH"] = ":".join([
        healib, os.environ.get("LD_LIBRARY_PATH", "")
    ])
    try:
        # Hack the ``LD_LIBRARY_PATH`` to import Xspec
        # Credit: https://stackoverflow.com/a/25457751/4856091
        print("sys.argv:", sys.argv)
        os.execv(sys.argv[0], sys.argv)
    except Exception:
        print("ERROR: failed to re-exec with new LD_LIBRARY_PATH")
        raise

sys.path.append(os.path.join(healib, "python"))
import xspec
print("Imported XSPEC!")

import argparse
import logging
from pprint import pprint

import yaml
import numpy as np
from astropy.io import fits
from astropy.wcs import WCS


logging.basicConfig(level=logging.INFO,
                    format="[%(levelname)s:%(lineno)d] %(message)s")
logger = logging.getLogger()


class Pie:
    """
    Pie region
    """
    def __init__(self, xc, yc, rin, rout, abegin, aend):
        self.xc = xc
        self.yc = yc
        self.rin = rin
        self.rout = rout
        self.abegin = abegin  # [deg] beginning angle
        self.aend = aend  # [deg] ending angle (may be > 360)
        # spectral model parameters
        self._modelpars = {}

    @staticmethod
    def cart2pol(x, y):
        rho = np.sqrt(x**2 + y**2)
        phi = 180 + np.rad2deg(np.arctan2(y, x))  # 0-360 [deg]
        return (rho, phi)

    def make_mask(self, shape):
        try:
            nrow, ncol = shape
        except TypeError:
            nrow = ncol = shape
        # HACK: to make the masks consistent with ``rand_position()``
        ix = self.xc - np.arange(ncol)
        iy = self.yc - np.arange(nrow)
        mx, my = np.meshgrid(ix, iy)
        rho, phi = self.cart2pol(mx, my)
        mask_rho = (rho >= self.rin) & (rho <= self.rout)
        mask_phi = (phi >= self.abegin) & (phi <= self.aend)
        if self.aend > 360:
            mask_phi |= (phi <= (self.aend-360))
        mask = mask_rho & mask_phi
        return mask

    def rand_position(self, n=None):
        if n is None:
            n = self.modelpar("counts")
        theta = np.random.uniform(low=self.abegin, high=self.aend, size=n)
        r = np.sqrt(np.random.uniform(low=self.rin**2, high=self.rout**2,
                                      size=n))
        x = r * np.cos(np.deg2rad(theta)) + self.xc
        y = r * np.sin(np.deg2rad(theta)) + self.yc
        return (x, y)

    def modelpar(self, key=None, value=None):
        if key is None:
            return self._modelpars
        elif value is None:
            return self._modelpars.get(key)
        else:
            self._modelpars[key] = value

    def set_model(self, nh, redshift):
        model = xspec.Model("wabs*apec")
        model.wabs.nH = nh
        model.apec.Redshift = redshift
        model.apec.kT = self.modelpar("temperature")
        model.apec.Abundanc = self.modelpar("abundance")
        self._model = model

    def rand_photons(self, n=None):
        if n is None:
            n = self.modelpar("counts")
        model = self._model
        mvalues = np.array(model.values(0), dtype=float)  # len: ebins
        p = mvalues / mvalues.sum()
        menergies = np.array(model.energies(0), dtype=float)  # len: ebins+1
        mebins = np.sqrt(menergies[1:] * menergies[:-1])
        photons = np.random.choice(mebins, size=n, p=p)
        return photons  # [keV]


class Regions:
    """
    Configured regions
    """
    def __init__(self, configs):
        self.configs = configs
        self.xc, self.yc = configs["center"]

    @property
    def rmax(self):
        rmax = 0
        for annulus in self.configs["regions"]:
            rin, rout = annulus["radius"]
            if rmax < rout:
                rmax = rout
        return rmax

    def make_mask(self, shape):
        try:
            nrow, ncol = shape
        except TypeError:
            nrow = ncol = shape
        ix = np.arange(ncol) - self.xc
        iy = np.arange(nrow) - self.yc
        mx, my = np.meshgrid(ix, iy)
        rho = np.sqrt(mx**2 + my**2)
        mask = (rho <= self.rmax)
        return mask

    @property
    def regions(self):
        reg_all = []
        for annulus in self.configs["regions"]:
            reg_annulus = []
            rin, rout = annulus["radius"]
            abegin = annulus["angle"]
            aend = abegin[1:] + [abegin[0]+360]
            npie = len(abegin)
            temperature = annulus["temperature"]
            abundance = annulus["abundance"]
            weight = annulus.get("weight", [1]*npie)
            for i in range(npie):
                pie = Pie(xc=self.xc, yc=self.yc, rin=rin, rout=rout,
                          abegin=abegin[i], aend=aend[i])
                pie.modelpar("temperature", temperature[i])
                pie.modelpar("abundance", abundance[i])
                pie.modelpar("weight", weight[i])
                reg_annulus.append(pie)
            reg_all.append(reg_annulus)
        return reg_all


def pixel2world(x, y, wcs):
    pix = np.column_stack([x, y])
    world = wcs.wcs_pix2world(pix, 0)
    ra = world[:, 0]
    dec = world[:, 1]
    return (ra, dec)  # [deg]


def main():
    parser = argparse.ArgumentParser(
        description="Make/simulate X-ray photon list for Suzaku simulation")
    parser.add_argument("config", help="configuration file in YAML format")
    args = parser.parse_args()

    configs = yaml.load(open(args.config))
    logger.info("Load configuration file: %s" % args.config)
    logger.info("Configurations:")
    pprint(configs)

    # Update XSPEC settings
    emin, emax = configs["erange"]  # [keV]
    ebins = configs["ebins"]
    xspec.AllModels.setEnergies("%.1f %.1f %d log" % (emin, emax, ebins))
    logger.info("Energy range: [%.1f, %.1f] [keV]" % (emin, emax))
    logger.info("Energy: %d logarithmic channels" % ebins)

    with fits.open(configs["image"]) as f:
        header = f[0].header
        image = f[0].data
    shape = image.shape
    logger.info("Image size: %dx%d" % (shape[1], shape[0]))

    wcs = WCS(header)
    regions = Regions(configs)
    reg_all = regions.regions
    mask_all = regions.make_mask(shape=shape)
    weight_all = np.sum(image[mask_all])

    counts_all = configs["counts"]
    logger.info("Total counts: %d" % counts_all)

    logger.info("nH: %.4f [1e22 cm^-2]" % configs["nh"])
    logger.info("Redshift: %.5f" % configs["redshift"])
    exposure = configs["exposure"] * 1e3  # [s]
    logger.info("Exposure time: %.1f [s]" % exposure)

    logger.info("Determining photons counts in each region ...")
    counts_sum = 0
    for i, annulus in enumerate(reg_all):
        for j, pie in enumerate(annulus):
            label = "annu#%d/pie#%d" % (i+1, j+1)
            mask = pie.make_mask(shape=shape)
            pixels = np.sum(mask)
            weight = np.sum(image[mask]) * pie.modelpar("weight")
            counts = int(counts_all * weight / weight_all)
            counts_sum += counts
            pie.modelpar("pixels", pixels)
            pie.modelpar("counts", counts)
            logger.info("%s: %d pixels, %d photons" % (label, pixels, counts))

    logger.info("Determined counts sum: %d" % counts_sum)
    logger.info("Adjusting total counts -> %d" % counts_all)
    for i, annulus in enumerate(reg_all):
        for j, pie in enumerate(annulus):
            label = "annu#%d/pie#%d" % (i+1, j+1)
            counts_old = pie.modelpar("counts")
            counts_new = round(counts_old * counts_all / counts_sum)
            pie.modelpar("counts", counts_new)
            logger.info("%s: adjusted photon counts: %d -> %d" %
                        (label, counts_old, counts_new))

    # Output files
    temp_map = np.zeros_like(image)
    abund_map = np.zeros_like(image)
    counts_map = np.zeros_like(image)
    weights_map = np.zeros_like(image)
    photonlist = []

    for i, annulus in enumerate(reg_all):
        for j, pie in enumerate(annulus):
            label = "annu#%d/pie#%d" % (i+1, j+1)
            pie.set_model(nh=configs["nh"], redshift=configs["redshift"])
            mask = pie.make_mask(shape=shape)
            temp = pie.modelpar("temperature")
            abund = pie.modelpar("abundance")
            counts = pie.modelpar("counts")
            logger.info("%s: kT=%.2f, Z=%.2f, %d photons" %
                        (label, temp, abund, counts))

            logger.info("%s: sampling photon positions ..." % label)
            x, y = pie.rand_position(n=counts)
            ra, dec = pixel2world(x, y, wcs=wcs)
            logger.info("%s: sampling photon energies ..." % label)
            energies = pie.rand_photons(n=counts)
            time = np.random.uniform(low=0, high=exposure, size=counts)
            photons = np.column_stack([time, energies, ra, dec])
            photonlist.append(photons)

            logger.info("%s: spatially binning photons ..." % label)
            rbins = np.arange(shape[0]+1, dtype=int)
            cbins = np.arange(shape[1]+1, dtype=int)
            hist2d, __, __ = np.histogram2d(y, x, bins=(rbins, cbins))
            counts_map += hist2d

            temp_map[mask] = temp
            abund_map[mask] = abund
            weights_map[mask] = pie.modelpar("weight")

    logger.info("Creating output FITS header ...")
    header_out = fits.Header()
    header_out.extend(wcs.to_header(), update=True)
    header_out["CREATOR"] = os.path.basename(sys.argv[0])
    header_out.add_history(" ".join(sys.argv))
    logger.info("Creating photons table ...")
    photons = np.row_stack(photonlist)
    photons = photons[photons[:, 0].argsort()]  # sort by time (1st column)
    hdu = fits.BinTableHDU.from_columns([
        fits.Column(name="PHOTON_TIME", format="D", unit="s",
                    array=photons[:, 0]),
        fits.Column(name="PHOTON_ENERGY", format="E", unit="keV",
                    array=photons[:, 1]),
        fits.Column(name="RA", format="E", unit="deg", array=photons[:, 2]),
        fits.Column(name="DEC", format="E", unit="deg", array=photons[:, 3]),
    ], header=header_out)
    hdu.name = "PHOTON_LIST"
    outfile = configs["outfiles"]["photons_table"]
    hdu.writeto(outfile, overwrite=configs["clobber"])
    logger.info("Wrote photons table to: %s" % outfile)

    data = np.stack([counts_map, weights_map], axis=0)
    hdu = fits.PrimaryHDU(data=data, header=header_out)
    outfile = configs["outfiles"]["counts_map"]
    hdu.writeto(outfile, overwrite=configs["clobber"])
    logger.info("Wrote counts/weights map to: %s" % outfile)
    #
    hdu = fits.PrimaryHDU(data=temp_map, header=header_out)
    outfile = configs["outfiles"]["temperature_map"]
    hdu.writeto(outfile, overwrite=configs["clobber"])
    logger.info("Wrote temperature map to: %s" % outfile)
    #
    hdu = fits.PrimaryHDU(data=abund_map, header=header_out)
    outfile = configs["outfiles"]["abundance_map"]
    hdu.writeto(outfile, overwrite=configs["clobber"])
    logger.info("Wrote abundance map to: %s" % outfile)


if __name__ == "__main__":
    main()