path: root/astro/spectrum/crosstalk_deprojection.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
#
# References:
# [1] Definition of RMF and ARF file formats
#     https://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/docs/memos/cal_gen_92_002/cal_gen_92_002.html
# [2] The OGIP Spectral File Format
#     https://heasarc.gsfc.nasa.gov/docs/heasarc/ofwg/docs/summary/ogip_92_007_summary.html
# [3] CIAO: Auxiliary Response File
#     http://cxc.harvard.edu/ciao/dictionary/arf.html
# [4] CIAO: Redistribution Matrix File
#     http://cxc.harvard.edu/ciao/dictionary/rmf.html
# [5] astropy - FITS format code
#     http://docs.astropy.org/en/stable/io/fits/usage/table.html#column-creation
# [6] XSPEC - Spectral Fitting
#     https://heasarc.gsfc.nasa.gov/docs/xanadu/xspec/manual/XspecSpectralFitting.html
# [7] Direct X-ray Spectra Deprojection
#     https://www-xray.ast.cam.ac.uk/papers/dsdeproj/
#     Sanders & Fabian 2007, MNRAS, 381, 1381
#
#
# Weitian LI
# Created: 2016-03-26
# Updated: 2016-06-07
#
# Change log:
# 2016-06-07:
#   * Explain the errors/uncertainties calculation approach
# 2016-04-20:
#   * Add argument 'add_history' to some methods (to avoid many duplicated
#     histories due to Monte Carlo)
#   * Rename 'reset_header_keywords()' to 'fix_header_keywords()',
#     and add mandatory spectral keywords if missing
#   * Add method 'fix_header()' to class 'Crosstalk' and 'Deprojection',
#     and fix the headers before write spectra
# 2016-04-19:
#   * Ignore numpy error due to division by zero
#   * Update tool description and sample configuration
#   * Add two other main methods: `main_deprojection()' and `main_crosstalk()'
#   * Add argument 'group_squeeze' to some methods for better performance
#   * Rename from 'correct_crosstalk.py' to 'crosstalk_deprojection.py'
# 2016-04-18:
#   * Implement deprojection function: class Deprojection
#   * Support spectral grouping (supply the grouping specification)
#   * Add grouping, estimate_errors, copy, randomize, etc. methods
#   * Utilize the Monte Carlo techniques to estimate the final spectral errors
#   * Collect all ARFs and RMFs within dictionaries
# 2016-04-06:
#   * Fix `RMF: get_rmfimg()' for XMM EPIC RMF
# 2016-04-02:
#   * Interpolate ARF in order to match the spectral channel energies
#   * Add version and date information
#   * Update documentations
#   * Update header history contents
# 2016-04-01:
#   * Greatly update the documentations (e.g., description, sample config)
#   * Add class `RMF'
#   * Add method `get_energy()' for class `ARF'
#   * Split out class `SpectrumSet' from `Spectrum'
#   * Implement background subtraction
#   * Add config `subtract_bkg' and corresponding argument
#
# XXX/FIXME:
#   * Deprojection: account for ARF differences across different regions
#
# TODO:
#   * Split classes ARF, RMF, Spectrum, and SpectrumSet to a separate module
#

__version__ = "0.5.3"
__date__    = "2016-06-07"


"""
Correct the crosstalk effect of XMM spectra by subtracting the photons
that scattered from the surrounding regions due to the finite PSF, and
by compensating the photons that scattered to the surrounding regions,
according to the generated crosstalk ARFs by SAS `arfgen'.

After the crosstalk effect being corrected, the deprojection is performed
to deproject the crosstalk-corrected spectra to derive the spectra with
both the crosstalk effect and projection effect corrected.


Sample config file (in `ConfigObj' syntax):
-----------------------------------------------------------
# operation mode: deprojection, crosstalk, or both (default)
mode = both
# supply a *groupped* spectrum (from which the "GROUPING" and "QUALITY"
# are used to group all the following spectra)
grouping = spec_grp.pi
# whether to subtract the background before crosstalk correction
subtract_bkg = True
# whether to fix the negative channel values due to spectral subtractions
fix_negative = False
# Monte Carlo times for spectral error estimation
mc_times = 5000
# show progress details and verbose information
verbose = True
# overwrite existing files
clobber = False

# NOTE:
# ONLY specifiy ONE set of projected spectra (i.e., from the same detector
# of one observation), since ALL the following specified spectra will be
# used for the deprojection.

[reg1]
...

[reg2]
outfile = deprojcc_reg2.pi
spec = reg2.pi
arf = reg2.arf
rmf = reg2.rmf
bkg = reg2_bkg.pi
  [[cross_in]]
    [[[in1]]]
    spec = reg1.pi
    arf = reg1.arf
    rmf = reg1.rmf
    bkg = reg1_bkg.pi
    cross_arf = reg_1-2.arf
    [[[in2]]]
    spec = reg3.pi
    arf = reg3.arf
    rmf = reg3.rmf
    bkg = reg3_bkg.pi
    cross_arf = reg_3-2.arf
  [[cross_out]]
  cross_arf = reg_2-1.arf, reg_2-3.arf

[...]
...
-----------------------------------------------------------
"""

WARNING = """
********************************* WARNING ************************************
The generated spectra are substantially modified (e.g., scale, add, subtract),
therefore, take special care when interpretating the fitting results,
especially the metal abundances and normalizations.
******************************************************************************
"""


import sys
import os
import argparse
from datetime import datetime
from copy import copy

import numpy as np
import scipy as sp
import scipy.interpolate
from astropy.io import fits
from configobj import ConfigObj


def group_data(data, grouping):
    """
    Group the data with respect to the supplied `grouping' specification
    (i.e., "GROUPING" columns of a spectrum).  The channel counts of the
    same group are summed up and assigned to the FIRST channel of this
    group, while the OTHRE channels are all set to ZERO.
    """
    data_grp = np.array(data).copy()
    for i in reversed(range(len(data))):
        if grouping[i] == 1:
            # the beginning channel of a group
            continue
        else:
            # other channels of a group
            data_grp[i-1] += data_grp[i]
            data_grp[i]    = 0
    assert np.isclose(sum(data_grp), sum(data))
    return data_grp


class ARF:  # {{{
    """
    Class to handle the ARF (ancillary/auxiliary response file),
    which contains the combined instrumental effective area
    (telescope/filter/detector) and the quantum efficiency (QE) as a
    function of energy averaged over time.
    The effective area is [cm^2], and the QE is [counts/photon]; they are
    multiplied together to create the ARF, resulting in [cm^2 counts/photon].

    **CAVEAT/NOTE**:
    Generally, the "ENERG_LO" and "ENERG_HI" columns of an ARF are *different*
    to the "E_MIN" and "E_MAX" columns of a RMF (which are corresponding
    to the spectrum channel energies).
    For the XMM EPIC *pn* and Chandra *ACIS*, the generated ARF does NOT have
    the same number of data points to that of spectral channels, i.e., the
    "ENERG_LO" and "ENERG_HI" columns of ARF is different to the "E_MIN" and
    "E_MAX" columns of RMF.
    Therefore it is necessary to interpolate and extrapolate the ARF curve
    in order to match the spectrum (or RMF "EBOUNDS" extension).
    As for the XMM EPIC *MOS1* and *MOS2*, the ARF data points match the
    spectral channels, i.e., the energy positions of each ARF data point and
    spectral channel are consistent.  Thus the interpolation is not needed.

    References:
    [1] CIAO: Auxiliary Response File
        http://cxc.harvard.edu/ciao/dictionary/arf.html
    [2] Definition of RMF and ARF file formats
        https://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/docs/memos/cal_gen_92_002/cal_gen_92_002.html
    """
    filename = None
    fitsobj  = None
    # only consider the "SPECTRUM" extension
    header   = None
    energ_lo = None
    energ_hi = None
    specresp = None
    # function of the interpolated ARF
    f_interp = None
    # energies of the spectral channels
    energy_channel = None
    # spectral channel grouping specification
    grouping       = None
    groupped       = False
    # groupped ARF channels with respect to the grouping
    specresp_grp   = None

    def __init__(self, filename):
        self.filename = filename
        self.fitsobj  = fits.open(filename)
        ext_specresp  = self.fitsobj["SPECRESP"]
        self.header   = ext_specresp.header
        self.energ_lo = ext_specresp.data["ENERG_LO"]
        self.energ_hi = ext_specresp.data["ENERG_HI"]
        self.specresp = ext_specresp.data["SPECRESP"]

    def get_data(self, groupped=False, group_squeeze=False, copy=True):
        if groupped:
            specresp = self.specresp_grp
            if group_squeeze:
                specresp = specresp[self.grouping == 1]
        else:
            specresp = self.specresp
        if copy:
            return specresp.copy()
        else:
            return specresp

    def get_energy(self, mean="geometric"):
        """
        Return the mean energy values of the ARF.

        Arguments:
          * mean: type of the mean energy:
                  + "geometric": geometric mean, i.e., e = sqrt(e_min*e_max)
                  + "arithmetic": arithmetic mean, i.e., e = 0.5*(e_min+e_max)
        """
        if mean == "geometric":
            energy = np.sqrt(self.energ_lo * self.energ_hi)
        elif mean == "arithmetic":
            energy = 0.5 * (self.energ_lo + self.energ_hi)
        else:
            raise ValueError("Invalid mean type: %s" % mean)
        return energy

    def interpolate(self, x=None, verbose=False):
        """
        Cubic interpolate the ARF curve using `scipy.interpolate'

        If the requested point is outside of the data range, the
        fill value of *zero* is returned.

        Arguments:
          * x: points at which the interpolation to be calculated.

        Return:
          If x is None, then the interpolated function is returned,
          otherwise, the interpolated data are returned.
        """
        if not hasattr(self, "f_interp") or self.f_interp is None:
            energy = self.get_energy()
            arf    = self.get_data(copy=False)
            if verbose:
                print("INFO: interpolating '%s' (this may take a while) ..." \
                        % self.filename, file=sys.stderr)
            f_interp = sp.interpolate.interp1d(energy, arf, kind="cubic",
                    bounds_error=False, fill_value=0.0, assume_sorted=True)
            self.f_interp = f_interp
        if x is not None:
            return self.f_interp(x)
        else:
            return self.f_interp

    def apply_grouping(self, energy_channel, grouping, verbose=False):
        """
        Group the ARF channels (INTERPOLATED with respect to the spectral
        channels) by the supplied grouping specification.

        Arguments:
          * energy_channel: energies of the spectral channel
          * grouping: spectral grouping specification

        Return: `self.specresp_grp'
        """
        if self.groupped:
            return
        if verbose:
            print("INFO: Grouping spectrum '%s' ..." % self.filename,
                    file=sys.stderr)
        self.energy_channel = energy_channel
        self.grouping = grouping
        # interpolate the ARF w.r.t the spectral channel energies
        arf_interp = self.interpolate(x=energy_channel, verbose=verbose)
        self.specresp_grp = group_data(arf_interp, grouping)
        self.groupped = True
