Able to read pha2 files and deal with ARF and RMF not specified in file

.gitignore

@@ -7,4 +7,4 @@
 *.rmf
 *.grmf
 *.gz
-data/tgcat/*
+

clarsach/__init__.py

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

clarsach/respond.py

@@ -2,13 +2,51 @@
 
 import numpy as np
 import astropy.io.fits as fits
+from astropy.units import si
 
 __all__ = ["RMF", "ARF"]
 
+CONST_HC = 12.398418573430595   # Copied from ISIS, [keV angs]
+KEV      = ['kev','keV']
+ANGS     = ['angs','angstrom','Angstrom','angstroms','Angstroms']
+ALLOWED_UNITS = KEV + ANGS
+
+def _Angs_keV(q):
+    """
+    Convert between keV and angs using hc converted to units of keV Angs
+
+    Parameters
+    ----------
+    q : numpy.ndarry
+        An array containing the array of interest (either keV or Angs)
+
+    Returns
+    -------
+    hc/q : numpy.ndarray
+        The converted values, monotonically increasing
+
+    sl   : slice
+        Slice used to ensure that the output is monotonically increasing
+        Is either [::1] or [::-1] (depending on order of input q)
+        This output is helpful for sorting data related to q
+    """
+    def _is_monotonically_increasing(x):
+        return all(x[1:] > x[:-1])
+
+    # Sometimes angs bins listed in reverse angstrom values (to match energies),
+    # in which case, no need to reverse
+    sl = slice(None, None, 1)
+    # Need to use reverse values if the bins are listed in increasing order
+    if _is_monotonically_increasing(q):
+        sl = slice(None, None, -1)
+
+    result = CONST_HC/q[sl]
+    return result, sl
+
+
 class RMF(object):
 
     def __init__(self, filename):
-
         self._load_rmf(filename)
         pass
 
@@ -256,7 +294,6 @@ def apply_rmf(self, spec):
 class ARF(object):
 
     def __init__(self, filename):
-
         self._load_arf(filename)
         pass
 

clarsach/spectrum.py

@@ -1,34 +1,36 @@
 import numpy as np
 import os
 
-from clarsach.respond import RMF, ARF
+from clarsach.respond import RMF, ARF, _Angs_keV
 from astropy.io import fits
+from astropy.units import si
 
 __all__ = ['XSpectrum']
 
-ALLOWED_UNITS      = ['keV','angs','angstrom','kev']
-ALLOWED_TELESCOPES = ['HETG','ACIS']
+KEV      = ['kev','keV']
+ANGS     = ['angs','angstrom','Angstrom','angstroms','Angstroms']
 
-CONST_HC    = 12.398418573430595   # Copied from ISIS, [keV angs]
-UNIT_LABELS = dict(zip(ALLOWED_UNITS, ['Energy (keV)', 'Wavelength (angs)']))
+ALLOWED_UNITS = KEV + ANGS
+ALLOWED_TELESCOPES = ['HETG','ACIS']
 
 # Not a very smart reader, but it works for HETG
 class XSpectrum(object):
-    def __init__(self, filename, telescope='HETG'):
+    def __init__(self, filename, telescope='HETG', row=None, arf=None, rmf=None, verbose=True):
         assert telescope in ALLOWED_TELESCOPES
 
         self.__store_path(filename)
 
         if telescope == 'HETG':
-            self._read_chandra(filename)
+            self._read_chandra(filename, arf=arf, rmf=rmf, row=row)
         elif telescope == 'ACIS':
-            self._read_chandra(filename)
+            self._read_chandra(filename, arf=arf, rmf=rmf)
 
-        if self.bin_unit != self.arf.e_unit:
-            print("Warning: ARF units and pha file units are not the same!!!")
+        if verbose:
+            if (self.arf is not None) and (self.bin_unit != self.arf.e_unit):
+                print("Warning: ARF units and pha file units are not the same!!!")
 
-        if self.bin_unit != self.rmf.energ_unit:
-            print("Warning: RMF units and pha file units are not the same!!!")
+            if (self.rmf is not None) and (self.bin_unit != self.rmf.energ_unit):
+                print("Warning: RMF units and pha file units are not the same!!!")
 
         return
 
@@ -62,51 +64,113 @@ def apply_resp(self, mflux, exposure=None):
         -------
         count_model : numpy.ndarray
             The model spectrum in units of counts/bin
+
+        If no ARF file exists, it will return the model flux after applying the RMF
+        If no RMF file exists, it will return the model flux after applying the ARF (with a warning)
+        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)
         else:
             mrate = mflux
 
-        count_model = self.rmf.apply_rmf(mrate)
-
-        return count_model
+        if self.rmf is not None:
+            count_model = self.rmf.apply_rmf(mrate)
+            return count_model
+        else:
+            print("Caution: no response file specified")
+            return mrate
 
     @property
     def bin_mid(self):
         return 0.5 * (self.bin_lo + self.bin_hi)
 
-    @property
-    def is_monotonically_increasing(self):
-        return all(self.bin_lo[1:] > self.bin_lo[:-1])
+    def _read_chandra(self, filename, arf=None, rmf=None, row=None):
+        TG_PART = {1:'HEG', 2:'MEG'}
+        this_dir = os.path.dirname(os.path.abspath(filename))
+        ff   = fits.open(filename)
+        data = ff[1].data
+
+        if row is not None:
+            assert row > 0
+            self.bin_lo = data['BIN_LO'][row-1]
+            self.bin_hi = data['BIN_HI'][row-1]
+            self.counts = data['COUNTS'][row-1]
+            tgp, tgm    = data['TG_PART'][row-1], data['TG_M'][row-1]
+            self.name   = "%s m=%d" % (TG_PART[tgp], tgm)
+        else:
+            self.bin_lo = data['BIN_LO']
+            self.bin_hi = data['BIN_HI']
+            self.counts = data['COUNTS']
+
+        # Deal with ARF and RMF
+        # Allow user to override file choice at the start, otherwise read filenames from header
+        # By default, the arf and rmf will be set to None (see kwargs in __init__ function call)
+        # If the filename specified is 'none', the ARF or RMF will be set to None
+        # If the filename is not 'none', the appropriate file will be loaded automatically
+        if rmf is not None:
+            self.rmf_file = rmf
+        else:
+            try:
+                fname = ff[1].header['RESPFILE']
+                if fname == 'none': self.rmf_file = None
+                else: self.rmf_file = this_dir + "/" + fname
+            except:
+                self.rmf_file = None
+
+        if arf is not None:
+            self.arf_file = arf
+        else:
+            try:
+                fname = ff[1].header['ANCRFILE']
+                if fname == 'none': self.arf_file = None
+                else: self.arf_file = this_dir + "/" + fname
+            except:
+                self.arf_file = None
+
+        self._assign_arf(self.arf_file)
+        self._assign_rmf(self.rmf_file)
+        self.bin_unit = data.columns['BIN_LO'].unit
+        self.exposure = ff[1].header['EXPOSURE']  # seconds
+        ff.close()
+
+    def _assign_arf(self, arf_inp):
+        if isinstance(arf_inp, str):
+            self.arf = ARF(arf_inp)
+        else:
+            self.arf = arf_inp
 
-    def _change_units(self, unit):
+    def _assign_rmf(self, rmf_inp):
+        if isinstance(rmf_inp, str):
+            self.rmf = RMF(rmf_inp)
+        else:
+            self.rmf = rmf_inp
+
+    def _return_in_units(self, unit):
         assert unit in ALLOWED_UNITS
         if unit == self.bin_unit:
             return (self.bin_lo, self.bin_hi, self.bin_mid, self.counts)
         else:
             # Need to use reverse values if the bins are listed in increasing order
-            if self.is_monotonically_increasing:
-                sl  = slice(None, None, -1)
-                print("is monotonically increasing")
-            # Sometimes its listed in reverse angstrom values (to match energies),
-            # in which case, no need to reverse
-            else:
-                sl  = slice(None, None, 1)
-                print("is NOT monotonically increasing")
-            new_lo  = CONST_HC/self.bin_hi[sl]
-            new_hi  = CONST_HC/self.bin_lo[sl]
+            new_lo, sl = _Angs_keV(self.bin_hi)
+            new_hi, sl = _Angs_keV(self.bin_lo)
             new_mid = 0.5 * (new_lo + new_hi)
             new_cts = self.counts[sl]
             return (new_lo, new_hi, new_mid, new_cts)
 
-    def hard_set_units(self, unit):
-        new_lo, new_hi, new_mid, new_cts = self._change_units(unit)
-        self.bin_lo = new_lo
-        self.bin_hi = new_hi
+    def _setbins_to_keV(self):
+        assert self.bin_unit in ANGS
+        new_bhi, sl = _Angs_keV(self.bin_lo)
+        new_blo, sl = _Angs_keV(self.bin_hi)
+        new_cts  = self.counts[sl]
+
+        # Now hard set everything
+        self.bin_lo = new_blo
+        self.bin_hi = new_bhi
         self.counts = new_cts
-        self.bin_unit = unit
+        self.bin_unit = si.keV
+        return
 
         return
 
@@ -118,30 +182,3 @@ def plot(self, ax, xunit='keV', **kwargs):
         ax.step(lo, cts, where='post', **kwargs)
         ax.set_xlabel(UNIT_LABELS[xunit])
         ax.set_ylabel('Counts')
-
-        return ax
-
-    def _read_chandra(self, filename):
-        this_dir = os.path.dirname(os.path.abspath(filename))
-        ff   = fits.open(filename)
-        data = ff[1].data
-        hdr = ff[1].header
-
-        self.bin_lo   = data['BIN_LO']
-        self.bin_hi   = data['BIN_HI']
-        self.bin_unit = data.columns['BIN_LO'].unit
-        self.counts   = data['COUNTS']
-
-        self.rmf_file = this_dir + "/" + hdr['RESPFILE']
-        self.arf_file = this_dir + "/" + hdr['ANCRFILE']
-        self.rmf = RMF(self.rmf_file)
-        self.arf = ARF(self.arf_file)
-
-        if "EXPOSURE" in list(hdr.keys()):
-            self.exposure = hdr['EXPOSURE']  # seconds
-        else:
-            self.exposure = 1.0
-
-        ff.close()
-
-        return