diff --git a/.gitignore b/.gitignore index 9a5a420..aa822a3 100644 --- a/.gitignore +++ b/.gitignore @@ -7,4 +7,4 @@ *.rmf *.grmf *.gz -data/tgcat/* \ No newline at end of file + diff --git a/clarsach/__init__.py b/clarsach/__init__.py index bbcfbde..65e605e 100644 --- a/clarsach/__init__.py +++ b/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 * diff --git a/clarsach/respond.py b/clarsach/respond.py index 3ebac4f..7d5b936 100644 --- a/clarsach/respond.py +++ b/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 diff --git a/clarsach/spectrum.py b/clarsach/spectrum.py index fd2d9e4..41a0317 100644 --- a/clarsach/spectrum.py +++ b/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,6 +64,10 @@ 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: @@ -69,44 +75,102 @@ def apply_resp(self, mflux, exposure=None): 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 \ No newline at end of file diff --git a/data/arfs/aciss_heg1_cy19.garf b/data/arfs/aciss_heg1_cy19.garf deleted file mode 100644 index 01512ca..0000000 --- a/data/arfs/aciss_heg1_cy19.garf +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:e974a45ecbbd439c9196ec13e3b01cb8ea17e830b7147163bbb996f6a39704b7 -size 256320 diff --git a/data/arfs/aciss_hetg0_cy19.arf b/data/arfs/aciss_hetg0_cy19.arf deleted file mode 100644 index c6cd159..0000000 --- a/data/arfs/aciss_hetg0_cy19.arf +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:f72c57747babc12dd8206ec6d949ec34617bf538ea99223793796a755d3adb02 -size 40320 diff --git a/data/fake_acis.pha b/data/fake_acis.pha deleted file mode 100644 index 01735de..0000000 --- a/data/fake_acis.pha +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:55ddb584833ceb37064453f530949aaa8eac9f9c9d6b0a32a61d7b865fc17f69 -size 43200 diff --git a/data/fake_heg_m1.pha b/data/fake_heg_m1.pha deleted file mode 100644 index 10a13f6..0000000 --- a/data/fake_heg_m1.pha +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:1c3bacc9ae8038352ad8c4272baf444a910d974553edffc93441503063892155 -size 239040 diff --git a/data/fake_heg_p1.pha b/data/fake_heg_p1.pha deleted file mode 100644 index ae2fbc0..0000000 --- a/data/fake_heg_p1.pha +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:d2cbe0b0b506674fa4563944639b5ba44f22d89df0f020b52b5a092151b68686 -size 239040 diff --git a/data/fake_meg_m1.pha b/data/fake_meg_m1.pha deleted file mode 100644 index 3d3860c..0000000 --- a/data/fake_meg_m1.pha +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:4e14cc41a5efb622ba0c9843f95a695b1d1631060a89494e0dc934bea1e84ce1 -size 273600 diff --git a/data/fake_meg_p1.pha b/data/fake_meg_p1.pha deleted file mode 100644 index 55083af..0000000 --- a/data/fake_meg_p1.pha +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:3dadc5be214311c127222e4bd797d50889156f8a79be7a946f9dd47fd4ea1b5a -size 273600 diff --git a/data/rmfs/aciss_heg-1_cy19.grmf b/data/rmfs/aciss_heg-1_cy19.grmf deleted file mode 100644 index 6d5038c..0000000 --- a/data/rmfs/aciss_heg-1_cy19.grmf +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:e751bb311b07b694f43d305a187a22e67284d46d715ef2a707a80a760c3edb46 -size 7151040 diff --git a/data/rmfs/aciss_heg1_cy19.grmf b/data/rmfs/aciss_heg1_cy19.grmf deleted file mode 100644 index 1948722..0000000 --- a/data/rmfs/aciss_heg1_cy19.grmf +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:1190c43cda91d89285c400fa44ef0641bca2ab80307ef5647d9d639f597e2fd0 -size 7151040 diff --git a/data/rmfs/aciss_hetg0_cy19.rmf b/data/rmfs/aciss_hetg0_cy19.rmf deleted file mode 100644 index fdd3680..0000000 --- a/data/rmfs/aciss_hetg0_cy19.rmf +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:66005031e16e255ca4ec666d92fada40076d9dac1d109fa09d24c15eb0718467 -size 1624320