# class ARF }}}


class RMF:  # {{{
    """
    Class to handle the RMF (redistribution matrix file),
    which maps from energy space into detector pulse height (or position)
    space.  Since detectors are not perfect, this involves a spreading of
    the observed counts by the detector resolution, which is expressed as
    a matrix multiplication.
    For X-ray spectral analysis, the RMF encodes the probability R(E,p)
    that a detected photon of energy E will be assisgned to a given
    channel value (PHA or PI) of p.

    The standard Legacy format [2] for the RMF uses a binary table in which
    each row contains R(E,p) for a single value of E as a function of p.
    Non-zero sequences of elements of R(E,p) are encoded using a set of
    variable length array columns.  This format is compact but hard to
    manipulate and understand.

    **CAVEAT/NOTE**:
    + See also the above ARF CAVEAT/NOTE.
    + The "EBOUNDS" extension contains the `CHANNEL', `E_MIN' and `E_MAX'
      columns.  This `CHANNEL' is the same as that of a spectrum.  Therefore,
      the energy values determined from the `E_MIN' and `E_MAX' columns are
      used to interpolate and extrapolate the ARF curve.
    + The `ENERG_LO' and `ENERG_HI' columns of the "MATRIX" extension are
      the same as that of a ARF.

    References:
    [1] CIAO: Redistribution Matrix File
        http://cxc.harvard.edu/ciao/dictionary/rmf.html
    [2] Definition of RMF and ARF file formats
        https://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/docs/memos/cal_gen_92_002/cal_gen_92_002.html
    """
    filename    = None
    fitsobj     = None
    ## extension "MATRIX"
    hdr_matrix  = None
    energ_lo    = None
    energ_hi    = None
    n_grp       = None
    f_chan      = None
    n_chan      = None
    # raw squeezed RMF matrix data
    matrix      = None
    ## extension "EBOUNDS"
    hdr_ebounds = None
    channel     = None
    e_min       = None
    e_max       = None
    ## converted 2D RMF matrix/image from the squeezed binary table
    #  size: len(energ_lo) x len(channel)
    rmfimg      = None

    def __init__(self, filename):
        self.filename    = filename
        self.fitsobj     = fits.open(filename)
        ## "MATRIX" extension
        ext_matrix       = self.fitsobj["MATRIX"]
        self.hdr_matrix  = ext_matrix.header
        self.energ_lo    = ext_matrix.data["ENERG_LO"]
        self.energ_hi    = ext_matrix.data["ENERG_HI"]
        self.n_grp       = ext_matrix.data["N_GRP"]
        self.f_chan      = ext_matrix.data["F_CHAN"]
        self.n_chan      = ext_matrix.data["N_CHAN"]
        self.matrix      = ext_matrix.data["MATRIX"]
        ## "EBOUNDS" extension
        ext_ebounds      = self.fitsobj["EBOUNDS"]
        self.hdr_ebounds = ext_ebounds.header
        self.channel     = ext_ebounds.data["CHANNEL"]
        self.e_min       = ext_ebounds.data["E_MIN"]
        self.e_max       = ext_ebounds.data["E_MAX"]

    def get_energy(self, mean="geometric"):
        """
        Return the mean energy values of the RMF "EBOUNDS".

        Arguments:
          * mean: type of the mean energy:
                  + "geometric": geometric mean, i.e., e = sqrt(e_min*e_max)
                  + "arithmetic": arithmetic mean, i.e., e = 0.5*(e_min+e_max)
        """
        if mean == "geometric":
            energy = np.sqrt(self.e_min * self.e_max)
        elif mean == "arithmetic":
            energy = 0.5 * (self.e_min + self.e_max)
        else:
            raise ValueError("Invalid mean type: %s" % mean)
        return energy

    def get_rmfimg(self):
        """
        Convert the RMF data in squeezed binary table (standard Legacy format)
        to a 2D image/matrix.
        """
        def _make_rmfimg_row(n_channel, dtype, f_chan, n_chan, mat_row):
            # make sure that `f_chan' and `n_chan' are 1-D numpy array
            f_chan = np.array(f_chan).reshape(-1)
            f_chan -= 1  # FITS indices are 1-based
            n_chan = np.array(n_chan).reshape(-1)
            idx = np.concatenate([ np.arange(f, f+n) \
                    for f, n in zip(f_chan, n_chan) ])
            rmfrow = np.zeros(n_channel, dtype=dtype)
            rmfrow[idx] = mat_row
            return rmfrow
        #
        if self.rmfimg is None:
            # Make the 2D RMF matrix/image
            n_energy  = len(self.energ_lo)
            n_channel = len(self.channel)
            rmf_dtype = self.matrix[0].dtype
            rmfimg = np.zeros(shape=(n_energy, n_channel), dtype=rmf_dtype)
            for i in np.arange(n_energy)[self.n_grp > 0]:
                rmfimg[i, :] = _make_rmfimg_row(n_channel, rmf_dtype,
                        self.f_chan[i], self.n_chan[i], self.matrix[i])
            self.rmfimg = rmfimg
        return self.rmfimg

    def write_rmfimg(self, outfile, clobber=False):
        rmfimg = self.get_rmfimg()
        # merge headers
        header = self.hdr_matrix.copy(strip=True)
        header.extend(self.hdr_ebounds.copy(strip=True))
        outfits = fits.PrimaryHDU(data=rmfimg, header=header)
        outfits.writeto(outfile, checksum=True, clobber=clobber)
# class RMF }}}


