aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rwxr-xr-xastro/spectrum/crosstalk_deprojection.py196
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",