diff --git a/docs/api.rst b/docs/api.rst index cebc3df21..a2b32fced 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -263,6 +263,16 @@ I/O Functionality .. automodule:: stingray.io :members: +Mission-specific I/O +-------------------- + +.. automodule:: stingray.mission_support.missions + :members: + +.. automodule:: stingray.mission_support.rxte + :members: + + Other Utility Functions ----------------------- diff --git a/docs/changes/816.feature.rst b/docs/changes/816.feature.rst new file mode 100644 index 000000000..f40b9f688 --- /dev/null +++ b/docs/changes/816.feature.rst @@ -0,0 +1 @@ +Add support for XTE data diff --git a/stingray/base.py b/stingray/base.py index 49c4dc4fb..e8727c5aa 100644 --- a/stingray/base.py +++ b/stingray/base.py @@ -2660,6 +2660,7 @@ def analyze_segments(self, func, segment_size, fraction_step=1, **kwargs): warnings.warn( f"Segment {i} ({tst}--{tsp}) has one data point or less. Skipping it " ) + continue lc_filt = self[st:sp] lc_filt.gti = np.asanyarray([[tst, tsp]]) diff --git a/stingray/io.py b/stingray/io.py index 8e3f48419..e3bb10de9 100644 --- a/stingray/io.py +++ b/stingray/io.py @@ -16,9 +16,22 @@ import stingray.utils as utils from stingray.loggingconfig import setup_logger -from .utils import assign_value_if_none, is_string, order_list_of_arrays, is_sorted +from .utils import ( + assign_value_if_none, + is_string, + order_list_of_arrays, + is_sorted, + make_dictionary_lowercase, +) from .gti import get_gti_from_all_extensions, load_gtis +from .mission_support import ( + read_mission_info, + rough_calibration, + get_rough_conversion_function, + mission_specific_event_interpretation, +) + # Python 3 import pickle @@ -93,50 +106,6 @@ def pi_to_energy(pis, rmf_file): return es -def rough_calibration(pis, mission): - """Make a rough conversion between PI channel and energy. - - Only works for NICER, NuSTAR, IXPE, and XMM. - - Parameters - ---------- - pis: float or array of floats - PI channels in data - mission: str - Mission name - - Returns - ------- - energies : float or array of floats - Energy values - - Examples - -------- - >>> rough_calibration(0, 'nustar') - 1.6 - >>> rough_calibration(0.0, 'ixpe') - 0.0 - >>> # It's case-insensitive - >>> rough_calibration(1200, 'XMm') - 1.2 - >>> rough_calibration(10, 'asDf') - Traceback (most recent call last): - ... - ValueError: Mission asdf not recognized - >>> rough_calibration(100, 'nicer') - 1.0 - """ - if mission.lower() == "nustar": - return pis * 0.04 + 1.6 - elif mission.lower() == "xmm": - return pis * 0.001 - elif mission.lower() == "nicer": - return pis * 0.01 - elif mission.lower() == "ixpe": - return pis / 375 * 15 - raise ValueError(f"Mission {mission.lower()} not recognized") - - def get_file_extension(fname): """Get the extension from the file name. @@ -198,69 +167,6 @@ def high_precision_keyword_read(hdr, keyword): return None -def _patch_mission_info(info, mission=None): - """Add some information that is surely missing in xselect.mdb. - - Examples - -------- - >>> info = {'gti': 'STDGTI'} - >>> new_info = _patch_mission_info(info, mission=None) - >>> assert new_info['gti'] == info['gti'] - >>> new_info = _patch_mission_info(info, mission="xmm") - >>> new_info['gti'] - 'STDGTI,GTI0' - """ - if mission is None: - return info - if mission.lower() == "xmm" and "gti" in info: - info["gti"] += ",GTI0" - return info - - -def read_mission_info(mission=None): - """Search the relevant information about a mission in xselect.mdb.""" - curdir = os.path.abspath(os.path.dirname(__file__)) - fname = os.path.join(curdir, "datasets", "xselect.mdb") - - # If HEADAS is defined, search for the most up-to-date version of the - # mission database - if os.getenv("HEADAS"): - hea_fname = os.path.join(os.getenv("HEADAS"), "bin", "xselect.mdb") - if os.path.exists(hea_fname): - fname = hea_fname - if mission is not None: - mission = mission.lower() - - db = {} - with open(fname) as fobj: - for line in fobj.readlines(): - line = line.strip() - if mission is not None and not line.lower().startswith(mission): - continue - if line.startswith("!") or line == "": - continue - allvals = line.split() - string = allvals[0] - value = allvals[1:] - if len(value) == 1: - value = value[0] - - data = string.split(":")[:] - if mission is None: - if data[0] not in db: - db[data[0]] = {} - previous_db_step = db[data[0]] - else: - previous_db_step = db - data = data[1:] - for key in data[:-1]: - if key not in previous_db_step: - previous_db_step[key] = {} - previous_db_step = previous_db_step[key] - previous_db_step[data[-1]] = value - return _patch_mission_info(db, mission) - - def _case_insensitive_search_in_list(string, list_of_strings): """Search for a string in a list of strings, in a case-insensitive way. @@ -357,16 +263,22 @@ def get_key_from_mission_info(info, key, default, inst=None, mode=None): >>> get_key_from_mission_info(info, "ghghg", "BU", inst="C", mode="M1") 'BU' """ - filt_info = copy.deepcopy(info) - if inst is not None and inst in filt_info: - filt_info.update(info[inst]) - filt_info.pop(inst) - if mode is not None and mode in filt_info: - filt_info.update(info[inst][mode]) - filt_info.pop(mode) + filt_info = make_dictionary_lowercase(info, recursive=True) + key = key.lower() + if inst is not None: + inst = inst.lower() + if inst in filt_info: + filt_info.update(filt_info[inst]) + filt_info.pop(inst) + if mode is not None: + mode = mode.lower() + if mode in filt_info: + filt_info.update(filt_info[mode]) + filt_info.pop(mode) if key in filt_info: return filt_info[key] + return default @@ -644,6 +556,10 @@ def load_events_and_gtis( mission_key = "TELESCOP" mission = probe_header[mission_key].lower() + mission_specific_processing = mission_specific_event_interpretation(mission) + + mission_specific_processing(hdulist) + db = read_mission_info(mission) instkey = get_key_from_mission_info(db, "instkey", "INSTRUME") instr = mode = None @@ -760,11 +676,15 @@ def load_events_and_gtis( returns.gti_list = gti_list returns.pi_list = pi returns.cal_pi_list = cal_pi + if "energy" in additional_data and np.any(additional_data["energy"] > 0.0): returns.energy_list = additional_data["energy"] else: try: - returns.energy_list = rough_calibration(cal_pi, mission) + func = get_rough_conversion_function( + mission, instrument=instr, epoch=t_start / 86400 + mjdref + ) + returns.energy_list = func(cal_pi, detector_id=detector_id) logger.info( f"A default calibration was applied to the {mission} data. " "See io.rough_calibration for details. " diff --git a/stingray/mission_support/__init__.py b/stingray/mission_support/__init__.py new file mode 100644 index 000000000..7515b2dd1 --- /dev/null +++ b/stingray/mission_support/__init__.py @@ -0,0 +1 @@ +from .missions import * diff --git a/stingray/mission_support/missions.py b/stingray/mission_support/missions.py new file mode 100644 index 000000000..8b1bb98a8 --- /dev/null +++ b/stingray/mission_support/missions.py @@ -0,0 +1,194 @@ +"""This module contains functions to interpret data from different missions. + +The key functions are: + +- `read_mission_info`: Search the relevant information about a mission in xselect.mdb. +- `get_rough_conversion_function`: Get a rough PI-Energy conversion function for a mission. +- `mission_specific_event_interpretation`: Get the mission-specific FITS interpretation + function. This function will take a FITS :class:`astropy.io.fits.HDUList` object and + modify it in place to make the read into Stingray easier. +- `rough_calibration` (obsolete): Make a rough conversion between PI channel and energy. + +Whenever a given mission needs complicate processing, its functions can be made available +for specific missions in their own separate modules. For example, the RXTE mission has its +own module, ``rxte.py``, which contains functions to interpret RXTE data. +""" + +import os +import warnings +from .rxte import rxte_calibration_func, rxte_pca_event_file_interpretation + + +def rough_calibration(pis, mission): + """Make a rough conversion between PI channel and energy. + + Only works for NICER, NuSTAR, IXPE, and XMM. + + Parameters + ---------- + pis: float or array of floats + PI channels in data + mission: str + Mission name + + Returns + ------- + energies : float or array of floats + Energy values + + Examples + -------- + >>> rough_calibration(0, 'nustar') + 1.62 + >>> rough_calibration(0.0, 'ixpe') + 0.0 + >>> # It's case-insensitive + >>> rough_calibration(1200, 'XMm') + 1.2 + >>> rough_calibration(10, 'asDf') + Traceback (most recent call last): + ... + ValueError: Mission asdf not recognized + >>> rough_calibration(100, 'nicer') + 1.0 + """ + if mission.lower() == "nustar": + return pis * 0.04 + 1.62 + elif mission.lower() == "xmm": + return pis * 0.001 + elif mission.lower() == "nicer": + return pis * 0.01 + elif mission.lower() == "ixpe": + return pis / 375 * 15 + raise ValueError(f"Mission {mission.lower()} not recognized") + + +def _patch_mission_info(info, mission=None): + """Add some information that is surely missing in xselect.mdb. + + Examples + -------- + >>> info = {'gti': 'STDGTI', 'ecol': 'PHA'} + >>> new_info = _patch_mission_info(info, mission=None) + >>> assert new_info['gti'] == info['gti'] + >>> new_info = _patch_mission_info(info, mission="xmm") + >>> new_info['gti'] + 'STDGTI,GTI0' + >>> new_info = _patch_mission_info(info, mission="xte") + >>> new_info['ecol'] + 'PHA' + """ + if mission is None: + return info + if mission.lower() == "xmm" and "gti" in info: + info["gti"] += ",GTI0" + if mission.lower() == "xte" and "ecol" in info: + info["ecol"] = "PHA" + info["ccol"] = "PCUID" + return info + + +def read_mission_info(mission=None): + """Search the relevant information about a mission in xselect.mdb.""" + curdir = os.path.abspath(os.path.dirname(__file__)) + fname = os.path.join(curdir, "..", "datasets", "xselect.mdb") + + # If HEADAS is defined, search for the most up-to-date version of the + # mission database + if os.getenv("HEADAS"): + hea_fname = os.path.join(os.getenv("HEADAS"), "bin", "xselect.mdb") + if os.path.exists(hea_fname): + fname = hea_fname + if mission is not None: + mission = mission.lower() + + db = {} + with open(fname) as fobj: + for line in fobj.readlines(): + line = line.strip() + if mission is not None and not line.lower().startswith(mission): + continue + if line.startswith("!") or line == "": + continue + allvals = line.split() + string = allvals[0] + value = allvals[1:] + if len(value) == 1: + value = value[0] + + data = string.split(":")[:] + if mission is None: + if data[0] not in db: + db[data[0]] = {} + previous_db_step = db[data[0]] + else: + previous_db_step = db + data = data[1:] + for key in data[:-1]: + if key not in previous_db_step: + previous_db_step[key] = {} + previous_db_step = previous_db_step[key] + previous_db_step[data[-1]] = value + return _patch_mission_info(db, mission) + + +def _wrap_function_ignoring_kwargs(func): + def func_wrapper(pi, **kwargs): + return func(pi) + + return func_wrapper + + +SIMPLE_CONVERSION_FUNCTIONS = { + "nustar": lambda pi: pi * 0.04 + 1.62, + "xmm": lambda pi: pi * 0.001, + "nicer": lambda pi: pi * 0.01, + "ixpe": lambda pi: pi / 375 * 15, +} + + +def get_rough_conversion_function(mission, instrument=None, epoch=None): + """Get a rough PI-Energy conversion function for a mission. + + The function should accept a PI channel and return the corresponding energy. + Additional keyword arguments (e.g. epoch, detector) can be passed to the function. + + Parameters + ---------- + mission : str + Mission name + instrument : str + Instrument onboard the mission + epoch : float + Epoch of the observation in MJD (important for missions updating their calibration). + + Returns + ------- + function + Conversion function + """ + + if mission.lower() in SIMPLE_CONVERSION_FUNCTIONS: + return _wrap_function_ignoring_kwargs(SIMPLE_CONVERSION_FUNCTIONS[mission.lower()]) + + if mission.lower() == "xte": + func = rxte_calibration_func(instrument, epoch) + return func + + raise ValueError(f"Mission {mission.lower()} not recognized") + + +def mission_specific_event_interpretation(mission): + """Get the mission-specific FITS interpretation function. + + This function will read a FITS :class:`astropy.io.fits.HDUList` object and modify it + in place to make the read into Stingray easier. + """ + + if mission.lower() == "xte": + return rxte_pca_event_file_interpretation + + def _empty(x): + return x + + return _empty diff --git a/stingray/mission_support/rxte.py b/stingray/mission_support/rxte.py new file mode 100644 index 000000000..d0842ec36 --- /dev/null +++ b/stingray/mission_support/rxte.py @@ -0,0 +1,337 @@ +import re +import numpy as np +from scipy.interpolate import interp1d +from astropy.time import Time +from astropy.table import Table + +c_match = re.compile(r"C\[(.*)\]") + +_EDGE_TIMES = [ + "1995-12-30T00:00:00.0", # launch + "1996-03-21T18:33:00.0", + "1996-04-15T23:05:00.0", + "1999-03-22T17:37:00.0", + "2000-05-13T00:00:00.0", + "2012-01-05T00:00:00.0", # decommissioning +] + +_EDGE_EPOCHS = Time(_EDGE_TIMES, format="isot", scale="utc").mjd + + +def _split_chan_spec(chan, sep="~"): + """Split a channel specification into a tuple of integers. + + If the specification is, e.g., ``10-13``, it will be interpreted as + ``(10, 13)``. The same thing will be omitted if the starting digit(s) of the first int are + omitted, e.g. ``10-3`` will still be interpreted as ``(10, 13)``. + If there is only one integer, it will be duplicated. + + Examples + -------- + >>> _split_chan_spec("10~13") + (10, 13) + >>> _split_chan_spec("10-3", sep="-") + (10, 13) + >>> _split_chan_spec("10") + (10, 10) + + """ + if sep in chan: + cs = chan.split(sep) + else: + cs = (chan, chan) + if len(cs[1]) < len(cs[0]): + c0 = cs[0] + c1 = cs[0][: -len(cs[1])] + cs[1] + cs = (c0, c1) + return (int(cs[0]), int(cs[1])) + + +def _split_C_string(string, sep="~"): + """Interpret the C string in the TEVTB2 key and return a list of channel tuples.""" + channel_list = [] + if ":" in string: + cs_range = string.split(":") + return [(c, c) for c in range(int(cs_range[0]), int(cs_range[1]) + 1)] + + for chan in string.split(","): + cs = _split_chan_spec(chan, sep) + channel_list.append(cs) + return channel_list + + +def _decode_energy_channels(tevtb2): + """Understand the channel information TEVTB2 key. + + Parameters + ---------- + tevtb2 : str + The TEVTB2 key from the FITS header. + + Returns + ------- + chans: list of tuples + A list of tuples, each containing the start and stop channel of a group of channels. + + Examples + -------- + >>> tevtb2 = '(M[1]{1},C[43,44~45]{6})' + >>> chans = _decode_energy_channels(tevtb2) + >>> assert chans == [(43, 43), (44, 45)] + >>> _decode_energy_channels('(M[1]{1})') + Traceback (most recent call last): + ... + ValueError: No C line found in the TEVTB2 key. + """ + dll_fmt_string_split = re.split(",(?=[A-Z])", tevtb2) + for line in dll_fmt_string_split: + if not line.startswith("C"): + continue + line = c_match.match(line).group(1) + break + else: + raise ValueError("No C line found in the TEVTB2 key.") + + return _split_C_string(line) + + +def pca_calibration_func(epoch): + """Return the appropriate calibration function for RXTE for a given observing epoch. + + This function has signature ``func(pha, detector_id)`` and gives the energy corresponding + to the PHA channel for the given detector (array values allowed). + + Internally, this is done by pre-allocating some arrays with the energy values for each + PHA channel and detector group (1-4 and 0, due to a damage that PCU 0 incurred in 2000), + and then returning a function that looks up the energy for each channel. + + This does not require any interpolation, as the calibration is tabulated for each channel, + and it is pretty efficient given the very small number of channels supported by the PCA (255). + + Parameters + ---------- + epoch : float + The epoch of the observation in MJD. + + Returns + ------- + conversion_function : callable + A function that converts PHA channel to energy. This function accepts + two arguments: the PHA channel and the PCU number. + + Examples + -------- + >>> conversion_function = pca_calibration_func(50082) + >>> float(conversion_function(10, 0)) + 3.04 + >>> conversion_function = pca_calibration_func(55930) + >>> float(conversion_function(10, 0)) + 4.53 + >>> float(conversion_function(10, 3)) + 4.49 + >>> assert np.array_equal(conversion_function(10, [0, 3]), [4.53, 4.49]) + >>> assert np.array_equal(conversion_function([10, 11], [0, 3]), [4.53, 4.90]) + """ + caltable = Table.read( + """ + Abs STD2 E1 E2 E3 E4 E5_0 E5_1234 + 0-4 0 1.51 1.61 1.94 2.13 1.95 2.06 + 5 1 1.76 1.91 2.29 2.54 2.38 2.47 + 6 2 2.02 2.22 2.64 2.96 2.81 2.87 + 7 3 2.27 2.52 2.99 3.37 3.24 3.28 + 8 4 2.53 2.83 3.35 3.79 3.67 3.68 + 9 5 2.78 3.13 3.70 4.21 4.09 4.09 + 10 6 3.04 3.44 4.05 4.63 4.53 4.49 + 11 7 3.30 3.75 4.41 5.04 4.96 4.90 + 12 8 3.55 4.05 4.76 5.46 5.39 5.31 + 13 9 3.81 4.36 5.12 5.88 5.82 5.71 + 14 10 4.07 4.67 5.47 6.30 6.25 6.12 + 15 11 4.32 4.98 5.82 6.72 6.68 6.53 + 16 12 4.58 5.28 6.18 7.14 7.11 6.94 + 17 13 4.84 5.59 6.54 7.56 7.55 7.35 + 18 14 5.10 5.90 6.89 7.98 7.98 7.76 + 19 15 5.35 6.21 7.25 8.40 8.42 8.17 + 20 16 5.61 6.52 7.60 8.82 8.85 8.57 + 21 17 5.87 6.83 7.96 9.24 9.28 8.98 + 22 18 6.13 7.14 8.32 9.67 9.72 9.40 + 23 19 6.39 7.45 8.68 10.09 10.16 9.81 + 24 20 6.65 7.76 9.03 10.51 10.59 10.22 + 25 21 6.91 8.07 9.39 10.93 11.03 10.63 + 26 22 7.17 8.38 9.75 11.36 11.47 11.04 + 27 23 7.43 8.69 10.11 11.78 11.90 11.45 + 28 24 7.69 9.00 10.47 12.21 12.34 11.87 + 29 25 7.95 9.31 10.83 12.63 12.78 12.28 + 30 26 8.21 9.63 11.19 13.06 13.22 12.69 + 31 27 8.47 9.94 11.55 13.48 13.66 13.11 + 32 28 8.73 10.25 11.91 13.91 14.10 13.52 + 33 29 8.99 10.56 12.27 14.34 14.54 13.93 + 34 30 9.25 10.88 12.63 14.76 14.98 14.35 + 35 31 9.52 11.19 12.99 15.19 15.42 14.76 + 36 32 9.78 11.50 13.36 15.62 15.86 15.18 + 37 33 10.04 11.82 13.72 16.05 16.30 15.60 + 38 34 10.30 12.13 14.08 16.47 16.74 16.01 + 39 35 10.57 12.45 14.44 16.90 17.19 16.43 + 40 36 10.83 12.76 14.81 17.33 17.63 16.85 + 41 37 11.09 13.08 15.17 17.76 18.07 17.26 + 42 38 11.36 13.39 15.54 18.19 18.52 17.68 + 43 39 11.62 13.71 15.90 18.62 18.96 18.10 + 44 40 11.89 14.03 16.26 19.05 19.41 18.52 + 45 41 12.15 14.34 16.63 19.49 19.85 18.94 + 46 42 12.42 14.66 17.00 19.92 20.30 19.36 + 47 43 12.68 14.98 17.36 20.35 20.75 19.78 + 48 44 12.95 15.29 17.73 20.78 21.19 20.20 + 49 45 13.21 15.61 18.09 21.22 21.64 20.62 + 50 46 13.48 15.93 18.46 21.65 22.09 21.04 + 51 47 13.74 16.25 18.83 22.08 22.54 21.46 + 52 48 14.01 16.57 19.20 22.52 22.98 21.88 + 53 49 14.28 16.89 19.56 22.95 23.43 22.30 + 54-5 50 14.81 17.52 20.30 23.82 24.33 23.15 + 56-7 51 15.35 18.16 21.04 24.70 25.24 24.00 + 58-9 52 15.88 18.81 21.78 25.57 26.14 24.85 + 60-1 53 16.42 19.45 22.52 26.45 27.05 25.70 + 62-3 54 16.96 20.09 23.26 27.33 27.95 26.55 + 64-5 55 17.50 20.74 24.01 28.21 28.86 27.40 + 66-7 56 18.04 21.39 24.75 29.09 29.78 28.26 + 68-9 57 18.58 22.03 25.50 29.97 30.69 29.12 + 70-1 58 19.12 22.68 26.25 30.86 31.61 29.97 + 72-3 59 19.67 23.33 27.00 31.74 32.52 30.83 + 74-5 60 20.21 23.99 27.75 32.63 33.44 31.70 + 76-7 61 20.76 24.64 28.50 33.52 34.36 32.56 + 78-9 62 21.30 25.29 29.26 34.42 35.29 33.43 + 80-1 63 21.85 25.95 30.02 35.31 36.21 34.29 + 82-3 64 22.40 26.61 30.77 36.21 37.14 35.16 + 84-5 65 22.95 27.27 31.53 37.10 38.07 36.03 + 86-7 66 23.50 27.93 32.29 38.00 39.01 36.91 + 88-9 67 24.06 28.59 33.06 38.91 39.94 37.78 + 90-1 68 24.61 29.25 33.82 39.81 40.88 38.66 + 92-3 69 25.17 29.91 34.59 40.72 41.82 39.53 + 94-5 70 25.72 30.58 35.35 41.62 42.76 40.41 + 96-7 71 26.28 31.25 36.12 42.53 43.70 41.30 + 98-9 72 26.84 31.92 36.89 43.44 44.64 42.18 + 100-1 73 27.40 32.59 37.67 44.36 45.59 43.06 + 102-3 74 27.96 33.26 38.44 45.27 46.54 43.95 + 104-5 75 28.52 33.93 39.22 46.19 47.49 44.84 + 106-7 76 29.09 34.61 39.99 47.11 48.45 45.73 + 108-9 77 29.65 35.28 40.77 48.03 49.41 46.62 + 110-1 78 30.22 35.96 41.55 48.96 50.36 47.52 + 112-3 79 30.79 36.64 42.34 49.88 51.33 48.41 + 114-5 80 31.36 37.32 43.12 50.81 52.29 49.31 + 116-7 81 31.93 38.00 43.91 51.74 53.25 50.21 + 118-9 82 32.50 38.68 44.70 52.67 54.22 51.12 + 120-1 83 33.07 39.37 45.49 53.61 55.19 52.02 + 122-3 84 33.64 40.06 46.28 54.54 56.17 52.93 + 124-5 85 34.22 40.74 47.07 55.48 57.14 53.83 + 126-7 86 34.80 41.43 47.87 56.42 58.12 54.74 + 128-9 87 35.38 42.13 48.66 57.37 59.10 55.66 + 130-1 88 35.95 42.82 49.46 58.31 60.08 56.57 + 132-3 89 36.54 43.51 50.26 59.26 61.07 57.49 + 134-5 90 37.12 44.21 51.06 60.21 62.06 58.40 + 136-8 91 37.99 45.26 52.27 61.64 63.54 59.78 + 139-41 92 38.87 46.31 53.48 63.07 65.04 61.17 + 142-4 93 39.75 47.36 54.70 64.51 66.54 62.56 + 145-7 94 40.64 48.42 55.92 65.95 68.04 63.96 + 148-50 95 41.53 49.49 57.14 67.40 69.55 65.36 + 151-3 96 42.42 50.55 58.37 68.86 71.07 66.76 + 154-6 97 43.32 51.62 59.60 70.32 72.59 68.17 + 157-9 98 44.21 52.70 60.84 71.78 74.12 69.59 + 160-2 99 45.12 53.78 62.08 73.25 75.65 71.01 + 163-5 100 46.02 54.86 63.33 74.73 77.19 72.43 + 166-8 101 46.93 55.95 64.58 76.21 78.74 73.86 + 169-71 102 47.84 57.04 65.84 77.70 80.30 75.30 + 172-4 103 48.76 58.13 67.10 79.19 81.86 76.74 + 175-7 104 49.68 59.23 68.37 80.69 83.43 78.18 + 178-80 105 50.60 60.33 69.64 82.20 85.00 79.63 + 181-3 106 51.53 61.44 70.91 83.71 86.58 81.09 + 184-6 107 52.46 62.55 72.19 85.23 88.17 82.55 + 187-9 108 53.39 63.67 73.48 86.75 89.76 84.02 + 190-2 109 54.33 64.79 74.77 88.28 91.37 85.49 + 193-5 110 55.27 65.92 76.07 89.81 92.98 86.97 + 196-8 111 56.22 67.05 77.37 91.36 94.59 88.46 + 199-201 112 57.17 68.18 78.68 92.91 96.22 89.95 + 202-4 113 58.12 69.32 79.99 94.46 97.85 91.45 + 205-7 114 59.08 70.47 81.30 96.02 99.49 92.95 + 208-10 115 60.04 71.62 82.63 97.59 101.14 94.46 + 211-3 116 61.00 72.77 83.96 99.17 102.79 95.97 + 214-6 117 61.97 73.93 85.29 100.75 104.46 97.49 + 217-9 118 62.95 75.10 86.63 102.34 106.13 99.02 + 220-2 119 63.93 76.27 87.98 103.93 107.81 100.55 + 223-5 120 64.91 77.44 89.33 105.54 109.50 102.09 + 226-8 121 65.90 78.62 90.69 107.15 111.19 103.64 + 229-31 122 66.89 79.81 92.05 108.76 112.90 105.19 + 232-4 123 67.89 81.00 93.42 110.39 114.61 106.75 + 235-7 124 68.89 82.20 94.80 112.02 116.33 108.32 + 238-41 125 70.23 83.80 96.64 114.21 118.65 110.42 + 242-5 126 71.58 85.42 98.50 116.41 120.97 112.53 + 246-9 127 72.94 87.04 100.37 118.63 123.32 114.65 + 250-5 128 74.99 89.50 103.19 121.98 126.87 117.86""", + format="ascii", + ) + abs_chan = caltable["Abs"] + chans = [_split_chan_spec(chan, sep="-") for chan in abs_chan] + + col_idx = np.searchsorted(_EDGE_EPOCHS, epoch) + if col_idx == 5: + col_1234 = "E5_1234" + col_0 = "E5_0" + else: + col_1234 = col_0 = f"E{col_idx}" + + # Create a step function for each group of PCUs + energies_0 = np.zeros(256) + energies_1234 = np.zeros(256) + for chan_sep, energy_0, energy_1234 in zip(chans, caltable[col_0], caltable[col_1234]): + energies_0[chan_sep[0] : chan_sep[1] + 1] = energy_0 + energies_1234[chan_sep[0] : chan_sep[1] + 1] = energy_1234 + + def func(chan, detector_id=0): + if detector_id == 0: + return energies_0[int(chan)] + return energies_1234[int(chan)] + + return np.vectorize(func) + + +def rxte_calibration_func(instrument, epoch): + """Return the calibration function for RXTE at a given epoch. + + Examples + -------- + >>> calibration_func = rxte_calibration_func("PCa", 50082) + >>> assert calibration_func(10) == pca_calibration_func(50082)(10) + >>> rxte_calibration_func("HEXTE", 55930) + Traceback (most recent call last): + ... + ValueError: Unknown XTE instrument: HEXTE + """ + if instrument.lower() == "pca": + return pca_calibration_func(epoch) + raise ValueError(f"Unknown XTE instrument: {instrument}") + + +def rxte_pca_event_file_interpretation(hdulist): + """Interpret the FITS header of an RXTE event file. + + At the moment, only science event files are supported. In these files, + the energy channels are stored in a column named PHA. However, this is not + the PHA column that can be directly used to convert to energy. These are + channels that get changed on a per-observation basis, and can be converted + to the "absolute" PHA channels (the ones tabulated in `pca_calibration_func`) + by using the TEVTB2 keyword. This function changes the content of the PHA column by + putting in the mean "absolute" PHA channel corresponding to each local PHA + channel. + + Parameters + ---------- + hdulist : `astropy.io.fits.HDUList` + The FITS file to interpret. + """ + if "XTE_SE" not in hdulist: + raise ValueError( + "No XTE_SE extension found. At the moment, only science events " + "are supported by Stingray for XTE." + ) + tevtb2 = hdulist["XTE_SE"].header["TEVTB2"] + local_chans = np.asarray([int(np.mean(ch)) for ch in _decode_energy_channels(tevtb2)]) + + # channels = _decode_energy_channels(tevtb2) + hdulist["XTE_SE"].data["PHA"] = local_chans[hdulist["XTE_SE"].data["PHA"]] diff --git a/stingray/mission_support/tests/__init__.py b/stingray/mission_support/tests/__init__.py new file mode 100644 index 000000000..4869c602a --- /dev/null +++ b/stingray/mission_support/tests/__init__.py @@ -0,0 +1,4 @@ +# Licensed under MIT license - see LICENSE.rst +""" +This packages contains affiliated package tests. +""" diff --git a/stingray/mission_support/tests/test_mission_support.py b/stingray/mission_support/tests/test_mission_support.py new file mode 100644 index 000000000..9a1f27fd5 --- /dev/null +++ b/stingray/mission_support/tests/test_mission_support.py @@ -0,0 +1,26 @@ +import pytest +import os +from astropy.io import fits +from stingray.mission_support import ( + rxte_pca_event_file_interpretation, + get_rough_conversion_function, +) + +curdir = os.path.abspath(os.path.dirname(__file__)) +datadir = os.path.join(curdir, "..", "..", "tests", "data") + + +class TestAll(object): + def test_wrong_mission_raises(self): + with pytest.raises(ValueError, match="Mission blah not recognized"): + get_rough_conversion_function("blah") + + +class TestXTE(object): + @classmethod + def setup_class(cls): + cls.wrongfile = os.path.join(datadir, "monol_testA.evt") + + def test_wrong_file_raises(self): + with pytest.raises(ValueError, match="No XTE_SE extension found."): + rxte_pca_event_file_interpretation(fits.open(self.wrongfile)) diff --git a/stingray/tests/data/xte_gx_test.evt.gz b/stingray/tests/data/xte_gx_test.evt.gz new file mode 100644 index 000000000..ce0a268ad Binary files /dev/null and b/stingray/tests/data/xte_gx_test.evt.gz differ diff --git a/stingray/tests/data/xte_test.evt.gz b/stingray/tests/data/xte_test.evt.gz new file mode 100644 index 000000000..976bfb574 Binary files /dev/null and b/stingray/tests/data/xte_test.evt.gz differ diff --git a/stingray/tests/test_io.py b/stingray/tests/test_io.py index 417c0dd6b..5f3603e31 100644 --- a/stingray/tests/test_io.py +++ b/stingray/tests/test_io.py @@ -89,7 +89,7 @@ def test_event_file_read_additional_warns_uncal(self): assert np.any(["Column energy not found" in r.message.args[0] for r in record]) # This is the default calibration for nustar data, as returned # from rough_calibration - assert np.allclose(vals.energy_list, vals.pi_list * 0.04 + 1.6) + assert np.allclose(vals.energy_list, vals.pi_list * 0.04 + 1.62) def test_event_file_read_additional_energy_cal(self): """Test event file reading.""" @@ -107,6 +107,22 @@ def test_event_file_read_xmm(self): load_events_and_gtis(fname, additional_columns=["PRIOR"]) assert np.any(["Trying first extension" in r.message.args[0] for r in record]) + def test_event_file_read_rxte(self): + """Test event file reading.""" + fname = os.path.join(datadir, "xte_test.evt.gz") + data = load_events_and_gtis(fname) + assert data.mission.lower() == "xte" + + def test_event_file_read_rxte_gx(self): + """Test event file reading.""" + fname = os.path.join(datadir, "xte_gx_test.evt.gz") + data = load_events_and_gtis(fname) + assert data.mission.lower() == "xte" + assert np.min(data.pi_list) == 2 + assert np.max(data.pi_list) == 255 + assert np.allclose(data.energy_list.min(), 2.06) + assert np.allclose(data.energy_list.max(), 117.86) + def test_event_file_read_no_mission(self): """Test event file reading.""" fname = os.path.join(datadir, "nomission.evt") diff --git a/stingray/utils.py b/stingray/utils.py index ee08796a0..9c59c6901 100644 --- a/stingray/utils.py +++ b/stingray/utils.py @@ -141,9 +141,43 @@ def mad(data, c=0.6745, axis=None): "find_nearest", "check_isallfinite", "heaviside", + "make_dictionary_lowercase", ] +def make_dictionary_lowercase(dictionary, recursive=False): + """Make all keys of a dictionary lowercase. + + Optionally, if some values are dictionaries, they can be made lowercase too. + + Parameters + ---------- + dictionary : dict + The dictionary to be made lowercase + + Other Parameters + ---------------- + recursive : bool + If ``True``, make all keys of nested dictionaries lowercase too. + + Examples + -------- + >>> d1 = {"A": 1, "B": 2, "C": {"D": 3, "E": {"F": 4}}} + >>> d2 = make_dictionary_lowercase(d1) + >>> assert d2 == {"a": 1, "b": 2, "c": {"D": 3, "E": {"F": 4}}} + >>> d3 = make_dictionary_lowercase(d1, recursive=True) + >>> assert d3 == {"a": 1, "b": 2, "c": {"d": 3, "e": {"f": 4}}} + """ + new_dict = {} + for key, value in dictionary.items(): + if recursive and isinstance(value, dict): + value = make_dictionary_lowercase(value, recursive=True) + + new_dict[key.lower()] = value + + return new_dict + + def force_array(x): """Convert an input to a numpy array.