class Spectrum:  # {{{
    """
    Class that deals with the X-ray spectrum file (usually *.pi).
    """
    filename  = None
    # FITS object return by `fits.open()'
    fitsobj   = None
    # header of "SPECTRUM" extension
    header    = None
    # "SPECTRUM" extension data
    channel   = None
    # name of the spectrum data column (i.e., type, "COUNTS" or "RATE")
    spec_type = None
    # unit of the spectrum data ("count" for "COUNTS", "count/s" for "RATE")
    spec_unit = None
    # spectrum data
    spec_data = None
    # estimated spectral errors for each channel/group
    spec_err  = None
    # statistical errors for each channel/group
    stat_err  = None
    # grouping and quality
    grouping  = None
    quality   = None
    # whether the spectral data being groupped
    groupped  = False
    # several important keywords
    EXPOSURE  = None
    BACKSCAL  = None
    AREASCAL  = None
    RESPFILE  = None
    ANCRFILE  = None
    BACKFILE  = None
    # numpy dtype and FITS format code of the spectrum data
    spec_dtype       = None
    spec_fits_format = None
    # output filename for writing the spectrum if no filename provided
    outfile = None

    def __init__(self, filename, outfile=None):
        self.filename = filename
        self.fitsobj  = fits.open(filename)
        ext_spec      = self.fitsobj["SPECTRUM"]
        self.header   = ext_spec.header.copy(strip=True)
        colnames      = ext_spec.columns.names
        if "COUNTS" in colnames:
            self.spec_type = "COUNTS"
        elif "RATE" in colnames:
            self.spec_type = "RATE"
        else:
            raise ValueError("Invalid spectrum file")
        self.channel          = ext_spec.data.columns["CHANNEL"].array
        col_spec_data         = ext_spec.data.columns[self.spec_type]
        self.spec_data        = col_spec_data.array.copy()
        self.spec_unit        = col_spec_data.unit
        self.spec_dtype       = col_spec_data.dtype
        self.spec_fits_format = col_spec_data.format
        # grouping and quality
        if "GROUPING" in colnames:
            self.grouping = ext_spec.data.columns["GROUPING"].array
        if "QUALITY"  in colnames:
            self.quality  = ext_spec.data.columns["QUALITY"].array
        # keywords
        self.EXPOSURE = self.header.get("EXPOSURE")
        self.BACKSCAL = self.header.get("BACKSCAL")
        self.AREASCAL = self.header.get("AREASCAL")
        self.RESPFILE = self.header.get("RESPFILE")
        self.ANCRFILE = self.header.get("ANCRFILE")
        self.BACKFILE = self.header.get("BACKFILE")
        # output filename
        self.outfile  = outfile

    def get_data(self, group_squeeze=False, copy=True):
        """
        Get the spectral data (i.e., self.spec_data).

        Arguments:
          * group_squeeze: whether squeeze the spectral data according to
                           the grouping (i.e., exclude the channels that
                           are not the first channel of the group, which
                           also have value of ZERO).
                           This argument is effective only the grouping
                           being applied.
        """
        if group_squeeze and self.groupped:
            spec_data = self.spec_data[self.grouping == 1]
        else:
            spec_data = self.spec_data
        if copy:
            return spec_data.copy()
        else:
            return spec_data

    def get_channel(self, copy=True):
        if copy:
            return self.channel.copy()
        else:
            return self.channel

    def set_data(self, spec_data, group_squeeze=True):
        """
        Set the spectral data of this spectrum to the supplied data.
        """
        if group_squeeze and self.groupped:
            assert sum(self.grouping == 1) == len(spec_data)
            self.spec_data[self.grouping == 1] = spec_data
        else:
            assert len(self.spec_data) == len(spec_data)
            self.spec_data = spec_data.copy()

    def add_stat_err(self, stat_err, group_squeeze=True):
        """
        Add the "STAT_ERR" column as the statistical errors of each spectral
        group, which are estimated by utilizing the Monte Carlo techniques.
        """
        self.stat_err = np.zeros(self.spec_data.shape,
                dtype=self.spec_data.dtype)
        if group_squeeze and self.groupped:
            assert sum(self.grouping == 1) == len(stat_err)
            self.stat_err[self.grouping == 1] = stat_err
        else:
            assert len(self.stat_err) == len(stat_err)
            self.stat_err = stat_err.copy()
        self.header["POISSERR"] = False

    def apply_grouping(self, grouping=None, quality=None):
        """
        Apply the spectral channel grouping specification to the spectrum.

        NOTE:
        * The spectral data (i.e., self.spec_data) is MODIFIED!
        * The spectral data within the same group are summed up.
        * The self grouping is overwritten if `grouping' is supplied, as well
          as the self quality.
        """
        if grouping is not None:
            self.grouping = grouping
        if quality is not None:
            self.quality = quality
        self.spec_data = group_data(self.spec_data, self.grouping)
        self.groupped = True

    def estimate_errors(self, gehrels=True):
        """
        Estimate the statistical errors of each spectral group (after
        applying grouping) for the source spectrum (and background spectrum).

        If `gehrels=True', the statistical error for a spectral group with
        N photons is given by `1 + sqrt(N + 0.75)'; otherwise, the error
        is given by `sqrt(N)'.

        Results: `self.spec_err'
        """
        eps = 1.0e-10
        if gehrels:
            self.spec_err = 1.0 + np.sqrt(self.spec_data + 0.75)
        else:
            self.spec_err = np.sqrt(self.spec_data)
        # replace the zeros with a very small value (because
        # `np.random.normal' requires `scale' > 0)
        self.spec_err[self.spec_err <= 0.0] = eps

    def copy(self):
        """
        Return a copy of this object, with the `np.ndarray' properties are
        copied.
        """
        new = copy(self)
        for k, v in self.__dict__.items():
            if isinstance(v, np.ndarray):
                setattr(new, k, v.copy())
        return new

    def randomize(self):
        """
        Randomize the spectral data according to the estimated spectral
        group errors by assuming the normal distribution.

        NOTE: this method should be called AFTER the `copy()' method.
        """
        if self.spec_err is None:
            raise ValueError("No valid 'spec_err' presents")
        if self.groupped:
            idx = self.grouping == 1
            self.spec_data[idx] = np.random.normal(self.spec_data[idx],
                    self.spec_err[idx])
        else:
            self.spec_data = np.random.normal(self.spec_data, self.spec_err)
        return self

    def fix_header_keywords(self,
            reset_kw=["ANCRFILE", "RESPFILE", "BACKFILE"]):
        """
        Reset the keywords to "NONE" to avoid confusion or mistakes,
        and also add mandatory spectral keywords if missing.

        Reference:
        [1] The OGIP Spectral File Format, Sec. 3.1.1
            https://heasarc.gsfc.nasa.gov/docs/heasarc/ofwg/docs/summary/ogip_92_007_summary.html
        """
        default_keywords = {
                ## Mandatory keywords
                #"EXTNAME"  : "SPECTRUM",
                "TELESCOP" : "NONE",
                "INSTRUME" : "NONE",
                "FILTER"   : "NONE",
                #"EXPOSURE" : <integration_time (s)>,
                "BACKFILE" : "NONE",
                "CORRFILE" : "NONE",
                "CORRSCAL" : 1.0,
                "RESPFILE" : "NONE",
                "ANCRFILE" : "NONE",
                "HDUCLASS" : "OGIP",
                "HDUCLAS1" : "SPECTRUM",
                "HDUVERS"  : "1.2.1",
                "POISSERR" : True,
                #"CHANTYPE" : "PI",
                #"DETCHANS" : <total_number_of_detector_channels>,
                ## Optional keywords for further information
                "BACKSCAL" : 1.0,
                "AREASCAL" : 1.0,
                # Type of spectral data:
                #   (1) "TOTAL": gross spectrum (source+bkg);
                #   (2) "NET": background-subtracted spectrum
                #   (3) "BKG" background spectrum
                #"HDUCLAS2" : "NET",
                # Details of the type of data:
                #   (1) "COUNT": data stored as counts
                #   (2) "RATE": data stored as counts/s
                "HDUCLAS3" : { "COUNTS":"COUNT",
                               "RATE":"RATE" }.get(self.spec_type),
        }
        # add mandatory keywords if missing
        for kw, value in default_keywords.items():
            if kw not in self.header:
                self.header[kw] = value
        # reset the specified keywords
        for kw in reset_kw:
            self.header[kw] = default_keywords.get(kw)

    def write(self, filename=None, clobber=False):
        """
        Create a new "SPECTRUM" table/extension and replace the original
        one, then write to output file.
        """
        if filename is None:
            filename = self.outfile
        columns = [
                fits.Column(name="CHANNEL", format="I", array=self.channel),
                fits.Column(name=self.spec_type, format=self.spec_fits_format,
                            unit=self.spec_unit, array=self.spec_data),
        ]
        if self.grouping is not None:
            columns.append(fits.Column(name="GROUPING",
                                       format="I", array=self.grouping))
        if self.quality  is not None:
            columns.append(fits.Column(name="QUALITY",
                                       format="I", array=self.quality))
        if self.stat_err is not None:
            columns.append(fits.Column(name="STAT_ERR", unit=self.spec_unit,
                                       format=self.spec_fits_format,
                                       array=self.stat_err))
        ext_spec_cols = fits.ColDefs(columns)
        ext_spec = fits.BinTableHDU.from_columns(ext_spec_cols,
                                                 header=self.header)
        self.fitsobj["SPECTRUM"] = ext_spec
        self.fitsobj.writeto(filename, clobber=clobber, checksum=True)
# class Spectrum }}}


