Add calibration with RMF file

docs/changes/804.feature.rst

@@ -0,0 +1 @@
+Add function to calibrate event lists based on RMF file

stingray/events.py

@@ -14,7 +14,7 @@
 from .base import StingrayTimeseries
 from .filters import get_deadtime_mask
 from .gti import generate_indices_of_boundaries
-from .io import load_events_and_gtis
+from .io import load_events_and_gtis, pi_to_energy
 from .lightcurve import Lightcurve
 from .utils import simon, njit
 from .utils import histogram
@@ -142,6 +142,9 @@ class EventList(StingrayTimeseries):
     ephem : str
         The JPL ephemeris used to barycenter the data, if any (e.g. DE430)
 
+    rmf_file : str, default None
+        The file name of the RMF file to use for calibration.
+
     **other_kw :
         Used internally. Any other keyword arguments will be ignored
 
@@ -207,6 +210,7 @@ def __init__(
         ephem=None,
         timeref=None,
         timesys=None,
+        rmf_file=None,
         **other_kw,
     ):
         if ncounts is not None:
@@ -215,6 +219,12 @@ def __init__(
                 DeprecationWarning,
             )
 
+        if rmf_file is not None:
+            if pi is None:
+                warnings.warn("PI channels must be provided to calibrate the energy")
+            else:
+                energy = pi_to_energy(pi, rmf_file)
+
         StingrayTimeseries.__init__(
             self,
             time=time,
@@ -232,6 +242,7 @@ def __init__(
             ephem=ephem,
             timeref=timeref,
             timesys=timesys,
+            rmf_file=rmf_file,
             **other_kw,
         )
 
@@ -559,7 +570,7 @@ def join(self, other, strategy="infer"):
         return self._join_timeseries(other, strategy=strategy, ignore_meta=["header", "ncounts"])
 
     @classmethod
-    def read(cls, filename, fmt=None, **kwargs):
+    def read(cls, filename, fmt=None, rmf_file=None, **kwargs):
         r"""Read a :class:`EventList` object from file.
 
         Currently supported formats are
@@ -589,6 +600,11 @@ def read(cls, filename, fmt=None, **kwargs):
 
         Other parameters
         ----------------
+        rmf_file : str, default None
+            The file name of the RMF file to use for energy calibration. Defaults to
+            None, which implies no channel->energy conversion at this stage (or a default
+            calibration applied to selected missions).
+
         kwargs : dict
             Any further keyword arguments to be passed to `load_events_and_gtis`
             for reading in event lists in OGIP/HEASOFT format
@@ -620,10 +636,26 @@ def read(cls, filename, fmt=None, **kwargs):
                 for key in evtdata.additional_data:
                     if not hasattr(evt, key.lower()):
                         setattr(evt, key.lower(), evtdata.additional_data[key])
+            if rmf_file is not None:
+                evt.convert_pi_to_energy(rmf_file)
             return evt
 
         return super().read(filename=filename, fmt=fmt)
 
+    def convert_pi_to_energy(self, rmf_file):
+        """Calibrate the energy column of the event list.
+
+        Defines the ``energy`` attribute of the event list by converting the
+        PI channels to energy using the provided RMF file.
+
+        Parameters
+        ----------
+        rmf_file : str
+            The file name of the RMF file to use for calibration.
+        """
+
+        self.energy = pi_to_energy(self.pi, rmf_file)
+
     def get_energy_mask(self, energy_range, use_pi=False):
         """Get a mask corresponding to events with a given energy range.
 

stingray/io.py

@@ -11,6 +11,7 @@
 from astropy.table import Table
 from astropy.logger import AstropyUserWarning
 import matplotlib.pyplot as plt
+from astropy.io import fits as pf
 
 import stingray.utils as utils
 from stingray.loggingconfig import setup_logger
@@ -38,10 +39,64 @@
 logger = setup_logger()
 
 
+def read_rmf(rmf_file):
+    """Load RMF info.
+
+    .. note:: Preliminary: only EBOUNDS are read.
+
+    Parameters
+    ----------
+    rmf_file : str
+        The rmf file used to read the calibration.
+
+    Returns
+    -------
+    pis : array-like
+        the PI channels
+    e_mins : array-like
+        the lower energy bound of each PI channel
+    e_maxs : array-like
+        the upper energy bound of each PI channel
+    """
+
+    with pf.open(rmf_file, checksum=True, memmap=False) as lchdulist:
+        lchdulist.verify("warn")
+        lctable = lchdulist["EBOUNDS"].data
+        pis = np.array(lctable.field("CHANNEL"))
+        e_mins = np.array(lctable.field("E_MIN"))
+        e_maxs = np.array(lctable.field("E_MAX"))
+
+    return pis, e_mins, e_maxs
+
+
+def pi_to_energy(pis, rmf_file):
+    """Read the energy channels corresponding to the given PI channels.
+
+    Parameters
+    ----------
+    pis : array-like
+        The channels to lookup in the rmf
+
+    Other Parameters
+    ----------------
+    rmf_file : str
+        The rmf file used to read the calibration.
+    """
+    calp, cal_emin, cal_emax = read_rmf(rmf_file)
+    es = np.zeros(len(pis), dtype=float)
+    for ic, c in enumerate(calp):
+        good = pis == c
+        if not np.any(good):
+            continue
+        es[good] = (cal_emin[ic] + cal_emax[ic]) / 2
+
+    return es
+
+
 def rough_calibration(pis, mission):
     """Make a rough conversion between PI channel and energy.
 
-    Only works for NICER, NuSTAR, and XMM.
+    Only works for NICER, NuSTAR, IXPE, and XMM.
 
     Parameters
     ----------
@@ -710,6 +765,13 @@ def load_events_and_gtis(
     else:
         try:
             returns.energy_list = rough_calibration(cal_pi, mission)
+            logger.info(
+                f"A default calibration was applied to the {mission} data. "
+                "See io.rough_calibration for details. "
+                "Use the `rmf_file` argument in `EventList.read`, or calibrate with "
+                "`EventList.convert_pi_to_energy(rmf_file)`, if you want to apply a specific "
+                "response matrix"
+            )
         except ValueError:
             returns.energy_list = None
     returns.instr = instr.lower()