diff --git a/hendrics/read_events.py b/hendrics/read_events.py index 4fa2afb1..68c649c4 100644 --- a/hendrics/read_events.py +++ b/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]