class SpectrumSet(Spectrum):  # {{{
    """
    This class handles a set of spectrum, including the source spectrum,
    RMF, ARF, and the background spectrum.

    **NOTE**:
    The "COUNTS" column data are converted from "int32" to "float32",
    since this spectrum will be subtracted/compensated according to the
    ratios of ARFs.
    """
    # ARF object for this spectrum
    arf = None
    # RMF object for this spectrum
    rmf = None
    # background Spectrum object for this spectrum
    bkg = None
    # inner and outer radius of the region from which the spectrum extracted
    radius_inner = None
    radius_outer = None
    # total angular range of the spectral region
    angle = None

    # numpy dtype and FITS format code to which the spectrum data be
    # converted if the data is "COUNTS"
    #_spec_dtype       = np.float32
    #_spec_fits_format = "E"
    _spec_dtype       = np.float64
    _spec_fits_format = "D"

    def __init__(self, filename, outfile=None, arf=None, rmf=None, bkg=None):
        super().__init__(filename, outfile)
        # convert spectrum data type if necessary
        if self.spec_data.dtype != self._spec_dtype:
            self.spec_data        = self.spec_data.astype(self._spec_dtype)
            self.spec_dtype       = self._spec_dtype
            self.spec_fits_format = self._spec_fits_format
        if arf is not None:
            if isinstance(arf, ARF):
                self.arf = arf
            else:
                self.arf = ARF(arf)
        if rmf is not None:
            if isinstance(rmf, RMF):
                self.rmf = rmf
            else:
                self.rmf = RMF(rmf)
        if bkg is not None:
            if isinstance(bkg, Spectrum):
                self.bkg = bkg
            else:
                self.bkg = Spectrum(bkg)
            # convert background spectrum data type if necessary
            if self.bkg.spec_data.dtype != self._spec_dtype:
                self.bkg.spec_data = self.bkg.spec_data.astype(self._spec_dtype)
                self.bkg.spec_dtype = self._spec_dtype
                self.bkg.spec_fits_format = self._spec_fits_format

    def get_energy(self, mean="geometric"):
        """
        Get the energy values of each channel if RMF present.

        NOTE:
        The "E_MIN" and "E_MAX" columns of the RMF is required to calculate
        the spectrum channel energies.
        And the channel energies are generally different to the "ENERG_LO"
        and "ENERG_HI" of the corresponding ARF.
        """
        if self.rmf is None:
            return None
        else:
            return self.rmf.get_energy(mean=mean)

    def get_arf(self, mean="geometric", groupped=True, copy=True):
        """
        Get the interpolated ARF data w.r.t the spectral channel energies
        if the ARF presents.

        Arguments:
          * groupped: (bool) whether to get the groupped ARF

        Return: (groupped) interpolated ARF data
        """
        if self.arf is None:
            return None
        else:
            return self.arf.get_data(groupped=groupped, copy=copy)

    def read_xflt(self):
        """
        Read the XFLT000# keywords from the header, check the validity (e.g.,
        "XFLT0001" should equals "XFLT0002", "XFLT0003" should equals 0).
        Sum all the additional XFLT000# pairs (e.g., ) which describes the
        regions angluar ranges.
        """
        eps = 1.0e-6
        xflt0001 = float(self.header["XFLT0001"])
        xflt0002 = float(self.header["XFLT0002"])
        xflt0003 = float(self.header["XFLT0003"])
        # XFLT000# validity check
        assert np.isclose(xflt0001, xflt0002)
        assert abs(xflt0003) < eps
        # outer radius of the region
        self.radius_outer = xflt0001
        # angular regions
        self.angle = 0.0
        num = 4
        while True:
            try:
                angle_begin = float(self.header["XFLT%04d" % num])
                angle_end   = float(self.header["XFLT%04d" % (num+1)])
                num += 2
            except KeyError:
                break
            self.angle += (angle_end - angle_begin)
        # if NO additional XFLT000# keys exist, assume "annulus" region
        if self.angle < eps:
            self.angle = 360.0

    def scale(self):
        """
        Scale the spectral data (and spectral group errors if present) of
        the source spectrum (and background spectra if present) according
        to the region angular size to make it correspond to the whole annulus
        region (i.e., 360 degrees).

        NOTE: the spectral data and errors (i.e., `self.spec_data', and
        `self.spec_err') is MODIFIED!
        """
        self.spec_data *= (360.0 / self.angle)
        if self.spec_err is not None:
            self.spec_err *= (360.0 / self.angle)
        # also scale the background spectrum if present
        if self.bkg:
            self.bkg.spec_data *= (360.0 / self.angle)
            if self.bkg.spec_err is not None:
                self.bkg.spec_err *= (360.0 / self.angle)

    def apply_grouping(self, grouping=None, quality=None, verbose=False):
        """
        Apply the spectral channel grouping specification to the source
        spectrum, the ARF (which is used during the later spectral
        manipulations), and the background spectrum (if presents).

        NOTE:
        * The spectral data (i.e., self.spec_data) is MODIFIED!
        * The spectral data within the same group are summed up.
        * The self grouping is overwritten if `grouping' is supplied, as well
          as the self quality.
        """
        super().apply_grouping(grouping=grouping, quality=quality)
        # also group the ARF accordingly
        self.arf.apply_grouping(energy_channel=self.get_energy(),
                                grouping=self.grouping, verbose=verbose)
        # group the background spectrum if present
        if self.bkg:
            self.bkg.spec_data = group_data(self.bkg.spec_data, self.grouping)

    def estimate_errors(self, gehrels=True):
        """
        Estimate the statistical errors of each spectral group (after
        applying grouping) for the source spectrum (and background spectrum).

        If `gehrels=True', the statistical error for a spectral group with
        N photons is given by `1 + sqrt(N + 0.75)'; otherwise, the error
        is given by `sqrt(N)'.

        Results: `self.spec_err' (and `self.bkg.spec_err')
        """
        super().estimate_errors(gehrels=gehrels)
        eps = 1.0e-10
        # estimate the errors for background spectrum if present
        if self.bkg:
            if gehrels:
                self.bkg.spec_err = 1.0 + np.sqrt(self.bkg.spec_data + 0.75)
            else:
                self.bkg.spec_err = np.sqrt(self.bkg.spec_data)
            self.bkg.spec_err[self.bkg.spec_err <= 0.0] = eps

    def subtract_bkg(self, inplace=True, add_history=False, verbose=False):
        """
        Subtract the background contribution from the source spectrum.
        The `EXPOSURE' and `BACKSCAL' values are required to calculate
        the fraction/ratio for the background subtraction.

        Arguments:
          * inplace: whether replace the `spec_data' with the background-
                     subtracted spectrum data; If True, the attribute
                     `spec_bkg_subtracted' is also set to `True' when
                     the subtraction finished.
                     The keywords "BACKSCAL" and "AREASCAL" are set to 1.0.

        Return:
          background-subtracted spectrum data
        """
        ratio = (self.EXPOSURE / self.bkg.EXPOSURE) * \
                (self.BACKSCAL / self.bkg.BACKSCAL) * \
                (self.AREASCAL / self.bkg.AREASCAL)
        operation = "  SUBTRACT_BACKGROUND: %s - %s * %s" % \
                (self.filename, ratio, self.bkg.filename)
        if verbose:
            print(operation, file=sys.stderr)
        spec_data_subbkg = self.spec_data - ratio * self.bkg.get_data()
        if inplace:
            self.spec_data = spec_data_subbkg
            self.spec_bkg_subtracted = True
            self.BACKSCAL = 1.0
            self.AREASCAL = 1.0
            # update header
            self.header["BACKSCAL"] = 1.0
            self.header["AREASCAL"] = 1.0
            self.header["BACKFILE"] = "NONE"
            self.header["HDUCLAS2"] = "NET"  # background-subtracted spectrum
            # also record history
            if add_history:
                self.header.add_history(operation)
        return spec_data_subbkg

    def subtract(self, spectrumset, cross_arf, groupped=False,
            group_squeeze=False, add_history=False, verbose=False):
        """
        Subtract the photons that originate from the surrounding regions
        but were scattered into this spectrum due to the finite PSF.

        The background of this spectrum and the given spectrum should
        both be subtracted before applying this subtraction for crosstalk
        correction, as well as the below `compensate()' procedure.

        NOTE:
        1. The crosstalk ARF must be provided, since the `spectrumset.arf'
           is required to be its ARF without taking crosstalk into account:
               spec1_new = spec1 - spec2 * (cross_arf_2_to_1 / arf2)
        2. The ARF are interpolated to match the energies of spetral channels.
        """
        operation = "  SUBTRACT: %s - (%s/%s) * %s" % (self.filename,
                cross_arf.filename, spectrumset.arf.filename,
                spectrumset.filename)
        if verbose:
            print(operation, file=sys.stderr)
        energy = self.get_energy()
        if groupped:
            spectrumset.arf.apply_grouping(energy_channel=energy,
                    grouping=self.grouping, verbose=verbose)
            cross_arf.apply_grouping(energy_channel=energy,
                    grouping=self.grouping, verbose=verbose)
            arfresp_spec  = spectrumset.arf.get_data(groupped=True,
                    group_squeeze=group_squeeze)
            arfresp_cross = cross_arf.get_data(groupped=True,
                    group_squeeze=group_squeeze)
        else:
            arfresp_spec  = spectrumset.arf.interpolate(x=energy,
                    verbose=verbose)
            arfresp_cross = cross_arf.interpolate(x=energy, verbose=verbose)
        with np.errstate(divide="ignore", invalid="ignore"):
            arf_ratio = arfresp_cross / arfresp_spec
            # fix nan/inf values due to division by zero
            arf_ratio[ ~ np.isfinite(arf_ratio) ] = 0.0
        spec_data = self.get_data(group_squeeze=group_squeeze) - \
                spectrumset.get_data(group_squeeze=group_squeeze)*arf_ratio
        self.set_data(spec_data, group_squeeze=group_squeeze)
        # record history
        if add_history:
            self.header.add_history(operation)

    def compensate(self, cross_arf, groupped=False, group_squeeze=False,
            add_history=False, verbose=False):
        """
        Compensate the photons that originate from this regions but were
        scattered into the surrounding regions due to the finite PSF.

        formula:
            spec1_new = spec1 + spec1 * (cross_arf_1_to_2 / arf1)
        """
        operation = "  COMPENSATE: %s + (%s/%s) * %s" % (self.filename,
                cross_arf.filename, self.arf.filename, self.filename)
        if verbose:
            print(operation, file=sys.stderr)
        energy = self.get_energy()
        if groupped:
            cross_arf.apply_grouping(energy_channel=energy,
                    grouping=self.grouping, verbose=verbose)
            arfresp_this  = self.arf.get_data(groupped=True,
                    group_squeeze=group_squeeze)
            arfresp_cross = cross_arf.get_data(groupped=True,
                    group_squeeze=group_squeeze)
        else:
            arfresp_this  = self.arf.interpolate(x=energy, verbose=verbose)
            arfresp_cross = cross_arf.interpolate(x=energy, verbose=verbose)
        with np.errstate(divide="ignore", invalid="ignore"):
            arf_ratio = arfresp_cross / arfresp_this
            # fix nan/inf values due to division by zero
            arf_ratio[ ~ np.isfinite(arf_ratio) ] = 0.0
        spec_data = self.get_data(group_squeeze=group_squeeze) + \
                self.get_data(group_squeeze=group_squeeze) * arf_ratio
        self.set_data(spec_data, group_squeeze=group_squeeze)
        # record history
        if add_history:
            self.header.add_history(operation)

    def fix_negative(self, add_history=False, verbose=False):
        """
        The subtractions may lead to negative counts, it may be necessary
        to fix these channels with negative values.
        """
        neg_counts = self.spec_data < 0
        N = len(neg_counts)
        neg_channels = np.arange(N, dtype=np.int)[neg_counts]
        if len(neg_channels) > 0:
            print("WARNING: %d channels have NEGATIVE counts" % \
                    len(neg_channels), file=sys.stderr)
        i = 0
        while len(neg_channels) > 0:
            i += 1
            if verbose:
                if i == 1:
                    print("*** Fixing negative channels: iter %d..." % i,
                            end="", file=sys.stderr)
                else:
                    print("%d..." % i, end="", file=sys.stderr)
            for ch in neg_channels:
                neg_val = self.spec_data[ch]
                if ch < N-2:
                    self.spec_data[ch] = 0
                    self.spec_data[(ch+1):(ch+3)] -= 0.5 * np.abs(neg_val)
                else:
                    # just set to zero if it is the last 2 channels
                    self.spec_data[ch] = 0
            # update negative channels indices
            neg_counts = self.spec_data < 0
            neg_channels = np.arange(N, dtype=np.int)[neg_counts]
        if i > 0:
            print("FIXED!", file=sys.stderr)
            # record history
            if add_history:
                self.header.add_history("  FIXED NEGATIVE CHANNELS")

    def set_radius_inner(self, radius_inner):
        """
        Set the inner radius of the spectral region.
        """
        assert radius_inner < self.radius_outer
        self.radius_inner = radius_inner

    def copy(self):
        """
        Return a copy of this object.
        """
        new = super().copy()
        if self.bkg:
            new.bkg = self.bkg.copy()
        return new

    def randomize(self):
        """
        Randomize the source (and background if present) spectral data
        according to the estimated spectral group errors by assuming the
        normal distribution.

        NOTE: this method should be called AFTER the `copy()' method.
        """
        super().randomize()
        if self.bkg:
            self.bkg.spec_data = np.random.normal(self.bkg.spec_data,
                                                  self.bkg.spec_err)
            self.bkg.spec_data[self.grouping == -1] = 0.0
        return self
