Merge pull request #162 from StingraySoftware/options_for_xte

Add option to avoid splitting by detector

hendrics/fspec.py

@@ -93,13 +93,21 @@ def average_periodograms(fspec_iterable, total=None):
 
 def _wrap_fun_cpds(arglist):
     f1, f2, outname, kwargs = arglist
-    return calc_cpds(f1, f2, outname=outname, **kwargs)
+    try:
+        return calc_cpds(f1, f2, outname=outname, **kwargs)
+    except Exception as e:
+        log.error(f"Error in {f1}/{f2}: {e}")
+        return None
 
 
 def _wrap_fun_pds(argdict):
     fname = argdict["fname"]
     argdict.pop("fname")
-    return calc_pds(fname, **argdict)
+    try:
+        return calc_pds(fname, **argdict)
+    except Exception as e:
+        log.error(f"Error in {fname}: {e}")
+        return None
 
 
 def sync_gtis(lc1, lc2):

hendrics/plot.py

@@ -791,10 +791,20 @@ def plot_lc(
         lc = lcdata.counts
         gti = lcdata.gti
         instr = lcdata.instr
+        if not hasattr(lcdata, "mjdref") or lcdata.mjdref is None:
+            lcdata.mjdref = 0
+
+        mjd = lcdata.mjdref + np.mean(lcdata.time) / 86400
+        mjd = np.round(mjd, 2)
 
         if fromstart:
-            time -= lcdata.gti[0, 0]
-            gti -= lcdata.gti[0, 0]
+            tref = lcdata.gti[0, 0]
+            mjd = lcdata.gti[0, 0] / 86400 + lcdata.mjdref
+        else:
+            tref = (mjd - lcdata.mjdref) * 86400
+
+        time -= tref
+        gti -= tref
 
         if instr == "PCA":
             # If RXTE, plot per PCU count rate
@@ -819,7 +829,8 @@ def plot_lc(
                 outqdpdata.append(lcdata.base[good])
             save_as_qdp(outqdpdata, filename=output_data_file, mode="a")
 
-    plt.xlabel("Time (s)")
+    plt.xlabel(f"Time (s from MJD {mjd}, MET {tref})")
+    print(f"Time (s from MJD {mjd}, MET {tref})")
     if instr == "PCA":
         plt.ylabel("light curve (Ct/bin/PCU)")
     else:

hendrics/read_events.py

@@ -9,6 +9,7 @@
 from stingray.io import load_events_and_gtis
 from stingray.events import EventList
 from stingray.gti import cross_two_gtis, cross_gtis, check_separate
+from stingray.gti import create_gti_from_condition, create_gti_mask
 from .io import load_events
 from .base import common_name
 from .base import hen_root
@@ -25,8 +26,10 @@ def treat_event_file(
     length_split=None,
     randomize_by=None,
     discard_calibration=False,
+    split_by_detector=True,
     additional_columns=None,
     fill_small_gaps=None,
+    bin_time_for_occultations=None,
 ):
     """Read data from an event file, with no external GTI information.
 
@@ -49,8 +52,14 @@ def treat_event_file(
         length_split is not None)
     discard_calibration: bool
         discard the automatic calibration done by Stingray (if any)
+    split_by_detector: bool, default True
+        split data from different detectors
     fill_small_gaps: float
         fill gaps between GTIs with random data, if shorter than this amount
+    bin_time_for_occultations: float
+        Create a light curve with this bin time and infer occultations not
+        recorded in GTIs. The rule is that the flux drops to zero and the average
+        counts per bin are significantly above 25 ct/s.
     """
     # gtistring = assign_value_if_none(gtistring, "GTI,GTI0,STDGTI")
     log.info("Opening %s" % filename)
@@ -70,6 +79,23 @@ def treat_event_file(
     if hasattr(events, "instr") and isinstance(events.instr, str):
         instr = events.instr.lower()
     gti = events.gti
+    if bin_time_for_occultations is not None and bin_time_for_occultations > 0:
+        lc = events.to_lc(bin_time_for_occultations)
+        meanrate = np.median(lc.counts)
+
+        if meanrate > 25:
+            good_gti = create_gti_mask(lc.time, lc.gti)
+            good = lc.counts > 0
+            new_bad = (~good) & good_gti
+            if np.any(new_bad):
+                warnings.warn(
+                    f"Found zero counts in the light curve at times{lc.time[new_bad]}"
+                )
+                gti = create_gti_from_condition(
+                    lc.time, (good_gti & good), safe_interval=bin_time_for_occultations
+                )
+                events.gti = gti
+
     lengths = np.array([g1 - g0 for (g0, g1) in gti])
     gti = gti[lengths >= min_length]
     events.gti = gti
@@ -83,7 +109,7 @@ def treat_event_file(
     if fill_small_gaps is not None:
         events = events.fill_bad_time_intervals(fill_small_gaps)
 
-    if detector_id is not None:
+    if detector_id is not None and split_by_detector:
         detectors = np.array(list(set(detector_id)))
     else:
         detectors = [None]
@@ -542,6 +568,12 @@ def main(args=None):
         default=False,
         action="store_true",
     )
+    parser.add_argument(
+        "--ignore-detectors",
+        help="Do not split by detector",
+        default=False,
+        action="store_true",
+    )
     parser.add_argument(
         "-l",
         "--length-split",
@@ -555,6 +587,15 @@ def main(args=None):
         help="Minimum length of GTIs to consider",
         default=0,
     )
+    parser.add_argument(
+        "--bin-time-for-occultations",
+        type=float,
+        help=(
+            "Create a light curve with this bin time and infer occultations not recorded in GTIs."
+            " (The flux drops to zero and the average count rate is significantly above 25 ct/s)"
+        ),
+        default=None,
+    )
     parser.add_argument("--gti-string", type=str, help="GTI string", default=None)
     parser.add_argument(
         "--randomize-by",
@@ -599,6 +640,8 @@ def main(args=None):
             "discard_calibration": args.discard_calibration,
             "additional_columns": args.additional,
             "fill_small_gaps": args.fill_small_gaps,
+            "split_by_detector": not args.ignore_detectors,
+            "bin_time_for_occultations": args.bin_time_for_occultations,
         }
 
         arglist = [[f, argdict] for f in files]

hendrics/tests/test_read_events.py

@@ -18,6 +18,9 @@
 from hendrics.fake import main
 import hendrics as hen
 from . import cleanup_test_dir
+from astropy.utils import minversion
+
+STINGRAY_LT_2_0 = not minversion(np, "2.0.0.dev")
 
 
 class TestMergeEvents:
@@ -202,6 +205,8 @@ class TestReadEvents:
     def setup_class(cls):
         curdir = os.path.abspath(os.path.dirname(__file__))
         cls.datadir = os.path.join(curdir, "data")
+        cls.fits_file_xte = os.path.join(cls.datadir, "xte_test.evt.gz")
+        cls.fits_file_xte_gx = os.path.join(cls.datadir, "xte_gx_test.evt.gz")
         cls.fits_fileA = os.path.join(cls.datadir, "monol_testA.evt")
         cls.fits_fileB = os.path.join(cls.datadir, "monol_testB.evt")
         cls.fits_file = os.path.join(cls.datadir, "monol_test_fake.evt")
@@ -252,6 +257,43 @@ def test_treat_event_file_xmm(self):
         assert "mjdref" in data
         assert "pi" in data and data["pi"].size > 0
 
+    def test_treat_event_file_xte_se(self):
+        treat_event_file(
+            self.fits_file_xte, split_by_detector=False, bin_time_for_occultations=1
+        )
+        new_filename = "xte_test_xte_pca_ev" + HEN_FILE_EXTENSION
+        assert os.path.exists(os.path.join(self.datadir, new_filename))
+        data = load_data(os.path.join(self.datadir, new_filename))
+        assert "instr" in data and data["instr"].lower() == "pca"
+        assert "gti" in data
+        assert "mjdref" in data
+        assert "pi" in data and data["pi"].size > 0
+        if not STINGRAY_LT_2_0:
+            assert "detector_id" in data and set(data["detector_id"]) == {2, 4}
+
+    @pytest.mark.skipif(STINGRAY_LT_2_0, reason="Stingray > 2.0.0 required")
+    def test_treat_event_file_xte_se_split(self):
+        treat_event_file(self.fits_file_xte, split_by_detector=True)
+        new_filename = "xte_test_xte_pca_det02_ev" + HEN_FILE_EXTENSION
+        assert os.path.exists(os.path.join(self.datadir, new_filename))
+        data = load_data(os.path.join(self.datadir, new_filename))
+        assert "instr" in data and data["instr"].lower() == "pca"
+        assert "gti" in data
+        assert "mjdref" in data
+        assert "pi" in data and data["pi"].size > 0
+
+    def test_treat_event_file_xte_gx(self):
+        treat_event_file(
+            self.fits_file_xte_gx, split_by_detector=False, bin_time_for_occultations=1
+        )
+        new_filename = "xte_gx_test_xte_pca_ev" + HEN_FILE_EXTENSION
+        assert os.path.exists(os.path.join(self.datadir, new_filename))
+        data = load_data(os.path.join(self.datadir, new_filename))
+        assert "instr" in data and data["instr"].lower() == "pca"
+        assert "gti" in data
+        assert "mjdref" in data
+        assert "pi" in data and data["pi"].size > 0
+
     def test_treat_event_file_xmm_gtisplit(self):
         treat_event_file(self.fits_file, gti_split=True)
         new_filename = "monol_test_fake_xmm_epn_det01_gti000_ev" + HEN_FILE_EXTENSION