From c1362d46a2ad60c42fd7b2c94a0aba026d54c7db Mon Sep 17 00:00:00 2001
From: Aaron LI <aaronly.me@outlook.com>
Date: Sat, 2 Apr 2016 12:34:11 +0800
Subject: correct_crosstalk.py: interpolate ARF when correcting crosstalk
---
python/correct_crosstalk.py | 266 ++++++++++++++++++++++++++++----------------
1 file changed, 169 insertions(+), 97 deletions(-)
(limited to 'python')
diff --git a/python/correct_crosstalk.py b/python/correct_crosstalk.py
index c514504..c3e9815 100755
--- a/python/correct_crosstalk.py
+++ b/python/correct_crosstalk.py
@@ -2,28 +2,46 @@
# -*- coding: utf-8 -*-
#
# References:
-# [?] astropy - FITS format code
+# [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
-# [?] XSPEC - Spectral Fitting
+# [5] XSPEC - Spectral Fitting
# https://heasarc.gsfc.nasa.gov/docs/xanadu/xspec/manual/XspecSpectralFitting.html
#
#
# Weitian LI
# Created: 2016-03-26
-# Updated: 2016-04-01
+# Updated: 2016-04-02
#
# ChangeLog:
+# 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
#
-# 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
+# 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"
"""
@@ -37,6 +55,8 @@ 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
@@ -73,6 +93,7 @@ bkg = reg2_bkg.pi
import numpy as np
import scipy as sp
+import scipy.interpolate
from astropy.io import fits
from configobj import ConfigObj
@@ -92,6 +113,9 @@ class ARF: # {{{
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
@@ -115,6 +139,8 @@ class ARF: # {{{
energ_lo = None
energ_hi = None
specresp = None
+ # function of the interpolated ARF
+ f_interp = None
def __init__(self, filename):
self.filename = filename
@@ -162,17 +188,19 @@ class ARF: # {{{
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 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 f_interp(x)
+ return self.f_interp(x)
else:
- return f_interp
+ return self.f_interp
# class ARF }}}
@@ -194,13 +222,13 @@ class RMF: # {{{
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.
+ + 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
@@ -302,9 +330,6 @@ class RMF: # {{{
class Spectrum: # {{{
"""
Class that deals with the X-ray spectrum file (usually *.pi).
-
- TODO:
- * to implement the grouping function (and quality columns)
"""
filename = None
# FITS object return by `fits.open()'
@@ -400,10 +425,6 @@ class SpectrumSet(Spectrum): # {{{
This class handles a set of spectrum, including the source spectrum,
RMF, ARF, and the background spectrum.
- TODO:
- * Subtract background spectrum before correcting crosstalk effects
- * Estimate channel errors by Monte Carlo simulations
-
**NOTE**:
The "COUNTS" column data are converted from "int32" to "float32",
since this spectrum will be subtracted/compensated according to the
@@ -438,22 +459,35 @@ class SpectrumSet(Spectrum): # {{{
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, copy=True):
+ def get_arf(self, mean="geometric", copy=True):
"""
- Get the corresponding ARF curve data for this spectrum.
+ 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:
- return self.arf.get_data(copy=copy)
+ energy = self.arf.get_energy(mean=mean)
+ resp = self.arf.get_data(copy=copy)
+ return [ energy, resp ]
- def subtract_bkg(self, inplace=True):
+ 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
@@ -471,13 +505,19 @@ class SpectrumSet(Spectrum): # {{{
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, spectrum, cross_arf, verbose=False):
+ 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.
@@ -487,17 +527,23 @@ class SpectrumSet(Spectrum): # {{{
correction, as well as the below `compensate()' procedure.
NOTE:
- The crosstalk ARF must be provided, since the `spectrum.arf' is
- required to be its ARF without taking crosstalk into account:
- spec1_new = spec1 - spec2 * (cross_arf_2_to_1 / arf2)
+ 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, spectrum.arf.filename, spectrum.filename)
+ cross_arf.filename, spectrumset.arf.filename,
+ spectrumset.filename)
if verbose:
print(operation, file=sys.stderr)
- arf_ratio = cross_arf.get_data() / spectrum.get_arf()
+ 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 -= spectrum.get_data() * arf_ratio
+ self.spec_data -= spectrumset.get_data() * arf_ratio
# record history
self.header.add_history(operation)
@@ -513,7 +559,11 @@ class SpectrumSet(Spectrum): # {{{
cross_arf.filename, self.arf.filename, self.filename)
if verbose:
print(operation, file=sys.stderr)
- arf_ratio = cross_arf.get_data() / self.get_arf()
+ 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
@@ -534,8 +584,11 @@ class SpectrumSet(Spectrum): # {{{
while len(neg_channels) > 0:
i += 1
if verbose:
- print("*** Fixing negative channels: iteration %d ..." % i,
- file=sys.stderr)
+ 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:
@@ -548,22 +601,24 @@ class SpectrumSet(Spectrum): # {{{
neg_counts = self.spec_data < 0
neg_channels = np.arange(N, dtype=np.int)[neg_counts]
if i > 0:
- print("*** Fixed negative channels ***", file=sys.stderr)
-# class Spectrum }}}
+ print("FIXED ***", file=sys.stderr)
+ # record history
+ self.header.add_history(" FIXED NEGATIVE CHANNELS")
+# class SpectrumSet }}}
class Crosstalk: # {{{
"""
Crosstalk correction.
"""
- # `Spectrum' object for the spectrum to be corrected
- spectrum = None
- # XXX: do NOT use list (e.g., []) here, otherwise, all the instances
- # will share these list properties.
- # `Spectrum' and `ARF' objects corresponding to the spectra from which
- # the photons were scattered into this spectrum.
- cross_in_spec = None
- cross_in_arf = None
+ # `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
@@ -575,17 +630,19 @@ class Crosstalk: # {{{
Arguments:
* config: a section of the whole config file (`ConfigObj' object)
"""
- self.cross_in_spec = []
- self.cross_in_arf = []
- self.cross_out_arf = []
+ self.cross_in_specset = []
+ self.cross_in_arf = []
+ self.cross_out_arf = []
# this spectrum to be corrected
- self.spectrum = Spectrum(filename=config["spec"],
- arffile=config["arf"], rmffile=config["rmf"],
- bkgfile=config["bkg"])
+ 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():
- spec = Spectrum(reg_in["spec"], reg_in["arf"])
- self.cross_in_spec.append(spec)
+ 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:
@@ -595,45 +652,72 @@ class Crosstalk: # {{{
# output filename
self.outfile = config["outfile"]
- def do_correction(self, fix_negative=False, verbose=False):
+ def do_correction(self, subtract_bkg=True, fix_negative=False,
+ verbose=False):
"""
- Perform the crosstalk correction.
+ 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.spectrum.header.add_history("Crosstalk Correction BEGIN")
- self.spectrum.header.add_history(" TOOL: %s @ %s" % (\
- os.path.basename(sys.argv[0]), datetime.utcnow().isoformat()))
+ 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 spec, cross_arf in zip(self.cross_in_spec, self.cross_in_arf):
- self.spectrum.subtract(spectrum=spec, cross_arf=cross_arf,
- verbose=verbose)
+ 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.spectrum.compensate(cross_arf=cross_arf, verbose=verbose)
+ 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.spectrum.fix_negative(verbose=verbose)
- self.spectrum.header.add_history("END Crosstalk Correction")
+ 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.spectrum.reset_header_keywords(
- keywords=["ANCRFILE", "RESPFILE", "BACKFILE"])
- self.spectrum.write(filename, clobber=clobber)
+ 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 of XMM spectra")
+ 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)")
+ "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",
@@ -644,29 +728,17 @@ def main():
config = ConfigObj(args.config)
- fix_negative = False
- if "fix_negative" in config.keys():
- fix_negative = config.as_bool("fix_negative")
- if args.fix_negative:
- fix_negative = args.fix_negative
-
- verbose = False
- if "verbose" in config.keys():
- verbose = config.as_bool("verbose")
- if args.verbose:
- verbose = args.verbose
-
- clobber = False
- if "clobber" in config.keys():
- clobber = config.as_bool("clobber")
- if args.clobber:
- clobber = args.clobber
+ 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(fix_negative=fix_negative, verbose=verbose)
+ crosstalk.do_correction(subtract_bkg=subtract_bkg,
+ fix_negative=fix_negative, verbose=verbose)
crosstalk.write(clobber=clobber)
--
cgit v1.2.2