# class SpectrumSet }}}


class Crosstalk:  # {{{
    """
    XMM-Newton PSF Crosstalk effect correction.
    """
    # `SpectrumSet' object for the spectrum to be corrected
    spectrumset = None
    # NOTE/XXX: do NOT use list (e.g., []) here, otherwise, all the
    #           instances will share these list properties.
    # `SpectrumSet' and `ARF' objects corresponding to the spectra from
    # which the photons were scattered into this spectrum.
    cross_in_specset = None
    cross_in_arf     = None
    # `ARF' objects corresponding to the regions to which the photons of
    # this spectrum were scattered into.
    cross_out_arf = None
    # grouping specification and quality data
    grouping = None
    quality = None
    # whether the spectrum is groupped
    groupped = False

    def __init__(self, config, arf_dict={}, rmf_dict={},
            grouping=None, quality=None):
        """
        Arguments:
          * config: a section of the whole config file (`ConfigObj' object)
        """
        self.cross_in_specset = []
        self.cross_in_arf     = []
        self.cross_out_arf    = []
        # this spectrum to be corrected
        self.spectrumset = SpectrumSet(filename=config["spec"],
                outfile=config["outfile"],
                arf=arf_dict.get(config["arf"], config["arf"]),
                rmf=rmf_dict.get(config.get("rmf"), config.get("rmf")),
                bkg=config.get("bkg"))
        # spectra and cross arf from which photons were scattered in
        for reg_in in config["cross_in"].values():
            specset = SpectrumSet(filename=reg_in["spec"],
                    arf=arf_dict.get(reg_in["arf"], reg_in["arf"]),
                    rmf=rmf_dict.get(reg_in.get("rmf"), reg_in.get("rmf")),
                    bkg=reg_in.get("bkg"))
            self.cross_in_specset.append(specset)
            self.cross_in_arf.append(arf_dict.get(reg_in["cross_arf"],
                                                  ARF(reg_in["cross_arf"])))
        # regions into which the photons of this spectrum were scattered into
        if "cross_out" in config.sections:
            cross_arf = config["cross_out"].as_list("cross_arf")
            for arffile in cross_arf:
                self.cross_out_arf.append(arf_dict.get(arffile, ARF(arffile)))
        # grouping and quality
        self.grouping = grouping
        self.quality = quality

    def apply_grouping(self, verbose=False):
        self.spectrumset.apply_grouping(grouping=self.grouping,
                quality=self.quality, verbose=verbose)
        # also group the related surrounding spectra
        for specset in self.cross_in_specset:
            specset.apply_grouping(grouping=self.grouping,
                quality=self.quality, verbose=verbose)
        self.groupped = True

    def estimate_errors(self, gehrels=True, verbose=False):
        if verbose:
            print("INFO: Estimating spectral errors ...")
        self.spectrumset.estimate_errors(gehrels=gehrels)
        # also estimate errors for the related surrounding spectra
        for specset in self.cross_in_specset:
            specset.estimate_errors(gehrels=gehrels)

    def do_correction(self, subtract_bkg=True, fix_negative=False,
            group_squeeze=True, add_history=False, verbose=False):
        """
        Perform the crosstalk correction.  The background contribution
        for each spectrum is subtracted first if `subtract_bkg' is True.
        The basic correction procedures are recorded to the header.
        """
        if add_history:
            self.spectrumset.header.add_history("Crosstalk Correction BEGIN")
            self.spectrumset.header.add_history("  TOOL: %s (v%s) @ %s" % (\
                    os.path.basename(sys.argv[0]), __version__,
                    datetime.utcnow().isoformat()))
        # background subtraction
        if subtract_bkg:
            if verbose:
                print("INFO: subtract background ...", file=sys.stderr)
            self.spectrumset.subtract_bkg(inplace=True,
                    add_history=add_history, verbose=verbose)
            # also apply background subtraction to the surrounding spectra
            for specset in self.cross_in_specset:
                specset.subtract_bkg(inplace=True,
                        add_history=add_history, verbose=verbose)
        # subtractions
        if verbose:
            print("INFO: apply subtractions ...", file=sys.stderr)
        for specset, cross_arf in zip(self.cross_in_specset,
                self.cross_in_arf):
            self.spectrumset.subtract(spectrumset=specset,
                    cross_arf=cross_arf, groupped=self.groupped,
                    group_squeeze=group_squeeze, add_history=add_history,
                    verbose=verbose)
        # compensations
        if verbose:
            print("INFO: apply compensations ...", file=sys.stderr)
        for cross_arf in self.cross_out_arf:
            self.spectrumset.compensate(cross_arf=cross_arf,
                    groupped=self.groupped, group_squeeze=group_squeeze,
                    add_history=add_history, verbose=verbose)
        # fix negative values in channels
        if fix_negative:
            if verbose:
                print("INFO: fix negative channel values ...", file=sys.stderr)
            self.spectrumset.fix_negative(add_history=add_history,
                    verbose=verbose)
        if add_history:
            self.spectrumset.header.add_history("END Crosstalk Correction")

    def fix_header(self):
        # fix header keywords
        self.spectrumset.fix_header_keywords(
                reset_kw=["RESPFILE", "ANCRFILE", "BACKFILE"])

    def copy(self):
        new = copy(self)
        # properly handle the copy of spectrumsets
        new.spectrumset = self.spectrumset.copy()
        new.cross_in_specset = [ specset.copy() \
                for specset in self.cross_in_specset ]
        return new

    def randomize(self):
        self.spectrumset.randomize()
        for specset in self.cross_in_specset:
            specset.randomize()
        return self

    def get_spectrum(self, copy=True):
        if copy:
            return self.spectrumset.copy()
        else:
            return self.spectrumset

    def write(self, filename=None, clobber=False):
        self.spectrumset.write(filename=filename, clobber=clobber)
