aboutsummaryrefslogtreecommitdiffstats
path: root/python
diff options
context:
space:
mode:
authorAaron LI <aaronly.me@outlook.com>2016-04-19 15:00:03 +0800
committerAaron LI <aaronly.me@outlook.com>2016-04-19 15:00:03 +0800
commitac62b7b58d39728f051c7fad7e89b9adbffe9a81 (patch)
tree06f5ed4c74f86243df2ae0265b2ad2c7353348c4 /python
parent2bacc105e03231940b87c3f6ddba6ec3675b60fc (diff)
downloadatoolbox-ac62b7b58d39728f051c7fad7e89b9adbffe9a81.tar.bz2
Major update and rewrite to correct_crosstalk.py/crosstalk_deprojection.py
* Implement deprojection function * Implement Monte Carlo function for spectral error estimation
Diffstat (limited to 'python')
-rwxr-xr-xpython/correct_crosstalk.py750
-rwxr-xr-xpython/crosstalk_deprojection.py1726
2 files changed, 1726 insertions, 750 deletions
diff --git a/python/correct_crosstalk.py b/python/correct_crosstalk.py
deleted file mode 100755
index 7d83912..0000000
--- a/python/correct_crosstalk.py
+++ /dev/null
@@ -1,750 +0,0 @@
-#!/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] CIAO: Auxiliary Response File
-# http://cxc.harvard.edu/ciao/dictionary/arf.html
-# [3] CIAO: Redistribution Matrix File
-# http://cxc.harvard.edu/ciao/dictionary/rmf.html
-# [4] astropy - FITS format code
-# http://docs.astropy.org/en/stable/io/fits/usage/table.html#column-creation
-# [5] XSPEC - Spectral Fitting
-# https://heasarc.gsfc.nasa.gov/docs/xanadu/xspec/manual/XspecSpectralFitting.html
-#
-#
-# Weitian LI
-# Created: 2016-03-26
-# Updated: 2016-04-06
-#
-# ChangeLog:
-# 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:
-# * SpectrumSet: estimate channel errors by Monte Carlo simulations
-#
-# TODO:
-# * Spectrum: implement the grouping function (and quality columns)
-# * Split classes ARF, RMF, Spectrum, and SpectrumSet to a separate module
-#
-
-__version__ = "0.3.0"
-__date__ = "2016-04-02"
-
-
-"""
-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'.
-
-
-Sample config file (in `ConfigObj' syntax):
------------------------------------------------------------
-verbose = True
-clobber = False
-# whether to subtract the background before crosstalk correction
-subtract_bkg = True
-# whether to fix the negative channel values due to spectral subtractions
-fix_negative = True
-
-[...]
-...
-
-[reg2]
-outfile = cc_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
-
-[...]
-...
------------------------------------------------------------
-"""
-
-
-import numpy as np
-import scipy as sp
-import scipy.interpolate
-from astropy.io import fits
-from configobj import ConfigObj
-
-import sys
-import os
-import argparse
-from datetime import datetime
-
-
-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
-
- 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, copy=True):
- if copy:
- return self.specresp.copy()
- else:
- return self.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: ARF interpolating (this may take a while) ...",
- 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
-# 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
- # several important keywords
- EXPOSURE = None
- BACKSCAL = None
- RESPFILE = None
- ANCRFILE = None
- BACKFILE = None
- # numpy dtype and FITS format code of the spectrum data
- spec_dtype = None
- spec_fits_format = None
-
- def __init__(self, filename):
- 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
- # 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")
-
- def get_data(self, copy=True):
- if copy:
- return self.spec_data.copy()
- else:
- return self.spec_data
-
- def get_channel(self, copy=True):
- if copy:
- return self.channel.copy()
- else:
- return self.channel
-
- def reset_header_keywords(self,
- keywords=["ANCRFILE", "RESPFILE", "BACKFILE"]):
- """
- Reset the keywords to "NONE" to avoid confusion or mistakes.
- """
- for kw in keywords:
- if kw in self.header:
- header[kw] = "NONE"
-
- def write(self, filename, clobber=False):
- """
- Create a new "SPECTRUM" table/extension and replace the original
- one, then write to output file.
- """
- ext_spec_cols = fits.ColDefs([
- 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)])
- 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
-
- # 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"
-
- def __init__(self, filename, arffile=None, rmffile=None, bkgfile=None):
- super(self.__class__, self).__init__(filename)
- # 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 arffile is not None:
- self.arf = ARF(arffile)
- if rmffile is not None:
- self.rmf = RMF(rmffile)
- if bkgfile is not None:
- self.bkg = Spectrum(bkgfile)
-
- 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", copy=True):
- """
- Get the corresponding ARF curve data if the ARF presents.
-
- Return:
- [ energy, resp ]
- where the "energy" and "resp" are the ARF energy values and
- spectral response array, respectively.
- """
- if self.arf is None:
- return None
- else:
- energy = self.arf.get_energy(mean=mean)
- resp = self.arf.get_data(copy=copy)
- return [ energy, resp ]
-
- def subtract_bkg(self, inplace=True, 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.
-
- 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
- # also record history
- self.header.add_history(operation)
- return spec_data_subbkg
-
- def subtract(self, spectrumset, cross_arf, 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()
- arfresp_spec = spectrumset.arf.interpolate(x=energy, verbose=verbose)
- arfresp_cross = cross_arf.interpolate(x=energy, verbose=verbose)
- arf_ratio = arfresp_cross / arfresp_spec
- # fix nan values due to division by zero
- arf_ratio[np.isnan(arf_ratio)] = 0.0
- self.spec_data -= spectrumset.get_data() * arf_ratio
- # record history
- self.header.add_history(operation)
-
- def compensate(self, cross_arf, 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()
- arfresp_this = self.arf.interpolate(x=energy, verbose=verbose)
- arfresp_cross = cross_arf.interpolate(x=energy, verbose=verbose)
- arf_ratio = arfresp_cross / arfresp_this
- # fix nan values due to division by zero
- arf_ratio[np.isnan(arf_ratio)] = 0.0
- self.spec_data += self.get_data() * arf_ratio
- # record history
- self.header.add_history(operation)
-
- def fix_negative(self, 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
- self.header.add_history(" FIXED NEGATIVE CHANNELS")
-# class SpectrumSet }}}
-
-
-class Crosstalk: # {{{
- """
- Crosstalk 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
- # output filename to which write the corrected spectrum
- outfile = None
-
- def __init__(self, config):
- """
- 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"],
- arffile=config["arf"], rmffile=config.get("rmf"),
- bkgfile=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"],
- arffile=reg_in["arf"], rmffile=reg_in.get("rmf"),
- bkgfile=reg_in.get("bkg"))
- self.cross_in_specset.append(specset)
- self.cross_in_arf.append(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(arffile))
- # output filename
- self.outfile = config["outfile"]
-
- def do_correction(self, subtract_bkg=True, fix_negative=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.
- """
- 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, verbose=verbose)
- # also apply background subtraction to the surrounding spectra
- for specset in self.cross_in_specset:
- specset.subtract_bkg(inplace=True, 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, 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,
- 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(verbose=verbose)
- self.spectrumset.header.add_history("END Crosstalk Correction")
-
- def write(self, filename=None, clobber=False):
- if filename is None:
- filename = self.outfile
- self.spectrumset.reset_header_keywords(
- keywords=["ANCRFILE", "BACKFILE"])
- self.spectrumset.write(filename, clobber=clobber)
-# class Crosstalk }}}
-
-
-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
-
-
-def main():
- 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("-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("-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()
-
- 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)
-
- for region in config.sections:
- if verbose:
- print("INFO: processing '%s' ..." % region, file=sys.stderr)
- crosstalk = Crosstalk(config.get(region))
- crosstalk.do_correction(subtract_bkg=subtract_bkg,
- fix_negative=fix_negative, verbose=verbose)
- crosstalk.write(clobber=clobber)
-
-
-if __name__ == "__main__":
- main()
-
-# vim: set ts=4 sw=4 tw=0 fenc=utf-8 ft=python: #
diff --git a/python/crosstalk_deprojection.py b/python/crosstalk_deprojection.py
new file mode 100755
index 0000000..dd099b4
--- /dev/null
+++ b/python/crosstalk_deprojection.py
@@ -0,0 +1,1726 @@
+#!/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] CIAO: Auxiliary Response File
+# http://cxc.harvard.edu/ciao/dictionary/arf.html
+# [3] CIAO: Redistribution Matrix File
+# http://cxc.harvard.edu/ciao/dictionary/rmf.html
+# [4] astropy - FITS format code
+# http://docs.astropy.org/en/stable/io/fits/usage/table.html#column-creation
+# [5] XSPEC - Spectral Fitting
+# https://heasarc.gsfc.nasa.gov/docs/xanadu/xspec/manual/XspecSpectralFitting.html
+# [6] 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-04-19
+#
+# ChangeLog:
+# 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.0"
+__date__ = "2016-04-19"
+
+
+"""
+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
+
+[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
+ 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 reset_header_keywords(self,
+ keywords=["ANCRFILE", "RESPFILE", "BACKFILE"]):
+ """
+ Reset the keywords to "NONE" to avoid confusion or mistakes.
+ """
+ for kw in keywords:
+ if kw in self.header:
+ self.header[kw] = "NONE"
+
+ 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, 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
+ # also record history
+ self.header.add_history(operation)
+ return spec_data_subbkg
+
+ def subtract(self, spectrumset, cross_arf, groupped=False,
+ group_squeeze=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
+ self.header.add_history(operation)
+
+ def compensate(self, cross_arf, groupped=False, group_squeeze=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
+ self.header.add_history(operation)
+
+ def fix_negative(self, 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
+ 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, 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.
+ """
+ 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, verbose=verbose)
+ # also apply background subtraction to the surrounding spectra
+ for specset in self.cross_in_specset:
+ specset.subtract_bkg(inplace=True, 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, 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,
+ 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(verbose=verbose)
+ self.spectrumset.header.add_history("END Crosstalk Correction")
+ # reset header keywords
+ self.spectrumset.reset_header_keywords(
+ keywords=["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, verbose=True):
+ #
+ # 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
+ 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_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 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 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.
+
+ TODO: investigate the errors calculation approach used here!
+ """
+ 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,
+ 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(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,
+ 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(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.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(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(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.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,
+ verbose=verbose)
+ 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,
+ 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: #