diff options
Diffstat (limited to 'python/correct_crosstalk.py')
| -rwxr-xr-x | python/correct_crosstalk.py | 340 |
1 files changed, 299 insertions, 41 deletions
diff --git a/python/correct_crosstalk.py b/python/correct_crosstalk.py index 4c7f820..9f86763 100755 --- a/python/correct_crosstalk.py +++ b/python/correct_crosstalk.py @@ -1,40 +1,66 @@ #!/usr/bin/env python3 -# -# Correct the crosstalk effect of XMM spectra by subtracting the -# scattered photons from surrounding regions, and by compensating -# the photons scattered to surrounding regions, according to the -# generated crosstalk ARFs. -# -# Sample config file (in `ConfigObj' syntax): -#---------------------- -# fix_negative = True -# verbose = True -# clobber = False -# -# [reg2] -# outfile = cc_reg2.pi -# spec = reg2.pi -# arf = reg2.arf -# [[cross_in]] -# [[[in1]]] -# spec = reg1.pi -# arf = reg1.arf -# cross_arf = reg_1-2.arf -# [[[in2]]] -# spec = reg3.pi -# arf = reg3.arf -# cross_arf = reg_3-2.arf -# [[cross_out]] -# cross_arf = reg_2-1.arf, reg_2-3.arf -#---------------------- +# -*- coding: utf-8 -*- # # Weitian LI # Created: 2016-03-26 -# Updated: 2016-03-28 +# Updated: 2016-04-01 +# +# ChangeLog: +# 2016-04-01: +# * Greatly update the documentations (e.g., description, sample config) +# * Add class `RMF' +# * Add method `get_energy()' for class `ARF' +# +# TODO/XXX: +# * Add background spectrum fields +# * Subtract background spectrum before correcting crosstalk effects +# * Estimate channel errors by Monte Carlo simulations +# * Split classes ARF, RMF, Spectrum to a separate module # -from astropy.io import fits + +""" +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 fix the negative channel values due to spectral subtractions +fix_negative = True + +[...] +... + +[reg2] +outfile = cc_reg2.pi +spec = reg2.pi +arf = reg2.arf + [[cross_in]] + [[[in1]]] + spec = reg1.pi + arf = reg1.arf + cross_arf = reg_1-2.arf + [[[in2]]] + spec = reg3.pi + arf = reg3.arf + 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 +from astropy.io import fits from configobj import ConfigObj import sys @@ -43,22 +69,45 @@ import argparse from datetime import datetime -class ARF: +class ARF: # {{{ """ - Deal with X-ray ARF file (.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**: + 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 - header = None + fitsobj = None + # only consider the "SPECTRUM" extension + header = None energ_lo = None energ_hi = None specresp = None def __init__(self, filename): self.filename = filename - self.fitsobj = fits.open(filename) - ext_specresp = self.fitsobj["SPECRESP"] - self.header = ext_specresp.header + 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"] @@ -69,27 +118,202 @@ class ARF: else: return self.specresp + def get_energy(self, mean="geometric"): + """ + Return the mean energy values of the ARF. -class Spectrum: + 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. + """ + 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) + if x is not None: + return f_interp(x) + else: + return 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 the above ARF CAVEAT/NOTE for the XMM EPIC pn. + + This class (currently) only deals with the "EBOUNDS" extension, which + contains the `CHANNEL', `E_MIN' and `E_MAX' columns. This `CHANNEL' + is the same as that of a spectrum. Therefore, the energy values can + be used to interpolate and extrapolate the ARF curve. + + The `ENERG_LO' and `ENERG_HI' columns of "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): + f_chan -= 1 # FITS indices are 1-based + 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 range(n_energy): + 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: # {{{ """ Deal with X-ray spectrum (.pi) + TODO/XXX: + * Add background spectrum fields + * Subtract background spectrum before correcting crosstalk effects + * Strip keywords ANCRFILE, RESPFILE, BACKFILE from the output spectrum + * Estimate channel errors by Monte Carlo simulations + NOTE: The "COUNTS" column data are converted from "int32" to "float32". """ filename = None - # FITS object return by `fits.open' + # FITS object return by `fits.open()' fitsobj = None # header of "SPECTRUM" extension header = None # "SPECTRUM" extension data channel = None - # name of the column containing the spectrum data, either "COUNTS" or "RATE" + # name of the column containing the spectrum data ("COUNTS" or "RATE") spec_colname = None # spectrum data spec_data = None # ARF object for this spectrum arf = None + # RMF object for this spectrum + rmf = None def __init__(self, filename, arffile): self.filename = filename @@ -114,7 +338,25 @@ class Spectrum: else: return self.spec_data + def get_channel(self, copy=True): + if copy: + return self.channel.copy() + else: + return self.channel + + def get_energy(self, mean="geometric"): + """ + Get the energy values of each channel if RMF present. + """ + if self.rmf is None: + return None + else: + return self.rmf.get_energy(mean=mean) + def get_arf(self, copy=True): + """ + Get the corresponding ARF curve data for this spectrum. + """ if self.arf is None: return None else: @@ -189,6 +431,15 @@ class Spectrum: if i > 0: print("*** Fixed negative channels ***", file=sys.stderr) + 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 @@ -202,9 +453,10 @@ class Spectrum: header=self.header) self.fitsobj["SPECTRUM"] = ext_spec self.fitsobj.writeto(filename, clobber=clobber, checksum=True) +# class Spectrum }}} -class Crosstalk: +class Crosstalk: # {{{ """ Crosstalk correction. """ @@ -245,6 +497,9 @@ class Crosstalk: self.outfile = config["outfile"] def do_correction(self, fix_negative=False, verbose=False): + """ + Perform the crosstalk correction. + """ self.spectrum.header.add_history("Crosstalk Correction BEGIN") self.spectrum.header.add_history(" TOOL: %s @ %s" % (\ os.path.basename(sys.argv[0]), datetime.utcnow().isoformat())) @@ -269,7 +524,10 @@ class Crosstalk: def write(self, filename=None, clobber=False): if filename is None: filename = self.outfile + self.spectrum.reset_header_keywords( + keywords=["ANCRFILE", "RESPFILE", "BACKFILE"]) self.spectrum.write(filename, clobber=clobber) +# class Crosstalk }}} def main(): |