# class Crosstalk }}}


class Deprojection:  # {{{
    """
    Perform the deprojection on a set of PROJECTED spectra with the
    assumption of spherical symmetry of the source object, and produce
    the DEPROJECTED spectra.

    NOTE:
    * Assumption of the spherical symmetry
    * Background should be subtracted before deprojection
    * ARF differences of different regions are taken into account

    Reference & Credit:
    [1] Direct X-ray Spectra Deprojection
        https://www-xray.ast.cam.ac.uk/papers/dsdeproj/
        Sanders & Fabian 2007, MNRAS, 381, 1381
    """
    spectra  = None
    grouping = None
    quality  = None

    def __init__(self, spectra, grouping=None, quality=None, verbose=False):
        """
        Arguments:
          * spectra: a set of spectra from the inner-most to the outer-most
                     regions (e.g., spectra after correcting crosstalk effect)
          * grouping: grouping specification for all the spectra
          * quality: quality column for the spectra
        """
        self.spectra = []
        for spec in spectra:
            if not isinstance(spec, SpectrumSet):
                raise ValueError("Not a 'SpectrumSet' object")
            spec.read_xflt()
            self.spectra.append(spec)
        self.spectra  = spectra
        self.grouping = grouping
        self.quality  = quality
        # sort spectra by `radius_outer'
        self.spectra.sort(key=lambda x: x.radius_outer)
        # set the inner radii
        radii_inner = [0.0] + [ x.radius_outer for x in self.spectra[:-1] ]
        for spec, rin in zip(self.spectra, radii_inner):
            spec.set_radius_inner(rin)
            if verbose:
                print("Deprojection: loaded spectrum: radius: (%s, %s)" % \
                        (spec.radius_inner, spec.radius_outer),
                        file=sys.stderr)
        # check EXPOSURE validity (all spectra must have the same exposures)
        exposures = [ spec.EXPOSURE for spec in self.spectra ]
        assert np.allclose(exposures[:-1], exposures[1:])

    def subtract_bkg(self, verbose=True):
        for spec in self.spectra:
            if not spec.bkg:
                raise ValueError("Spectrum '%s' has NO background" % \
                        spec.filename)
            spec.subtract_bkg(inplace=True, verbose=verbose)

    def apply_grouping(self, verbose=False):
        for spec in self.spectra:
            spec.apply_grouping(grouping=self.grouping, quality=self.quality,
                    verbose=verbose)

    def estimate_errors(self, gehrels=True):
        for spec in self.spectra:
            spec.estimate_errors(gehrels=gehrels)

    def scale(self):
        """
        Scale the spectral data according to the region angular size.
        """
        for spec in self.spectra:
            spec.scale()

    def do_deprojection(self, group_squeeze=True,
            add_history=False, verbose=False):
        #
        # TODO/XXX: How to apply ARF correction here???
        #
        num_spec = len(self.spectra)
        tmp_spec_data = self.spectra[0].get_data(group_squeeze=group_squeeze)
        spec_shape = tmp_spec_data.shape
        spec_dtype = tmp_spec_data.dtype
        spec_per_vol = [None] * num_spec
        #
        for shellnum in reversed(range(num_spec)):
            if verbose:
                print("DEPROJECTION: deprojecting shell %d ..." % shellnum,
                        file=sys.stderr)
            spec = self.spectra[shellnum]
            # calculate projected spectrum of outlying shells
            proj_spec = np.zeros(spec_shape, spec_dtype)
            for outer in range(shellnum+1, num_spec):
                vol = self.projected_volume(
                        r1=self.spectra[outer].radius_inner,
                        r2=self.spectra[outer].radius_outer,
                        R1=spec.radius_inner,
                        R2=spec.radius_outer)
                proj_spec += spec_per_vol[outer] * vol
            #
            this_spec = spec.get_data(group_squeeze=group_squeeze, copy=True)
            deproj_spec = this_spec - proj_spec
            # calculate the volume that this spectrum is from
            this_vol = self.projected_volume(
                    r1=spec.radius_inner, r2=spec.radius_outer,
                    R1=spec.radius_inner, R2=spec.radius_outer)
            # calculate the spectral data per unit volume
            spec_per_vol[shellnum] = deproj_spec / this_vol
        # set the spectral data to these deprojected values
        self.set_spec_data(spec_per_vol, group_squeeze=group_squeeze)
        # add history to header
        if add_history:
            self.add_history()

    def get_spec_data(self, group_squeeze=True, copy=True):
        """
        Extract the spectral data of each spectrum after deprojection
        performed.
        """
        return [ spec.get_data(group_squeeze=group_squeeze, copy=copy)
                 for spec in self.spectra ]

    def set_spec_data(self, spec_data, group_squeeze=True):
        """
        Set `spec_data' for each spectrum to the deprojected spectral data.
        """
        assert len(spec_data) == len(self.spectra)
        for spec, data in zip(self.spectra, spec_data):
            spec.set_data(data, group_squeeze=group_squeeze)

    def add_stat_err(self, stat_err, group_squeeze=True):
        """
        Add the "STAT_ERR" column to each spectrum.
        """
        assert len(stat_err) == len(self.spectra)
        for spec, err in zip(self.spectra, stat_err):
            spec.add_stat_err(err, group_squeeze=group_squeeze)

    def add_history(self):
        """
        Append a brief history about this tool to the header.
        """
        history = "Deprojected by %s (v%s) @ %s" % (
                os.path.basename(sys.argv[0]), __version__,
                datetime.utcnow().isoformat())
        for spec in self.spectra:
            spec.header.add_history(history)

    def fix_header(self):
        # fix header keywords
        for spec in self.spectra:
            spec.fix_header_keywords(
                    reset_kw=["RESPFILE", "ANCRFILE", "BACKFILE"])

    def write(self, filenames=[], clobber=False):
        """
        Write the deprojected spectra to output file.
        """
        if filenames == []:
            filenames = [ spec.outfile for spec in self.spectra ]
        for spec, outfile in zip(self.spectra, filenames):
            spec.write(filename=outfile, clobber=clobber)

    @staticmethod
    def projected_volume(r1, r2, R1, R2):
        """
        Calculate the projected volume of a spherical shell of radii r1 -> r2
        onto an annulus on the sky of radius R1 -> R2.

        This volume is the integral:
          Int(R=R1,R2) Int(x=sqrt(r1^2-R^2),sqrt(r2^2-R^2)) 2*pi*R dx dR
          =
          Int(R=R1,R2) 2*pi*R * (sqrt(r2^2-R^2) - sqrt(r1^2-R^2)) dR

        Note that the above integral is only half the total volume
        (i.e., front only).
        """
        def sqrt_trunc(x):
            if x > 0:
                return np.sqrt(x)
            else:
                return 0.0
        #
        p1 = sqrt_trunc(r1**2 - R2**2)
        p2 = sqrt_trunc(r1**2 - R1**2)
        p3 = sqrt_trunc(r2**2 - R2**2)
        p4 = sqrt_trunc(r2**2 - R1**2)
        return 2.0 * (2.0/3.0) * np.pi * ((p1**3 - p2**3) + (p4**3 - p3**3))
