Modify pyxsis to handle XMM RGS data

pyxsis/io/__init__.py

@@ -1 +1,2 @@
 from .chandra_hetg import load_chandra_hetg
+from .xmm_rgs import load_xmm_rgs

pyxsis/io/xmm_rgs.py

@@ -0,0 +1,53 @@
+import os
+import numpy as np
+
+from astropy.io import fits
+from astropy.units import Unit
+
+from .. import XBinSpectrum
+
+__all__ = ['load_xmm_rgs']
+
+def load_xmm_rgs(filename, arf=None, rmf=None):
+    """
+    Load XMM RGS spectral data from a file into a spectrum object.
+
+    **Inputs**
+
+    filename : str
+        The path to the FITS file
+
+    arf : str -or- pyxsis.xrayspectrum1d.ARF
+        Filename for the area response file (ARF) or a pre-loaded AreaResponse object
+        For RGS data, there seems to be no separate ARF file. So this will usually be None
+
+    rmf : str -or- pyxsis.xrayspectrum1d.RMF
+        Filename for the response matrix file (RMF) or a pre-loaded ResponseMatrix object
+
+    **Returns**
+
+    pyxsis XraySpectrum1D object representing the data in the input FITS file
+    """
+    this_dir = os.path.dirname(os.path.abspath(filename))
+
+    with fits.open(filename) as hdu:
+        header  = hdu[0].header
+        meta    = {'header': header}
+        data    = hdu[1].data
+        datahdr = hdu[1].header
+
+        bin_unit  = Unit(datahdr['TCUNI1'])
+        bin_cen   = datahdr['TCRVL1'] # center of the reference bin
+        bin_cen_i = datahdr['TLMIN1']-1 # reference bin (assuming 1 refers to first bin)
+        bin_width = datahdr['TCDLT1'] # width of each channel
+        bin_lo    = (np.arange(datahdr['TLMAX1']) * bin_width + bin_cen - bin_width / 2.0) * bin_unit
+        bin_hi    = (np.arange(datahdr['TLMAX1']) * bin_width + bin_cen + bin_width / 2.0) * bin_unit
+
+        counts   = data['COUNTS'] * Unit('count')
+        exposure = datahdr['EXPOSURE'] * Unit('second')
+
+        if rmf is None:
+            rmf = this_dir + "/" + hdu[1].header['RESPFILE']
+
+    return XBinSpectrum(bin_lo, bin_hi, counts, areascal=data['AREASCAL'],
+                          exposure=exposure, arf=arf, rmf=rmf)

pyxsis/plot.py

