diff --git a/pyxsis/io/__init__.py b/pyxsis/io/__init__.py index c985645..d800551 100644 --- a/pyxsis/io/__init__.py +++ b/pyxsis/io/__init__.py @@ -1 +1,2 @@ from .chandra_hetg import load_chandra_hetg +from .xmm_rgs import load_xmm_rgs diff --git a/pyxsis/io/xmm_rgs.py b/pyxsis/io/xmm_rgs.py new file mode 100644 index 0000000..9650f2b --- /dev/null +++ b/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) diff --git a/pyxsis/plot.py b/pyxsis/plot.py index 6034f5f..362ac69 100644 --- a/pyxsis/plot.py +++ b/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): diff --git a/pyxsis/xrayspectrum1d.py b/pyxsis/xrayspectrum1d.py index 07b9475..cbf7ad4 100644 --- a/pyxsis/xrayspectrum1d.py +++ b/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,10 +32,14 @@ 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. @@ -43,22 +47,27 @@ class XraySpectrum1D(Spectrum1D): 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,8 +476,8 @@ 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 """ @@ -472,7 +485,7 @@ def __init__(self, e_low, e_high, eff_area, 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