# class Deprojection }}}


# Helper functions {{{
def calc_median_errors(results):
    """
    Calculate the median and errors for the spectral data gathered
    through Monte Carlo simulations.

    NOTE:
    Errors calculation just use the quantiles,
    i.e., 1sigma ~= 68.3% = Q(84.15%) - Q(15.85%)
    """
    results = np.array(results)
    # `results' now has shape: (mc_times, num_spec, num_channel)
    # sort by the Monte Carlo simulation axis
    results.sort(0)
    mc_times = results.shape[0]
    medians  = results[ int(mc_times * 0.5) ]
    lowerpcs = results[ int(mc_times * 0.1585) ]
    upperpcs = results[ int(mc_times * 0.8415) ]
    errors   = np.sqrt(0.5 * ((medians-lowerpcs)**2 + (upperpcs-medians)**2))
    return (medians, errors)


def set_argument(name, default, cmdargs, config):
    value = default
    if name in config.keys():
        value = config.as_bool(name)
    value_cmd = vars(cmdargs)[name]
    if value_cmd != default:
        value = value_cmd  # command arguments overwrite others
    return value
# helper functions }}}


# main routine {{{
def main(config, subtract_bkg, fix_negative, mc_times,
        verbose=False, clobber=False):
    # collect ARFs and RMFs into dictionaries (avoid interpolation every time)
    arf_files = set()
    rmf_files = set()
    for region in config.sections:
        config_reg = config[region]
        arf_files.add(config_reg.get("arf"))
        rmf_files.add(config_reg.get("rmf"))
        for reg_in in config_reg["cross_in"].values():
            arf_files.add(reg_in.get("arf"))
            arf_files.add(reg_in.get("cross_arf"))
        if "cross_out" in config_reg.sections:
            for arf in config_reg["cross_out"].as_list("cross_arf"):
                arf_files.add(arf)
    arf_files = arf_files - set([None])
    arf_dict  = { arf: ARF(arf) for arf in arf_files }
    rmf_files = rmf_files - set([None])
    rmf_dict  = { rmf: RMF(rmf) for rmf in rmf_files }
    if verbose:
        print("INFO: arf_files:", arf_files, file=sys.stderr)
        print("INFO: rmf_files:", rmf_files, file=sys.stderr)

    # get the GROUPING and QUALITY data
    grouping_fits = fits.open(config["grouping"])
    grouping = grouping_fits["SPECTRUM"].data.columns["GROUPING"].array
    quality  = grouping_fits["SPECTRUM"].data.columns["QUALITY"].array
    # squeeze the groupped spectral data, etc.
    group_squeeze = True

    # crosstalk objects (BEFORE background subtraction)
    crosstalks_cleancopy = []
    # crosstalk-corrected spectra
    cc_spectra = []

    # correct crosstalk effects for each region first
    for region in config.sections:
        if verbose:
            print("INFO: processing '%s' ..." % region, file=sys.stderr)
        crosstalk = Crosstalk(config.get(region),
                arf_dict=arf_dict, rmf_dict=rmf_dict,
                grouping=grouping, quality=quality)
        crosstalk.apply_grouping(verbose=verbose)
        crosstalk.estimate_errors(verbose=verbose)
        # keep a (almost) clean copy of the crosstalk object
        crosstalks_cleancopy.append(crosstalk.copy())
        if verbose:
            print("INFO: doing crosstalk correction ...", file=sys.stderr)
        crosstalk.do_correction(subtract_bkg=subtract_bkg,
                fix_negative=fix_negative, group_squeeze=group_squeeze,
                add_history=True, verbose=verbose)
        cc_spectra.append(crosstalk.get_spectrum(copy=True))

    # load back the crosstalk-corrected spectra for deprojection
    if verbose:
        print("INFO: preparing spectra for deprojection ...", file=sys.stderr)
    deprojection = Deprojection(spectra=cc_spectra, grouping=grouping,
            quality=quality, verbose=verbose)
    if verbose:
        print("INFO: scaling spectra according the region angular size...",
                file=sys.stderr)
    deprojection.scale()
    if verbose:
        print("INFO: doing deprojection ...", file=sys.stderr)
    deprojection.do_deprojection(add_history=True, verbose=verbose)
    deproj_results = [ deprojection.get_spec_data(
            group_squeeze=group_squeeze, copy=True) ]

    # Monte Carlo for spectral group error estimation
    print("INFO: Monte Carlo to estimate spectral errors (%d times) ..." % \
            mc_times, file=sys.stderr)
    for i in range(mc_times):
        if i % 100 == 0:
            print("%d..." % i, end="", flush=True, file=sys.stderr)
        # correct crosstalk effects
        cc_spectra_copy = []
        for crosstalk in crosstalks_cleancopy:
            # copy and randomize
            crosstalk_copy = crosstalk.copy().randomize()
            crosstalk_copy.do_correction(subtract_bkg=subtract_bkg,
                    fix_negative=fix_negative, group_squeeze=group_squeeze,
                    add_history=False, verbose=False)
            cc_spectra_copy.append(crosstalk_copy.get_spectrum(copy=True))
        # deproject spectra
        deprojection_copy = Deprojection(spectra=cc_spectra_copy,
                grouping=grouping, quality=quality, verbose=False)
        deprojection_copy.scale()
        deprojection_copy.do_deprojection(add_history=False, verbose=False)
        deproj_results.append(deprojection_copy.get_spec_data(
                group_squeeze=group_squeeze, copy=True))
    print("DONE!", flush=True, file=sys.stderr)

    if verbose:
        print("INFO: Calculating the median and errors for each spectrum ...",
                file=sys.stderr)
    medians, errors = calc_median_errors(deproj_results)
    deprojection.set_spec_data(medians, group_squeeze=group_squeeze)
    deprojection.add_stat_err(errors, group_squeeze=group_squeeze)
    if verbose:
        print("INFO: Writing the crosstalk-corrected and deprojected " + \
              "spectra with estimated statistical errors ...", file=sys.stderr)
    deprojection.fix_header()
    deprojection.write(clobber=clobber)
# main routine }}}


# main_deprojection routine {{{
def main_deprojection(config, mc_times, verbose=False, clobber=False):
    """
    Only perform the spectral deprojection.
    """
    # collect ARFs and RMFs into dictionaries (avoid interpolation every time)
    arf_files = set()
    rmf_files = set()
    for region in config.sections:
        config_reg = config[region]
        arf_files.add(config_reg.get("arf"))
        rmf_files.add(config_reg.get("rmf"))
    arf_files = arf_files - set([None])
    arf_dict  = { arf: ARF(arf) for arf in arf_files }
    rmf_files = rmf_files - set([None])
    rmf_dict  = { rmf: RMF(rmf) for rmf in rmf_files }
    if verbose:
        print("INFO: arf_files:", arf_files, file=sys.stderr)
        print("INFO: rmf_files:", rmf_files, file=sys.stderr)

    # get the GROUPING and QUALITY data
    grouping_fits = fits.open(config["grouping"])
    grouping = grouping_fits["SPECTRUM"].data.columns["GROUPING"].array
    quality  = grouping_fits["SPECTRUM"].data.columns["QUALITY"].array
    # squeeze the groupped spectral data, etc.
    group_squeeze = True

    # load spectra for deprojection
    if verbose:
        print("INFO: preparing spectra for deprojection ...", file=sys.stderr)
    proj_spectra = []
    for region in config.sections:
        config_reg = config[region]
        specset = SpectrumSet(filename=config_reg["spec"],
                outfile=config_reg["outfile"],
                arf=arf_dict.get(config_reg["arf"], config_reg["arf"]),
                rmf=rmf_dict.get(config_reg["rmf"], config_reg["rmf"]),
                bkg=config_reg["bkg"])
        proj_spectra.append(specset)

    deprojection = Deprojection(spectra=proj_spectra, grouping=grouping,
            quality=quality, verbose=verbose)
    deprojection.apply_grouping(verbose=verbose)
    deprojection.estimate_errors()
    if verbose:
        print("INFO: scaling spectra according the region angular size ...",
                file=sys.stderr)
    deprojection.scale()

    # keep a (almost) clean copy of the input projected spectra
    proj_spectra_cleancopy = [ spec.copy() for spec in proj_spectra ]

    if verbose:
        print("INFO: subtract the background ...", file=sys.stderr)
    deprojection.subtract_bkg(verbose=verbose)
    if verbose:
        print("INFO: doing deprojection ...", file=sys.stderr)
    deprojection.do_deprojection(add_history=True, verbose=verbose)
    deproj_results = [ deprojection.get_spec_data(
            group_squeeze=group_squeeze, copy=True) ]

    # Monte Carlo for spectral group error estimation
    print("INFO: Monte Carlo to estimate spectral errors (%d times) ..." % \
            mc_times, file=sys.stderr)
    for i in range(mc_times):
        if i % 100 == 0:
            print("%d..." % i, end="", flush=True, file=sys.stderr)
        # copy and randomize the input projected spectra
        proj_spectra_copy = [ spec.copy().randomize()
                              for spec in proj_spectra_cleancopy ]
        # deproject spectra
        deprojection_copy = Deprojection(spectra=proj_spectra_copy,
                grouping=grouping, quality=quality, verbose=False)
        deprojection_copy.subtract_bkg(verbose=False)
        deprojection_copy.do_deprojection(add_history=False, verbose=False)
        deproj_results.append(deprojection_copy.get_spec_data(
                group_squeeze=group_squeeze, copy=True))
    print("DONE!", flush=True, file=sys.stderr)

    if verbose:
        print("INFO: Calculating the median and errors for each spectrum ...",
                file=sys.stderr)
    medians, errors = calc_median_errors(deproj_results)
    deprojection.set_spec_data(medians, group_squeeze=group_squeeze)
    deprojection.add_stat_err(errors, group_squeeze=group_squeeze)
    if verbose:
        print("INFO: Writing the deprojected spectra " + \
              "with estimated statistical errors ...", file=sys.stderr)
    deprojection.fix_header()
    deprojection.write(clobber=clobber)
