Add flux spectra

clarsach/__init__.py

@@ -1,5 +1,5 @@
 #!/usr/bin/env python
 '''Clarsach'''
-from . import spectrum
+from .spectrum import XSpectrum, calculate_flux_spectrum
 from .respond import RMF, ARF
-from .models import *
+from .models import Powerlaw

clarsach/respond.py

@@ -260,6 +260,7 @@ def _load_arf(self, filename):
         self.e_high = np.array(data.field("ENERG_HI"))
         self.e_unit = data.columns["ENERG_LO"].unit
         self.specresp = np.array(data.field("SPECRESP"))
+        self.exposure = hdr["EXPOSURE"]
 
         if "FRACEXPO" in data.columns.names:
             self.fracexpo = data["FRACEXPO"]
@@ -268,7 +269,7 @@ def _load_arf(self, filename):
 
         return
 
-    def apply_arf(self, spec):
+    def apply_arf(self, spec, apply_exp=True, apply_frac_exp=True):
         """
         Fold the spectrum through the ARF.
 
@@ -281,6 +282,12 @@ def apply_arf(self, spec):
         spec : numpy.ndarray
             The (model) spectrum to be folded
 
+        apply_exp : bool, default True
+            If True, apply exposure to ARF-corrected spectrum
+
+        apply_frac_exp : bool, default True
+            If True, apply fractional exposure to ARF-corrected spectrum
+
         Returns
         -------
         s_arf : numpy.ndarray
@@ -292,4 +299,11 @@ def apply_arf(self, spec):
                                                       "be of same size as the " \
                                                       "ARF array."
 
-        return np.array(spec) * self.specresp
+        spec_corr = np.array(spec) * self.specresp
+
+        if apply_exp:
+            spec_corr *= self.exposure
+        if apply_frac_exp:
+            spec_corr *= self.fracexpo
+
+        return spec_corr

clarsach/spectrum.py

@@ -4,16 +4,64 @@
 from clarsach.respond import RMF, ARF
 from astropy.io import fits
 
-__all__ = ['XSpectrum']
+__all__ = ['XSpectrum', 'calculate_flux_spectrum']
 
 ALLOWED_UNITS      = ['keV','angs','angstrom','kev']
 ALLOWED_TELESCOPES = ['HETG','ACIS']
 
 CONST_HC    = 12.398418573430595   # Copied from ISIS, [keV angs]
 UNIT_LABELS = dict(zip(ALLOWED_UNITS, ['Energy (keV)', 'Wavelength (angs)']))
 
+def calculate_flux_spectrum(spec, counts=None):
+    """
+    Convert a spectrum from units of photon counts to
+    units of photon flux.
+
+    **Warning**: You should do this for plotting *only*.
+    Converting from flux to counts, as is done in the detector,
+    involves a convolution with the response matrix of the detector.
+    Deconvolving this correctly is very hard, and not done in this
+    function. This is for visualization purposes only!
+
+    Parameters
+    ----------
+    spec : clarsach.XSpectrum object
+        An XSpectrum object containing the X-ray spectrum
+
+    counts : iterable
+        An array with counts to be converted to the flux
+        If this is None, use the `counts` attribute of the `spec`
+        object.
+
+    Returns
+    -------
+    flux_spectrum:
+        The flux spectrum
+
+    """
+
+    # make a flat spectrum so that I can integrate
+    # ARF and RMF only
+    flat_model = np.ones_like(spec.counts) * (spec.rmf.energ_hi -
+                                              spec.rmf.energ_lo)
+
+    # now apply ARF to the flat model
+    m_arf = spec.arf.apply_arf(flat_model)
+
+    # apply RMF to the flat model
+    m_rmf = spec.rmf.apply_rmf(m_arf)
+
+    # divide the observed counts by the flat model with the ARF/RMF applied
+    if counts is None:
+        flux_spec = spec.counts / m_rmf / spec.exposure
+    else:
+        flux_spec = counts / m_rmf / spec.exposure
+
+    return flux_spec
+
 # Not a very smart reader, but it works for HETG
 class XSpectrum(object):
+
     def __init__(self, filename, telescope='HETG'):
         assert telescope in ALLOWED_TELESCOPES
 
@@ -35,11 +83,30 @@ def __store_path(self, filename):
         return
 
     def apply_resp(self, mflux):
-        # Given a model flux spectrum, apply the response
-        mrate  = self.arf.apply_arf(mflux)  # phot/s per bin
-        mrate  *= self.exposure * self.arf.fracexpo  # phot per bin
-        result = self.rmf.apply_rmf(mrate)  # counts per bin
-        return result
+        """
+        Given a model flux spectrum, apply the response, both ARF and RMF.
+
+        Parameters
+        ----------
+        mflux : iterable
+            list or array with the model flux spectrum
+
+        Returns
+        -------
+        m_rmf: numpy.ndarray
+            The counts spectrum with responses applied
+
+        """
+        # For some instruments, the ARF could not exist
+        if self.arf is not None:
+            mrate  = self.arf.apply_arf(mflux, apply_exp=True,
+                                        apply_frac_exp=True)  # phot/s per bin
+        else:
+            mrate = mflux
+
+        m_rmf = self.rmf.apply_rmf(mrate)  # counts per bin
+
+        return m_rmf
 
     @property
     def bin_mid(self):
@@ -76,14 +143,24 @@ def hard_set_units(self, unit):
         self.counts = new_cts
         self.bin_unit = unit
 
-    def plot(self, ax, xunit='keV', **kwargs):
+    def plot(self, ax, xunit='keV', yunit="counts", **kwargs):
         lo, hi, mid, cts = self._change_units(xunit)
-        counts_err       = np.sqrt(cts)
-        ax.errorbar(mid, cts, yerr=counts_err,
-                    ls='', marker=None, color='k', capsize=0, alpha=0.5)
-        ax.step(lo, cts, where='post', **kwargs)
+        counts_err = np.sqrt(cts)
+        if yunit == "counts":
+            ax.errorbar(mid, cts, yerr=counts_err,
+                        ls='', marker=None, color='k', capsize=0, alpha=0.5)
+            ax.set_ylabel('Counts')
+            ax.step(lo, cts, where='post', **kwargs)
+
+        elif yunit == "flux":
+            flux = calculate_flux_spectrum(self)
+            flux_err = calculate_flux_spectrum(self, counts_err)
+            ax.errorbar(mid, flux, yerr=flux_err,
+                        ls='', marker=None, color='k', capsize=0, alpha=0.5)
+            ax.set_ylabel('Flux')
+            ax.step(lo, flux, where='post', **kwargs)
+
         ax.set_xlabel(UNIT_LABELS[xunit])
-        ax.set_ylabel('Counts')
 
     def _read_chandra(self, filename):
         this_dir = os.path.dirname(os.path.abspath(filename))
@@ -99,3 +176,5 @@ def _read_chandra(self, filename):
         self.arf = ARF(self.arf_file)
         self.exposure = ff[1].header['EXPOSURE']  # seconds
         ff.close()
+
+        self.flux = calculate_flux_spectrum(self, None)