Introduce safe intervals in read_events

hendrics/read_events.py

@@ -3,6 +3,7 @@
 
 import os
 import warnings
+from collections.abc import Iterable
 
 import numpy as np
 from stingray.events import EventList
@@ -33,6 +34,7 @@ def treat_event_file(
     additional_columns=None,
     fill_small_gaps=None,
     bin_time_for_occultations=None,
+    safe_interval=None,
 ):
     """Read data from an event file, with no external GTI information.
 
@@ -63,6 +65,9 @@ def treat_event_file(
         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.
+    safe_interval : float or [float, float]
+        A safe interval to exclude at both ends (if single float) or the start
+        and the end (if pair of values) of GTIs.
     """
     # gtistring = assign_value_if_none(gtistring, "GTI,GTI0,STDGTI")
     log.info(f"Opening {filename}")
@@ -85,25 +90,34 @@ def treat_event_file(
     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)
+            log.info("Filtering out occultations")
+            good_gti = create_gti_mask(lc.time, lc.gti, safe_interval=safe_interval)
             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]}")
+                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
+    elif safe_interval is not None:
+        if not isinstance(safe_interval, Iterable):
+            safe_interval = [safe_interval, safe_interval]
+        gti[:, 0] += safe_interval[0]
+        gti[:, 1] -= safe_interval[1]
 
-    lengths = np.array([g1 - g0 for (g0, g1) in gti])
+    lengths = gti[:, 1] - gti[:, 0]
     gti = gti[lengths >= min_length]
     events.gti = gti
     detector_id = events.detector_id
 
     if randomize_by is not None:
-        events.time += np.random.uniform(-randomize_by / 2, randomize_by / 2, events.time.size)
+        events.time += np.random.uniform(
+            -randomize_by / 2, randomize_by / 2, events.time.size
+        )
 
     if fill_small_gaps is not None:
         events = events.fill_bad_time_intervals(fill_small_gaps)
@@ -146,12 +160,15 @@ def treat_event_file(
                 label = "gti"
 
             for ig, g in enumerate(gti_chunks):
-                outfile_local = f"{outroot_local}_{label}{ig:03d}_ev" + HEN_FILE_EXTENSION
+                outfile_local = (
+                    f"{outroot_local}_{label}{ig:03d}_ev" + HEN_FILE_EXTENSION
+                )
 
                 good_gti = cross_two_gtis([g], gti)
                 if noclobber and os.path.exists(outfile_local):
                     warnings.warn(
-                        f"{outfile_local} exists, " + "and noclobber option used. Skipping"
+                        f"{outfile_local} exists, "
+                        + "and noclobber option used. Skipping"
                     )
                     return
                 good = np.logical_and(events.time >= g[0], events.time < g[1])
@@ -271,7 +288,9 @@ def join_eventlists(event_file1, event_file2, new_event_file=None, ignore_instr=
         Output event file
     """
     if new_event_file is None:
-        new_event_file = common_name(event_file1, event_file2) + "_ev" + HEN_FILE_EXTENSION
+        new_event_file = (
+            common_name(event_file1, event_file2) + "_ev" + HEN_FILE_EXTENSION
+        )
 
     events1 = load_events(event_file1)
     events2 = load_events(event_file2)
@@ -344,7 +363,11 @@ def join_many_eventlists(eventfiles, new_event_file=None, ignore_instr=False):
         if not np.isclose(events.mjdref, first_events.mjdref):
             warnings.warn(f"{event_file} has a different MJDREF")
             continue
-        if hasattr(events, "instr") and not events.instr == first_events.instr and not ignore_instr:
+        if (
+            hasattr(events, "instr")
+            and not events.instr == first_events.instr
+            and not ignore_instr
+        ):
             warnings.warn(f"{event_file} is from a different instrument")
             continue
         elif ignore_instr:
@@ -453,11 +476,14 @@ def main_join(args=None):
     import argparse
 
     description = (
-        "Read a cleaned event files and saves the relevant " "information in a standard format"
+        "Read a cleaned event files and saves the relevant "
+        "information in a standard format"
     )
     parser = argparse.ArgumentParser(description=description)
     parser.add_argument("files", help="Files to join", type=str, nargs="+")
-    parser.add_argument("-o", "--output", type=str, help="Name of output file", default=None)
+    parser.add_argument(
+        "-o", "--output", type=str, help="Name of output file", default=None
+    )
     parser.add_argument(
         "--ignore-instr",
         help="Ignore instrument names in channels",
@@ -501,7 +527,8 @@ def main_splitevents(args=None):
     parser.add_argument(
         "--overlap",
         type=float,
-        help="Overlap factor. 0 for no overlap, 0.5 for " "half-interval overlap, and so on.",
+        help="Overlap factor. 0 for no overlap, 0.5 for "
+        "half-interval overlap, and so on.",
         default=None,
     )
     parser.add_argument(
@@ -515,7 +542,9 @@ def main_splitevents(args=None):
 
     if args.split_at_mjd is not None:
         return split_eventlist_at_mjd(args.fname, mjd=args.split_at_mjd)
-    return split_eventlist(args.fname, max_length=args.length_split, overlap=args.overlap)
+    return split_eventlist(
+        args.fname, max_length=args.length_split, overlap=args.overlap
+    )
 
 
 def main(args=None):
@@ -526,7 +555,8 @@ def main(args=None):
     from .base import _add_default_args, check_negative_numbers_in_args
 
     description = (
-        "Read a cleaned event files and saves the relevant " "information in a standard format"
+        "Read a cleaned event files and saves the relevant "
+        "information in a standard format"
     )
     parser = argparse.ArgumentParser(description=description)
     parser.add_argument("files", help="List of files", nargs="+")
@@ -586,6 +616,13 @@ def main(args=None):
         "XMM)",
         default=None,
     )
+    parser.add_argument(
+        "--safe-interval",
+        nargs=2,
+        type=float,
+        default=[0, 0],
+        help="Interval at start and stop of GTIs used" + " for filtering",
+    )
     parser.add_argument(
         "--fill-small-gaps",
         type=float,
@@ -623,6 +660,7 @@ def main(args=None):
             "fill_small_gaps": args.fill_small_gaps,
             "split_by_detector": not args.ignore_detectors,
             "bin_time_for_occultations": args.bin_time_for_occultations,
+            "safe_interval": args.safe_interval,
         }
 
         arglist = [[f, argdict] for f in files]