@@ -135,7 +135,8 @@ def plot_unfold(ax, spectrum, xunit='keV', perbin=False,
     """
     # Models will always be in keV bin units
     # a non-model of ones (integrated)
-    no_mod  = np.ones_like(spectrum.arf.eff_area)
+    # This set of calcs will get the effective response for each bin (cm^2 ct / phot)
+    no_mod  = np.ones(len(spectrum.counts))
     eff_tmp = spectrum.apply_response(no_mod)
 
     def _bin_exp(exp, binning):

pyxsis/xrayspectrum1d.py

@@ -17,7 +17,7 @@ class XraySpectrum1D(Spectrum1D):
     """Spectrum properties specific to X-ray data
 
     **Attributes**
-    
+
     Inherits from specutils Spectrum1D object. The following inputs
     are stored as additional attributes.
 
@@ -32,33 +32,42 @@ class XraySpectrum1D(Spectrum1D):
 
     exposure : astropy.Quantity
         Exposure time for the dataset
-
+
+    areascal : numpy.ndarray
+        Scale factor for effective region size (or effective area).
+        A data column found in XMM spectra.
+
     arf : ARF object
         Telescope response file describing the effective area as a function of photon energy.
-    
+
     rmf : RMF object
         Telescope response file, a 2D matrix, describing the detector signal (pulse heights) distribution as a function photon energy.
 
     rest_value : astropy.units.Quantity, default 0 Angstrom
         See Spectrum1D rest_value input
     """
     def __init__(self, bin_lo, bin_hi, counts, exposure,
-                 arf=None, rmf=None, **kwargs):
+                 areascal=1.0, arf=None, rmf=None, **kwargs):
         """
         **Inputs**
 
         bin_lo : astropy.Quantity
             The left edges for bin values
-        
+
         bin_hi : astropy.Quantity
             The right edges for bin values
-        
+
         counts : astropy.Quantity
             Counts histogram for the X-ray spectrum
-        
+
         exposure : astropy.Quantity
             Exposure time for the dataset
 
+        areascal : numpy.ndarray
+            Scale factor for effective region size (or effective area). A data
+            column found in XMM spectra. This is used to adjust the output of
+            apply_response
+
         arf : ARF or string (default: None)
             Strings will be passed to ARF.__init__
             All other input types are stored as arf attribute
@@ -75,6 +84,7 @@ def __init__(self, bin_lo, bin_hi, counts, exposure,
         self.bin_lo = bin_lo
         self.bin_hi = bin_hi
         self.exposure = exposure
+        self.areascal = areascal
         self.assign_rmf(rmf)
         self.assign_arf(arf)
 
@@ -88,7 +98,7 @@ def assign_arf(self, arf_inp, **kwargs):
         Assign an auxiliary response file (ARF) object to the XraySpectrum1D object
 
         **Inputs**
-        
+
         arf_inp : string
             File name for the area response file (FITS file)
 
@@ -107,12 +117,12 @@ def assign_rmf(self, rmf_inp):
         Assign a redistribution matrix file (RMF) object to the XraySpectrum1D object
 
         **Input**
-        
+
         rmf_inp : string
             File name for the response matrix file (FITS file)
 
         **Returns**
-        
+
         Modifies the XraySpectrum1D.rmf attribute
         """
         if isinstance(rmf_inp, str):
@@ -132,7 +142,7 @@ def apply_response(self, mflux, exposure=None):
         bins as in the ARF (where the ARF exists)!
 
         **Inputs**
-        
+
         mflux : astropy.units.Quantity
             A list or array with the model flux values,
             typically with units of phot/s/cm^-2
@@ -146,7 +156,7 @@ def apply_response(self, mflux, exposure=None):
             default behaviour and manually set the exposure time to use.
 
         **Returns**
-        
+
         model_counts : numpy.ndarray
             The model spectrum in units of counts/bin
 
@@ -155,16 +165,19 @@ def apply_response(self, mflux, exposure=None):
         If no ARF and no RMF, it will return the model flux spectrum (with a warning)
         """
         if self.arf is not None:
-            mrate  = self.arf.apply_arf(mflux, exposure=exposure) # ct/bin with no RMF applied
+            mrate  = self.arf.apply_arf(mflux, exposure=exposure)
         else:
-            mrate = mflux # assumes ct/bin is the input with no RMF applied
+            print("Warning: No ARF assigned. Assuming detector area is")
+            print("contained in response matrix with units of s * cm^2 ")
+            mrate = mflux * u.cm**2 * u.second
 
         if self.rmf is not None:
             result = self.rmf.apply_rmf(mrate.value)
         else:
+            print("Warning: No RMF assigned.")
             result = mrate
 
-        return result
+        return result * self.areascal
 
 ## ----  Supporting response file objects
 
@@ -176,7 +189,7 @@ def __init__(self, energ_lo=None, energ_hi=None, matrix=None, energ_unit=None,
         Redistribution Matrix File (RMF) for an X-ray spectrum.
 
         **Attributes**
-        
+
         filename : str
             Stores the filename used with RMF.read
 
@@ -276,12 +289,12 @@ def __get_tlmin(self, h):
         Get the tlmin keyword for `F_CHAN`.
 
         **Inputs**
-        
+
         h : an astropy.io.fits.hdu.table.BinTableHDU object
             The extension containing the `F_CHAN` column
 
         **Returns**
-        
+
         tlmin : int
             The tlmin keyword
         """
@@ -372,12 +385,12 @@ def apply_rmf(self, model):
         sure to get the indices right.
 
         **Inputs**
-        
+
         model : numpy.ndarray
             The (model) spectrum to be folded
 
         **Returns**
-        
+
         counts : numpy.ndarray
             The (model) spectrum after folding, in
             counts/s/channel
@@ -443,7 +456,7 @@ def __init__(self, e_low, e_high, eff_area,
         Auxiliary Response File (ARF) for an X-ray spectrum.
 
         **Attributes**
-        
+
         filename : str
             Stores the filename used with ARF.read
 
@@ -463,16 +476,16 @@ def __init__(self, e_low, e_high, eff_area,
         fracexpo : float or numpy.ndarray
             Fractional exposure time for the spectrum
             (sometimes constant, sometimes dependent on spectral channel).
-            These values are stored for reference; generally, they are already
-            accounted for in the eff_area array.
+            These values are stored for reference; they are already
+            accounted for in the eff_area array for the HETG data, it seems.
 
         e_mid : Property that returns the middle of each energy bin
         """
         self.filename = None
         self.e_low    = e_low
         self.e_high   = e_high
         self.eff_area = eff_area
-        self.fracexpo = fracexpo
+        self.fracexpo = fracexpo # a feature of Chandra HETG ARFS
         self.exposure = exposure
 
     @property
@@ -485,7 +498,7 @@ def read(filename, block='SPECRESP'):
         Load an ARF object from FITS file.
 
         **Inputs**
-        
+
         filename : str
             Path to the FITS file
 
@@ -494,7 +507,7 @@ def read(filename, block='SPECRESP'):
             under an extension other than "SPECRESP"
 
         **Returns**
-        
+
         The ARF object that is represented by the FITS file
         """
         # open the FITS file and extract the SPECRESP block
@@ -539,7 +552,7 @@ def apply_arf(self, model, exposure=None):
         multiplication with the input spectrum.
 
         **Parameters**
-        
+
         model : numpy.ndarray
             The model spectrum to which the arf will be applied
 
@@ -552,7 +565,7 @@ def apply_arf(self, model, exposure=None):
             value.
 
         **Returns**
-        
+
         s_arf : numpy.ndarray
             The (model) spectrum after folding, in
             counts/s/channel