Add flux spectra

clarsach/spectrum.py

@@ -12,8 +12,56 @@
 
 __all__ = ['XSpectrum']
 
+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
 
@@ -94,3 +142,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)

clarsach/tests/test_spectrum.py

@@ -7,6 +7,7 @@
                          ('HETG',"data/fake_meg_m1.pha"),
                          ('HETG',"data/fake_meg_p1.pha"),
                          ('ACIS',"data/fake_acis.pha")])
+
 def test_load_xspectrum(ttype, filename):
     test = XSpectrum(filename, telescope=ttype)
     assert isinstance(test, XSpectrum)