Skip to content

Commit

Permalink
Fix handling of ROSAT RMF file with AstroPy
Browse files Browse the repository at this point in the history
This commit adds tests and fixes the code handling the ROSAT RMF
with the AstroPy I/O backend. There is a known issue with the
CIAO Crates library in CIAO 4.9 which means that the ROSAT RMF
can not be read in using the Crates I/O backend with Python 3.5
(it works okay with Python 2.7). The tests are therefore setup to
skip certain combinations of Python and backend.

The test values used as the "truth" here were calculated using
CIAO 4.9/Python 2.7 and validated against XSPEC 12.9.1p.

It leverages a fixture for cleaning the Sherpa ui state from
changeset 0d44a3e

An explanation of the fix for handling non-variable-length
RMF files with the AstroPy backend:

PR sherpa#358 added AstroPy support for RMF files where the data was not
encoded using FITS variable-length arrays for the matrix. This means
that the data can not just be concatenated together. Instead each
row in the matrix (which, in this case, is a regular nchannels by
nbins array) has to be subset (maybe multiple times).

The commit I added to do this - cfad850
worked in the Swift case, where the F_CHAN/N_CHAN values were set up
to effectively remove the "space saving" mechanism they provide (that
is, all elements of each row in the MATRIX column are used). It does
not work for cases where F_CHAN/N_CHAN are not set to 0/nbins for
each row. The code was embarassingly-confused over what it was meant
to be doing, and only worked for the Swift case because of the data's
simplified structure.
  • Loading branch information
DougBurke authored and olaurino committed Feb 9, 2018
1 parent e37c21d commit 649d733
Show file tree
Hide file tree
Showing 3 changed files with 516 additions and 29 deletions.
2 changes: 1 addition & 1 deletion sherpa-test-data
50 changes: 22 additions & 28 deletions sherpa/astro/io/pyfits_backend.py
Original file line number Diff line number Diff line change
Expand Up @@ -697,12 +697,14 @@ def get_rmf_data(arg, make_copy=False):
dtype=SherpaUInt)
data['n_chan'] = _require_vec(hdu, 'N_CHAN', fix_type=True,
dtype=SherpaUInt)

# Read MATRIX as-is -- we will flatten it below, because
# we need to remove all rows corresponding to n_grp[row] == 0
#
# TODO: I would expect this to error out if MATRIX is not available
#
data['matrix'] = None
if 'MATRIX' not in hdu.columns.names:
pass
else:
if 'MATRIX' in hdu.columns.names:
data['matrix'] = hdu.data.field('MATRIX')

data['header'] = _get_meta_data(hdu)
Expand Down Expand Up @@ -739,23 +741,11 @@ def get_rmf_data(arg, make_copy=False):
finally:
rmf.close()

# ## For every row i of the response matrix, such that
# ## n_grp[i] == 0, we need to remove that row from the
# ## n_chan, f_chan, and matrix arrays we are constructing
# ## to be passed up to the DataRMF data structure.

# ## This is trivial for n_chan and f_chan. For the matrix
# ## array this can be more work -- can't just remove all
# ## zeroes, because some rows where n_grp[row] > 0 might
# ## still have zeroes in the matrix. I add new code first
# ## to deal with the matrix, then simpler code to remove zeroes
# ## from n_chan and f_chan.

# Read in MATRIX column with structure as-is -- i.e., as an array of
# arrays. Then flatten it, but include only those arrays that come from
# rows where n_grp[row] > 0. Zero elements can only be included from
# rows where n_grp[row] > 0. SMD 05/23/13

# Remove any rows from the MATRIX and F_CHAN/N_CHAN where N_GRP is 0,
# since they do not add any data. This is to match the Crates backend.
# Note that crates uses the sherpa.astro.utils.resp_init routine,
# but it's not clear why, so it is not used here for now.
#
good = (data['n_grp'] > 0)
data['matrix'] = data['matrix'][good]

Expand All @@ -779,15 +769,18 @@ def get_rmf_data(arg, make_copy=False):
n_chan = data['n_chan'][good]

rowdata = []
for i, (ng, fcs, ncs) in enumerate(zip(n_grp, f_chan, n_chan)):
mrow = matrix[i]
for mrow, ng, ncs in zip(matrix, n_grp, n_chan):
# Need a RMF which ng>1 to test this with.
if ng == 1:
fcs = [fcs]
ncs = [ncs]
for fc, nc in zip(fcs, ncs):
# It appears that F_CHAN is 0-based
#
rowdata.append(mrow[fc: fc + nc])
start = 0
for nc in ncs:
# n_chan can be an unsigned integer. Adding a Python
# integer to a NumPy unsigned integer appears to return
# a float.
end = start + numpy.int(nc)
rowdata.append(mrow[start:end])
start = end

data['matrix'] = numpy.concatenate(rowdata)

Expand All @@ -796,7 +789,6 @@ def get_rmf_data(arg, make_copy=False):
# Flatten f_chan and n_chan vectors into 1D arrays as crates does
# according to group
#
# TODO: should the n_grp > 0 filter be applied before all this?
if data['f_chan'].ndim > 1 and data['n_chan'].ndim > 1:
f_chan = []
n_chan = []
Expand All @@ -806,6 +798,8 @@ def get_rmf_data(arg, make_copy=False):
f_chan.append(fch[i])
n_chan.append(nch[i])

# This automatically filters out rows where N_GRP is 0.
#
data['f_chan'] = numpy.asarray(f_chan, SherpaUInt)
data['n_chan'] = numpy.asarray(n_chan, SherpaUInt)
else:
Expand Down

0 comments on commit 649d733

Please sign in to comment.