# main_deprojection routine }}}


# main_crosstalk routine {{{
def main_crosstalk(config, subtract_bkg, fix_negative, mc_times,
        verbose=False, clobber=False):
    """
    Only perform the crosstalk correction.
    """
    # collect ARFs and RMFs into dictionaries (avoid interpolation every time)
    arf_files = set()
    rmf_files = set()
    for region in config.sections:
        config_reg = config[region]
        arf_files.add(config_reg.get("arf"))
        rmf_files.add(config_reg.get("rmf"))
        for reg_in in config_reg["cross_in"].values():
            arf_files.add(reg_in.get("arf"))
            arf_files.add(reg_in.get("cross_arf"))
        if "cross_out" in config_reg.sections:
            for arf in config_reg["cross_out"].as_list("cross_arf"):
                arf_files.add(arf)
    arf_files = arf_files - set([None])
    arf_dict  = { arf: ARF(arf) for arf in arf_files }
    rmf_files = rmf_files - set([None])
    rmf_dict  = { rmf: RMF(rmf) for rmf in rmf_files }
    if verbose:
        print("INFO: arf_files:", arf_files, file=sys.stderr)
        print("INFO: rmf_files:", rmf_files, file=sys.stderr)

    # get the GROUPING and QUALITY data
    if "grouping" in config.keys():
        grouping_fits = fits.open(config["grouping"])
        grouping = grouping_fits["SPECTRUM"].data.columns["GROUPING"].array
        quality  = grouping_fits["SPECTRUM"].data.columns["QUALITY"].array
        group_squeeze = True
    else:
        grouping = None
        quality  = None
        group_squeeze = False

    # crosstalk objects (BEFORE background subtraction)
    crosstalks_cleancopy = []
    # crosstalk-corrected spectra
    cc_spectra = []

    # correct crosstalk effects for each region first
    for region in config.sections:
        if verbose:
            print("INFO: processing '%s' ..." % region, file=sys.stderr)
        crosstalk = Crosstalk(config.get(region),
                arf_dict=arf_dict, rmf_dict=rmf_dict,
                grouping=grouping, quality=quality)
        if grouping is not None:
            crosstalk.apply_grouping(verbose=verbose)
        crosstalk.estimate_errors(verbose=verbose)
        # keep a (almost) clean copy of the crosstalk object
        crosstalks_cleancopy.append(crosstalk.copy())
        if verbose:
            print("INFO: doing crosstalk correction ...", file=sys.stderr)
        crosstalk.do_correction(subtract_bkg=subtract_bkg,
                fix_negative=fix_negative, group_squeeze=group_squeeze,
                add_history=True, verbose=verbose)
        crosstalk.fix_header()
        cc_spectra.append(crosstalk.get_spectrum(copy=True))

    # spectral data of the crosstalk-corrected spectra
    cc_results = []
    cc_results.append([ spec.get_data(group_squeeze=group_squeeze, copy=True)
                        for spec in cc_spectra ])

    # Monte Carlo for spectral group error estimation
    print("INFO: Monte Carlo to estimate spectral errors (%d times) ..." % \
            mc_times, file=sys.stderr)
    for i in range(mc_times):
        if i % 100 == 0:
            print("%d..." % i, end="", flush=True, file=sys.stderr)
        # correct crosstalk effects
        cc_spectra_copy = []
        for crosstalk in crosstalks_cleancopy:
            # copy and randomize
            crosstalk_copy = crosstalk.copy().randomize()
            crosstalk_copy.do_correction(subtract_bkg=subtract_bkg,
                    fix_negative=fix_negative, group_squeeze=group_squeeze,
                    add_history=False, verbose=False)
            cc_spectra_copy.append(crosstalk_copy.get_spectrum(copy=True))
        cc_results.append([ spec.get_data(group_squeeze=group_squeeze,
                                          copy=True)
                            for spec in cc_spectra_copy ])
    print("DONE!", flush=True, file=sys.stderr)

    if verbose:
        print("INFO: Calculating the median and errors for each spectrum ...",
                file=sys.stderr)
    medians, errors = calc_median_errors(cc_results)
    if verbose:
        print("INFO: Writing the crosstalk-corrected spectra " + \
                "with estimated statistical errors ...",
                file=sys.stderr)
    for spec, data, err in zip(cc_spectra, medians, errors):
        spec.set_data(data, group_squeeze=group_squeeze)
        spec.add_stat_err(err, group_squeeze=group_squeeze)
        spec.write(clobber=clobber)
# main_crosstalk routine }}}


if __name__ == "__main__":
    # arguments' default values
    default_mode = "both"
    default_mc_times = 5000
    # commandline arguments parser
    parser = argparse.ArgumentParser(
            description="Correct the crosstalk effects for XMM EPIC spectra",
            epilog="Version: %s (%s)" % (__version__, __date__))
    parser.add_argument("config", help="config file in which describes " +\
            "the crosstalk relations ('ConfigObj' syntax)")
    parser.add_argument("-m", "--mode", dest="mode", default=default_mode,
            help="operation mode (both | crosstalk | deprojection)")
    parser.add_argument("-B", "--no-subtract-bkg", dest="subtract_bkg",
            action="store_false", help="do NOT subtract background first")
    parser.add_argument("-N", "--fix-negative", dest="fix_negative",
            action="store_true", help="fix negative channel values")
    parser.add_argument("-M", "--mc-times", dest="mc_times",
            type=int, default=default_mc_times,
            help="Monte Carlo times for error estimation")
    parser.add_argument("-C", "--clobber", dest="clobber",
            action="store_true", help="overwrite output file if exists")
    parser.add_argument("-v", "--verbose", dest="verbose",
            action="store_true", help="show verbose information")
    args = parser.parse_args()
    # merge commandline arguments and config
    config       = ConfigObj(args.config)
    subtract_bkg = set_argument("subtract_bkg", True,  args, config)
    fix_negative = set_argument("fix_negative", False, args, config)
    verbose      = set_argument("verbose",      False, args, config)
    clobber      = set_argument("clobber",      False, args, config)
    # operation mode
    mode = config.get("mode", default_mode)
    if args.mode != default_mode:
        mode = args.mode
    # Monte Carlo times
    mc_times = config.as_int("mc_times")
    if args.mc_times != default_mc_times:
        mc_times = args.mc_times

    if mode.lower() == "both":
        print("MODE: CROSSTALK + DEPROJECTION", file=sys.stderr)
        main(config, subtract_bkg=subtract_bkg, fix_negative=fix_negative,
                mc_times=mc_times, verbose=verbose, clobber=clobber)
    elif mode.lower() == "deprojection":
        print("MODE: DEPROJECTION", file=sys.stderr)
        main_deprojection(config, mc_times=mc_times,
                verbose=verbose, clobber=clobber)
    elif mode.lower() == "crosstalk":
        print("MODE: CROSSTALK", file=sys.stderr)
        main_crosstalk(config, subtract_bkg=subtract_bkg,
                fix_negative=fix_negative, mc_times=mc_times,
                verbose=verbose, clobber=clobber)
    else:
        raise ValueError("Invalid operation mode: %s" % mode)
    print(WARNING)

#  vim: set ts=4 sw=4 tw=0 fenc=utf-8 ft=python: #