diff options
| -rwxr-xr-x | astro/spectrum/crosstalk_deprojection.py | 196 |
1 files changed, 104 insertions, 92 deletions
diff --git a/astro/spectrum/crosstalk_deprojection.py b/astro/spectrum/crosstalk_deprojection.py index 263cb15..617b8ff 100755 --- a/astro/spectrum/crosstalk_deprojection.py +++ b/astro/spectrum/crosstalk_deprojection.py @@ -1,29 +1,10 @@ #!/usr/bin/env python3 -# -*- coding: utf-8 -*- # -# References: -# [1] Definition of RMF and ARF file formats -# https://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/docs/memos/cal_gen_92_002/cal_gen_92_002.html -# [2] The OGIP Spectral File Format -# https://heasarc.gsfc.nasa.gov/docs/heasarc/ofwg/docs/summary/ogip_92_007_summary.html -# [3] CIAO: Auxiliary Response File -# http://cxc.harvard.edu/ciao/dictionary/arf.html -# [4] CIAO: Redistribution Matrix File -# http://cxc.harvard.edu/ciao/dictionary/rmf.html -# [5] astropy - FITS format code -# http://docs.astropy.org/en/stable/io/fits/usage/table.html#column-creation -# [6] XSPEC - Spectral Fitting -# https://heasarc.gsfc.nasa.gov/docs/xanadu/xspec/manual/XspecSpectralFitting.html -# [7] Direct X-ray Spectra Deprojection -# https://www-xray.ast.cam.ac.uk/papers/dsdeproj/ -# Sanders & Fabian 2007, MNRAS, 381, 1381 -# -# -# Weitian LI -# Created: 2016-03-26 -# Updated: 2016-06-07 +# Copyright (c) 2016-2017 Weitian LI +# 2016-03-26 # # Change log: +# 2017-11-12: Cleanups # 2016-06-07: # * Explain the errors/uncertainties calculation approach # 2016-04-20: @@ -60,16 +41,11 @@ # * Implement background subtraction # * Add config `subtract_bkg' and corresponding argument # -# XXX/FIXME: -# * Deprojection: account for ARF differences across different regions -# # TODO: +# * Deprojection: account for ARF differences across different regions # * Split classes ARF, RMF, Spectrum, and SpectrumSet to a separate module # -__version__ = "0.5.3" -__date__ = "2016-06-07" - """ Correct the crosstalk effect of XMM spectra by subtracting the photons @@ -82,6 +58,25 @@ to deproject the crosstalk-corrected spectra to derive the spectra with both the crosstalk effect and projection effect corrected. +References +---------- +[1] Definition of RMF and ARF file formats + https://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/docs/memos/cal_gen_92_002/cal_gen_92_002.html +[2] The OGIP Spectral File Format + https://heasarc.gsfc.nasa.gov/docs/heasarc/ofwg/docs/summary/ogip_92_007_summary.html +[3] CIAO: Auxiliary Response File + http://cxc.harvard.edu/ciao/dictionary/arf.html +[4] CIAO: Redistribution Matrix File + http://cxc.harvard.edu/ciao/dictionary/rmf.html +[5] astropy - FITS format code + http://docs.astropy.org/en/stable/io/fits/usage/table.html#column-creation +[6] XSPEC - Spectral Fitting + https://heasarc.gsfc.nasa.gov/docs/xanadu/xspec/manual/XspecSpectralFitting.html +[7] Direct X-ray Spectra Deprojection + https://www-xray.ast.cam.ac.uk/papers/dsdeproj/ + Sanders & Fabian 2007, MNRAS, 381, 1381 + + Sample config file (in `ConfigObj' syntax): ----------------------------------------------------------- # operation mode: deprojection, crosstalk, or both (default) @@ -135,14 +130,6 @@ bkg = reg2_bkg.pi ----------------------------------------------------------- """ -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 @@ -151,12 +138,23 @@ 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 +__version__ = "0.5.5" +__date__ = "2017-11-12" + +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. +****************************************************************************** +""" + + def group_data(data, grouping): """ Group the data with respect to the supplied `grouping' specification @@ -164,7 +162,7 @@ def group_data(data, grouping): 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() + data_grp = np.array(data) for i in reversed(range(len(data))): if grouping[i] == 1: # the beginning channel of a group @@ -186,7 +184,8 @@ class ARF: # {{{ 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**: + 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). @@ -200,7 +199,8 @@ class ARF: # {{{ spectral channels, i.e., the energy positions of each ARF data point and spectral channel are consistent. Thus the interpolation is not needed. - References: + References + ---------- [1] CIAO: Auxiliary Response File http://cxc.harvard.edu/ciao/dictionary/arf.html [2] Definition of RMF and ARF file formats @@ -279,9 +279,9 @@ class ARF: # {{{ 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", + print("INFO: interpolating ARF '%s' (may take a while) ..." % + self.filename, file=sys.stderr) + f_interp = scipy.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: @@ -331,7 +331,8 @@ class RMF: # {{{ variable length array columns. This format is compact but hard to manipulate and understand. - **CAVEAT/NOTE**: + 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, @@ -340,7 +341,8 @@ class RMF: # {{{ + The `ENERG_LO' and `ENERG_HI' columns of the "MATRIX" extension are the same as that of a ARF. - References: + References + ---------- [1] CIAO: Redistribution Matrix File http://cxc.harvard.edu/ciao/dictionary/rmf.html [2] Definition of RMF and ARF file formats @@ -412,8 +414,8 @@ class RMF: # {{{ 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) ]) + 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 @@ -636,7 +638,8 @@ class Spectrum: # {{{ Reset the keywords to "NONE" to avoid confusion or mistakes, and also add mandatory spectral keywords if missing. - Reference: + Reference + --------- [1] The OGIP Spectral File Format, Sec. 3.1.1 https://heasarc.gsfc.nasa.gov/docs/heasarc/ofwg/docs/summary/ogip_92_007_summary.html """ @@ -715,7 +718,8 @@ class SpectrumSet(Spectrum): # {{{ This class handles a set of spectrum, including the source spectrum, RMF, ARF, and the background spectrum. - **NOTE**: + NOTE + ---- The "COUNTS" column data are converted from "int32" to "float32", since this spectrum will be subtracted/compensated according to the ratios of ARFs. @@ -771,7 +775,8 @@ class SpectrumSet(Spectrum): # {{{ """ Get the energy values of each channel if RMF present. - NOTE: + 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" @@ -782,7 +787,7 @@ class SpectrumSet(Spectrum): # {{{ else: return self.rmf.get_energy(mean=mean) - def get_arf(self, mean="geometric", groupped=True, copy=True): + def get_arf(self, groupped=True, copy=True): """ Get the interpolated ARF data w.r.t the spectral channel energies if the ARF presents. @@ -835,7 +840,9 @@ class SpectrumSet(Spectrum): # {{{ 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 + NOTE + ---- + The spectral data and errors (i.e., `self.spec_data', and `self.spec_err') is MODIFIED! """ self.spec_data *= (360.0 / self.angle) @@ -853,7 +860,8 @@ class SpectrumSet(Spectrum): # {{{ spectrum, the ARF (which is used during the later spectral manipulations), and the background spectrum (if presents). - NOTE: + 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 @@ -904,11 +912,11 @@ class SpectrumSet(Spectrum): # {{{ 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) + 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() @@ -937,7 +945,8 @@ class SpectrumSet(Spectrum): # {{{ both be subtracted before applying this subtraction for crosstalk correction, as well as the below `compensate()' procedure. - NOTE: + 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) @@ -966,8 +975,8 @@ class SpectrumSet(Spectrum): # {{{ 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 + spec_data = (self.get_data(group_squeeze=group_squeeze) - + spectrumset.get_data(group_squeeze=group_squeeze)*arf_ratio) self.set_data(spec_data, group_squeeze=group_squeeze) # record history if add_history: @@ -1001,8 +1010,8 @@ class SpectrumSet(Spectrum): # {{{ 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 + spec_data = (self.get_data(group_squeeze=group_squeeze) + + self.get_data(group_squeeze=group_squeeze) * arf_ratio) self.set_data(spec_data, group_squeeze=group_squeeze) # record history if add_history: @@ -1017,8 +1026,8 @@ class SpectrumSet(Spectrum): # {{{ N = len(neg_counts) neg_channels = np.arange(N, dtype=int)[neg_counts] if len(neg_channels) > 0: - print("WARNING: %d channels have NEGATIVE counts" % \ - len(neg_channels), file=sys.stderr) + print("WARNING: %d channels have NEGATIVE counts" % + len(neg_channels), file=sys.stderr) i = 0 while len(neg_channels) > 0: i += 1 @@ -1067,7 +1076,9 @@ class SpectrumSet(Spectrum): # {{{ according to the estimated spectral group errors by assuming the normal distribution. - NOTE: this method should be called AFTER the `copy()' method. + NOTE + ---- + This method should be called *after* the `copy()' method. """ super().randomize() if self.bkg: @@ -1158,7 +1169,7 @@ class Crosstalk: # {{{ """ if add_history: self.spectrumset.header.add_history("Crosstalk Correction BEGIN") - self.spectrumset.header.add_history(" TOOL: %s (v%s) @ %s" % (\ + self.spectrumset.header.add_history(" TOOL: %s (v%s) @ %s" % ( os.path.basename(sys.argv[0]), __version__, datetime.utcnow().isoformat())) # background subtraction @@ -1205,8 +1216,8 @@ class Crosstalk: # {{{ 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 ] + new.cross_in_specset = [specset.copy() + for specset in self.cross_in_specset] return new def randomize(self): @@ -1232,12 +1243,14 @@ class Deprojection: # {{{ assumption of spherical symmetry of the source object, and produce the DEPROJECTED spectra. - NOTE: + NOTE + ---- * Assumption of the spherical symmetry * Background should be subtracted before deprojection * ARF differences of different regions are taken into account - Reference & Credit: + Credit + ------ [1] Direct X-ray Spectra Deprojection https://www-xray.ast.cam.ac.uk/papers/dsdeproj/ Sanders & Fabian 2007, MNRAS, 381, 1381 @@ -1270,9 +1283,8 @@ class Deprojection: # {{{ 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) + 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:]) @@ -1280,8 +1292,8 @@ class Deprojection: # {{{ def subtract_bkg(self, verbose=True): for spec in self.spectra: if not spec.bkg: - raise ValueError("Spectrum '%s' has NO background" % \ - spec.filename) + raise ValueError("Spectrum '%s' has NO background" % + spec.filename) spec.subtract_bkg(inplace=True, verbose=verbose) def apply_grouping(self, verbose=False): @@ -1423,9 +1435,10 @@ def calc_median_errors(results): Calculate the median and errors for the spectral data gathered through Monte Carlo simulations. - NOTE: - Errors calculation just use the quantiles, - i.e., 1sigma ~= 68.3% = Q(84.15%) - Q(15.85%) + NOTE + ---- + Errors are calculated from the percentiles, + i.e., 1sigma ~= 68.3% = P(84.15%) - P(15.85%) """ results = np.array(results) # `results' now has shape: (mc_times, num_spec, num_channel) @@ -1520,8 +1533,8 @@ def main(config, subtract_bkg, fix_negative, mc_times, 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) + 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) @@ -1550,7 +1563,7 @@ def main(config, subtract_bkg, fix_negative, mc_times, 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 " + \ + print("INFO: Writing the crosstalk-corrected and deprojected " + "spectra with estimated statistical errors ...", file=sys.stderr) deprojection.fix_header() deprojection.write(clobber=clobber) @@ -1619,8 +1632,8 @@ def main_deprojection(config, mc_times, verbose=False, clobber=False): 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) + 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) @@ -1643,7 +1656,7 @@ def main_deprojection(config, mc_times, verbose=False, clobber=False): 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 " + \ + print("INFO: Writing the deprojected spectra " + "with estimated statistical errors ...", file=sys.stderr) deprojection.fix_header() deprojection.write(clobber=clobber) @@ -1719,8 +1732,8 @@ def main_crosstalk(config, subtract_bkg, fix_negative, mc_times, 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) + 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) @@ -1743,9 +1756,8 @@ def main_crosstalk(config, subtract_bkg, fix_negative, mc_times, 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) + 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) @@ -1761,8 +1773,8 @@ if __name__ == "__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("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", |