Repository: StingraySoftware/HENDRICS

Issue:
HENfake problem
Hi, when I first tried to use command HENfake 

HENfake -e xmm_epn_det04_deorb.nc -o test.fit

the system prompted me the following message: 

Traceback (most recent call last):
 File "/home/gianluca/stingray/bin/HENfake", line 8, in <module>
   sys.exit(main())
 File "/home/gianluca/stingray/lib64/python3.6/site-packages/hendrics/fake.py", line 865, in main
   additional_columns=additional_columns,
 File "/home/gianluca/stingray/lib64/python3.6/site-packages/hendrics/fake.py", line 241, in generate_fake_fits_observation
   ev_list = event_list.time
AttributeError: 'tuple' object has no attribute 'time'

Any idea/hint ?

Thanks, 
gianluca 

Discussion:

Patch:
diff --git a/hendrics/__init__.py b/hendrics/__init__.py
index c348b16..9442938 100644
--- a/hendrics/__init__.py
+++ b/hendrics/__init__.py
@@ -24,48 +24,30 @@ if not _ASTROPY_SETUP_:
         HAS_NETCDF = False
         pass
 
-    import copy
-    import numpy as np
-    import stingray.utils
-    from stingray.events import EventList
-
-    # Patch stingray.utils
-    def _root_squared_mean(array):
-        import numpy as np
-
-        array = np.asarray(array)
-        return np.sqrt(np.sum(array ** 2)) / len(array)
-
-    stingray.utils._root_squared_mean = _root_squared_mean
-
-    class _MonkeyPatchedEventList(EventList):
-        def apply_mask(self, mask, inplace=False):  # pragma: no cover
-            """For compatibility with old stingray version.
-
-            Examples
-            --------
-            >>> evt = _MonkeyPatchedEventList(time=[0, 1, 2])
-            >>> newev0 = evt.apply_mask([True, True, False], inplace=False);
-            >>> newev1 = evt.apply_mask([True, True, False], inplace=True);
-            >>> np.allclose(newev0.time, [0, 1])
-            True
-            >>> np.allclose(newev1.time, [0, 1])
-            True
-            >>> evt is newev1
-            True
-            """
-            if inplace:
-                new_ev = self
-            else:
-                new_ev = copy.deepcopy(self)
-            for attr in "time", "energy", "pi", "cal_pi":
-                if hasattr(new_ev, attr) and getattr(new_ev, attr) is not None:
-                    setattr(new_ev, attr, getattr(new_ev, attr)[mask])
-            return new_ev
+    import stingray
+    from .compat import (
+        _MonkeyPatchedEventList,
+        filter_for_deadtime,
+        get_deadtime_mask,
+    )
+
+    from .compat import (
+        prange,
+        array_take,
+        HAS_NUMBA,
+        njit,
+        vectorize,
+        float32,
+        float64,
+        int32,
+        int64,
+    )
 
     try:
-        e = EventList(time=[1, 2, 3])
+        e = stingray.events.EventList(time=[1, 2, 3])
         e.energy = None
         e.apply_mask([True, True, False])
     except (TypeError, AttributeError):
         stingray.events.EventList = _MonkeyPatchedEventList
+        stingray.filters.filter_for_deadtime = filter_for_deadtime
+        stingray.filters.get_deadtime_mask = get_deadtime_mask
diff --git a/hendrics/base.py b/hendrics/base.py
index 0a59555..c87331d 100644
--- a/hendrics/base.py
+++ b/hendrics/base.py
@@ -5,7 +5,6 @@ import sys
 import copy
 import os
 import warnings
-from functools import wraps
 from collections.abc import Iterable
 from pathlib import Path
 import tempfile
@@ -36,62 +35,6 @@ try:
 except ImportError:
     HAS_SKIMAGE = False
 
-try:
-    from numba import jit, njit, prange, vectorize
-    from numba import float32, float64, int32, int64
-    from numba import types
-    from numba.extending import overload_method
-
-    @overload_method(types.Array, "take")  # pragma: no cover
-    def array_take(arr, indices):
-        if isinstance(indices, types.Array):
-
-            def take_impl(arr, indices):
-                n = indices.shape[0]
-                res = np.empty(n, arr.dtype)
-                for i in range(n):
-                    res[i] = arr[indices[i]]
-                return res
-
-            return take_impl
-
-    HAS_NUMBA = True
-except ImportError:
-    HAS_NUMBA = False
-
-    def njit(**kwargs):
-        """Dummy decorator in case jit cannot be imported."""
-
-        def true_decorator(func):
-            @wraps(func)
-            def wrapped(*args, **kwargs):
-                r = func(*args, **kwargs)
-                return r
-
-            return wrapped
-
-        return true_decorator
-
-    jit = njit
-
-    def prange(*args):
-        """Dummy decorator in case jit cannot be imported."""
-        return range(*args)
-
-    class vectorize(object):
-        def __init__(self, *args, **kwargs):
-            pass
-
-        def __call__(self, func):
-            wrapped_f = np.vectorize(func)
-
-            return wrapped_f
-
-    float32 = float64 = int32 = int64 = lambda x, y: None
-
-    array_take = np.take
-
-
 try:
     from tqdm import tqdm as show_progress
 except ImportError:
@@ -100,6 +43,18 @@ except ImportError:
         return a
 
 
+from . import (
+    prange,
+    array_take,
+    HAS_NUMBA,
+    njit,
+    vectorize,
+    float32,
+    float64,
+    int32,
+    int64,
+)
+
 __all__ = [
     "array_take",
     "njit",
diff --git a/hendrics/compat/__init__.py b/hendrics/compat/__init__.py
new file mode 100644
index 0000000..d38a34d
--- /dev/null
+++ b/hendrics/compat/__init__.py
@@ -0,0 +1,16 @@
+from .compatibility import (
+    _MonkeyPatchedEventList,
+    filter_for_deadtime,
+    get_deadtime_mask,
+)
+from .compatibility import (
+    prange,
+    array_take,
+    HAS_NUMBA,
+    njit,
+    vectorize,
+    float32,
+    float64,
+    int32,
+    int64,
+)
diff --git a/hendrics/compat/compatibility.py b/hendrics/compat/compatibility.py
new file mode 100644
index 0000000..d408ca6
--- /dev/null
+++ b/hendrics/compat/compatibility.py
@@ -0,0 +1,390 @@
+import copy
+from functools import wraps
+import numpy as np
+import stingray.utils
+from stingray.events import EventList
+
+try:
+    from numba import jit, njit, prange, vectorize
+    from numba import float32, float64, int32, int64
+    from numba import types
+    from numba.extending import overload_method
+
+    @overload_method(types.Array, "take")  # pragma: no cover
+    def array_take(arr, indices):
+        """Adapt np.take to arrays"""
+        if isinstance(indices, types.Array):
+
+            def take_impl(arr, indices):
+                n = indices.shape[0]
+                res = np.empty(n, arr.dtype)
+                for i in range(n):
+                    res[i] = arr[indices[i]]
+                return res
+
+            return take_impl
+
+    HAS_NUMBA = True
+except ImportError:
+    HAS_NUMBA = False
+
+    def njit(**kwargs):
+        """Dummy decorator in case jit cannot be imported."""
+
+        def true_decorator(func):
+            @wraps(func)
+            def wrapped(*args, **kwargs):
+                r = func(*args, **kwargs)
+                return r
+
+            return wrapped
+
+        return true_decorator
+
+    jit = njit
+
+    def prange(*args):
+        """Dummy decorator in case jit cannot be imported."""
+        return range(*args)
+
+    class vectorize(object):
+        def __init__(self, *args, **kwargs):
+            pass
+
+        def __call__(self, func):
+            wrapped_f = np.vectorize(func)
+
+            return wrapped_f
+
+    float32 = float64 = int32 = int64 = lambda x, y: None
+
+    def array_take(arr, indices):
+        """Adapt np.take to arrays"""
+        return np.take(arr, indices)
+
+    # Patch stingray.utils
+
+
+def _root_squared_mean(array):
+    array = np.asarray(array)
+    return np.sqrt(np.sum(array ** 2)) / len(array)
+
+
+stingray.utils._root_squared_mean = _root_squared_mean
+
+
+def _paralyzable_dead_time(event_list, dead_time):
+    """Apply paralyzable dead time to an event list.
+
+    Parameters
+    ----------
+    event_list : array of floats
+        Event times of arrival
+    dead_time: float
+        Dead time (single value)
+
+    Returns
+    -------
+    output_event_list : array of floats
+        Filtered event times
+    mask : array of bools
+        Final mask, showing all good events in the original event list.
+    """
+    mask = np.ones(len(event_list), dtype=bool)
+    dead_time_end = event_list + dead_time
+    bad = dead_time_end[:-1] > event_list[1:]
+    # Easy: paralyzable case. Here, events coming during dead time produce
+    # more dead time. So...
+    mask[1:][bad] = False
+
+    return event_list[mask], mask
+
+
+@njit()
+def _nonpar_core(event_list, dead_time_end, mask):
+    """Numba-compiled core of the non-paralyzable dead time calculation.
+
+    Parameters
+    ----------
+    event_list : array of floats
+        Event times of arrival
+    dead_time_end : array of floats
+        End of the dead time of each event
+    mask : array of bools
+        Final mask of good events. Initially, all entries must be ``True``
+
+    Return
+    ------
+    mask : array of bools
+        Final mask of good events
+    """
+    for i in range(1, len(event_list)):
+        if event_list[i] < dead_time_end[i - 1]:
+            dead_time_end[i] = dead_time_end[i - 1]
+            mask[i] = False
+    return mask
+
+
+def _non_paralyzable_dead_time(event_list, dead_time):
+    """Apply non-paralyzable dead time to an event list.
+
+    Parameters
+    ----------
+    event_list : array of floats
+        Event times of arrival
+    dead_time: float
+        Dead time (single value)
+
+    Returns
+    -------
+    output_event_list : array of floats
+        Filtered event times
+    mask : array of bools
+        Final mask, showing all good events in the original event list.
+    """
+    event_list_dbl = (event_list - event_list[0]).astype(np.double)
+    dead_time_end = event_list_dbl + np.double(dead_time)
+    mask = np.ones(event_list_dbl.size, dtype=bool)
+    mask = _nonpar_core(event_list_dbl, dead_time_end, mask)
+    return event_list[mask], mask
+
+
+class _DeadtimeFilterOutput(object):
+    uf_events = None
+    is_event = None
+    deadtime = None
+    mask = None
+    bkg = None
+
+
+def get_deadtime_mask(
+    ev_list,
+    deadtime,
+    bkg_ev_list=None,
+    dt_sigma=None,
+    paralyzable=False,
+    return_all=False,
+    verbose=False,
+):
+    """Filter an event list for a given dead time."""
+    additional_output = _DeadtimeFilterOutput()
+
+    # Create the total lightcurve, and a "kind" array that keeps track
+    # of the events classified as "signal" (True) and "background" (False)
+    if bkg_ev_list is not None:
+        tot_ev_list = np.append(ev_list, bkg_ev_list)
+        ev_kind = np.append(
+            np.ones(len(ev_list), dtype=bool),
+            np.zeros(len(bkg_ev_list), dtype=bool),
+        )
+        order = np.argsort(tot_ev_list)
+        tot_ev_list = tot_ev_list[order]
+        ev_kind = ev_kind[order]
+        del order
+    else:
+        tot_ev_list = ev_list
+        ev_kind = np.ones(len(ev_list), dtype=bool)
+
+    additional_output.uf_events = tot_ev_list
+    additional_output.is_event = ev_kind
+    additional_output.deadtime = deadtime
+    additional_output.uf_mask = np.ones(tot_ev_list.size, dtype=bool)
+    additional_output.bkg = tot_ev_list[np.logical_not(ev_kind)]
+
+    if deadtime <= 0.0:
+        if deadtime < 0:
+            raise ValueError("Dead time is less than 0. Please check.")
+        retval = [np.ones(ev_list.size, dtype=bool), additional_output]
+        return retval
+
+    nevents = len(tot_ev_list)
+    all_ev_kind = ev_kind.copy()
+
+    if dt_sigma is not None:
+        deadtime_values = ra.normal(deadtime, dt_sigma, nevents)
+        deadtime_values[deadtime_values < 0] = 0.0
+    else:
+        deadtime_values = np.zeros(nevents) + deadtime
+
+    initial_len = len(tot_ev_list)
+
+    # Note: saved_mask gives the mask that produces tot_ev_list_filt from
+    # tot_ev_list. The same mask can be used to also filter all other arrays.
+    if paralyzable:
+        tot_ev_list_filt, saved_mask = _paralyzable_dead_time(
+            tot_ev_list, deadtime_values
+        )
+
+    else:
+        tot_ev_list_filt, saved_mask = _non_paralyzable_dead_time(
+            tot_ev_list, deadtime_values
+        )
+    del tot_ev_list
+
+    ev_kind = ev_kind[saved_mask]
+    deadtime_values = deadtime_values[saved_mask]
+    final_len = tot_ev_list_filt.size
+    if verbose:
+        log.info(
+            "filter_for_deadtime: "
+            "{0}/{1} events rejected".format(
+                initial_len - final_len, initial_len
+            )
+        )
+
+    retval = saved_mask[all_ev_kind]
+
+    if return_all:
+        # uf_events: source and background events together
+        # ev_kind : kind of each event in uf_events.
+        # bkg : Background events
+        additional_output.uf_events = tot_ev_list_filt
+        additional_output.is_event = ev_kind
+        additional_output.deadtime = deadtime_values
+        additional_output.bkg = tot_ev_list_filt[np.logical_not(ev_kind)]
+        retval = [retval, additional_output]
+
+    return retval
+
+
+def filter_for_deadtime(event_list, deadtime, **kwargs):
+    """Filter an event list for a given dead time.
+
+    This function accepts either a list of times or a
+    `stingray.events.EventList` object.
+
+    For the complete optional parameter list, see `get_deadtime_mask`
+
+    Parameters
+    ----------
+    ev_list : array-like or class:`stingray.events.EventList`
+        The event list
+    deadtime: float
+        The (mean, if not constant) value of deadtime
+
+    Returns
+    -------
+    new_ev_list : class:`stingray.events.EventList` object
+        The filtered event list
+    additional_output : dict
+        See `get_deadtime_mask`
+
+    """
+    # Need to import here to avoid circular imports in the top module.
+    from stingray.events import EventList
+
+    local_retall = kwargs.pop("return_all", False)
+
+    if isinstance(event_list, EventList):
+        retval = event_list.apply_deadtime(
+            deadtime, return_all=local_retall, **kwargs
+        )
+    else:
+        mask, retall = get_deadtime_mask(
+            event_list, deadtime, return_all=True, **kwargs
+        )
+        retval = event_list[mask]
+        if local_retall:
+            retval = [retval, retall]
+
+    return retval
+
+
+class _MonkeyPatchedEventList(EventList):
+    def __init__(self, *args, **kwargs):
+        EventList.__init__(self, *args, **kwargs)
+        if not hasattr(self, "mission"):
+            self.mission = None
+        if not hasattr(self, "instr"):
+            self.instr = None
+        if not hasattr(self, "timesys"):
+            self.timesys = None
+        if not hasattr(self, "timeref"):
+            self.timeref = None
+        if not hasattr(self, "ephem"):
+            self.ephem = None
+
+    def apply_mask(self, mask, inplace=False):  # pragma: no cover
+        """For compatibility with old stingray version.
+
+        Examples
+        --------
+        >>> evt = _MonkeyPatchedEventList(time=[0, 1, 2])
+        >>> newev0 = evt.apply_mask([True, True, False], inplace=False);
+        >>> newev1 = evt.apply_mask([True, True, False], inplace=True);
+        >>> np.allclose(newev0.time, [0, 1])
+        True
+        >>> np.allclose(newev1.time, [0, 1])
+        True
+        >>> evt is newev1
+        True
+        """
+        if inplace:
+            new_ev = self
+        else:
+            new_ev = copy.deepcopy(self)
+        for attr in "time", "energy", "pi", "cal_pi":
+            if hasattr(new_ev, attr) and getattr(new_ev, attr) is not None:
+                setattr(new_ev, attr, getattr(new_ev, attr)[mask])
+        return new_ev
+
+    def apply_deadtime(self, deadtime, inplace=False, **kwargs):
+        """Apply deadtime filter to this event list.
+
+        Additional arguments in ``kwargs`` are passed to `get_deadtime_mask`
+
+        Parameters
+        ----------
+        deadtime : float
+            Value of dead time to apply to data
+        inplace : bool, default False
+            If True, apply the deadtime to the current event list. Otherwise,
+            return a new event list.
+
+        Returns
+        -------
+        new_event_list : `EventList` object
+            Filtered event list. if `inplace` is True, this is the input object
+            filtered for deadtime, otherwise this is a new object.
+        additional_output : object
+            Only returned if `return_all` is True. See `get_deadtime_mask` for
+            more details.
+
+        Examples
+        --------
+        >>> events = np.array([1, 1.05, 1.07, 1.08, 1.1, 2, 2.2, 3, 3.1, 3.2])
+        >>> events = EventList(events)
+        >>> events.pi=np.array([1, 2, 2, 2, 2, 1, 1, 1, 2, 1])
+        >>> events.energy=np.array([1, 2, 2, 2, 2, 1, 1, 1, 2, 1])
+        >>> events.mjdref = 10
+        >>> filt_events, retval = events.apply_deadtime(0.11, inplace=False,
+        ...                                             verbose=False,
+        ...                                             return_all=True)
+        >>> filt_events is events
+        False
+        >>> expected = np.array([1, 2, 2.2, 3, 3.2])
+        >>> np.allclose(filt_events.time, expected)
+        True
+        >>> np.allclose(filt_events.pi, 1)
+        True
+        >>> np.allclose(filt_events.energy, 1)
+        True
+        >>> np.allclose(events.pi, 1)
+        False
+        >>> filt_events = events.apply_deadtime(0.11, inplace=True,
+        ...                                     verbose=False)
+        >>> filt_events is events
+        True
+        """
+        local_retall = kwargs.pop("return_all", False)
+
+        mask, retall = get_deadtime_mask(
+            self.time, deadtime, return_all=True, **kwargs
+        )
+
+        new_ev = self.apply_mask(mask, inplace=inplace)
+
+        if local_retall:
+            new_ev = [new_ev, retall]
+
+        return new_ev
diff --git a/hendrics/efsearch.py b/hendrics/efsearch.py
index 79de737..1cbfb56 100644
--- a/hendrics/efsearch.py
+++ b/hendrics/efsearch.py
@@ -1851,8 +1851,8 @@ def main_accelsearch(args=None):
         results["fmax"] = fmax
         results["zmax"] = zmax
         if hasattr(events, "mission"):
-            results["mission"] = events.mission
-        results["instr"] = events.instr
+            results["mission"] = events.mission.replace(",", "+")
+        results["instr"] = events.instr.replace(",", "+")
         results["mjdref"] = np.double(events.mjdref)
         results["pepoch"] = events.mjdref + results["time"] / 86400.0
 
diff --git a/hendrics/fake.py b/hendrics/fake.py
index 6848502..7e93e0c 100644
--- a/hendrics/fake.py
+++ b/hendrics/fake.py
@@ -7,10 +7,12 @@ import copy
 import numpy as np
 import numpy.random as ra
 from astropy import log
+from astropy.io.fits import Header
 from astropy.logger import AstropyUserWarning
 from stingray.events import EventList
 from stingray.lightcurve import Lightcurve
 from stingray.utils import assign_value_if_none
+from stingray.filters import filter_for_deadtime
 from .io import get_file_format, load_lcurve
 from .io import load_events, save_events, HEN_FILE_EXTENSION
 from .base import _empty, r_in
@@ -19,169 +21,40 @@ from .lcurve import lcurve_from_fits
 from .base import njit, deorbit_events
 
 
-def _paralyzable_dead_time(event_list, dead_time):
-    mask = np.ones(len(event_list), dtype=bool)
-    dead_time_end = event_list + dead_time
-    bad = dead_time_end[:-1] > event_list[1:]
-    # Easy: paralyzable case. Here, events coming during dead time produce
-    # more dead time. So...
-    mask[1:][bad] = False
-
-    return event_list[mask], mask
-
-
-@njit()
-def _nonpar_core(event_list, dead_time_end, mask):
-    for i in range(1, len(event_list)):
-        if event_list[i] < dead_time_end[i - 1]:
-            dead_time_end[i] = dead_time_end[i - 1]
-            mask[i] = False
-    return mask
-
-
-def _non_paralyzable_dead_time(event_list, dead_time):
-    event_list_dbl = (event_list - event_list[0]).astype(np.double)
-    dead_time_end = event_list_dbl + np.double(dead_time)
-    mask = np.ones(event_list_dbl.size, dtype=bool)
-    mask = _nonpar_core(event_list_dbl, dead_time_end, mask)
-    return event_list[mask], mask
-
-
-def filter_for_deadtime(
-    event_list,
-    deadtime,
-    bkg_ev_list=None,
-    dt_sigma=None,
-    paralyzable=False,
-    additional_data=None,
-    return_all=False,
-):
-    """Filter an event list for a given dead time.
-
-    Parameters
-    ----------
-    ev_list : array-like
-        The event list
-    deadtime: float
-        The (mean, if not constant) value of deadtime
-
-    Returns
-    -------
-    new_ev_list : array-like
-        The filtered event list
-    additional_output : dict
-        Object with all the following attributes. Only returned if
-        `return_all` is True
-    uf_events : array-like
-        Unfiltered event list (events + background)
-    is_event : array-like
-        Boolean values; True if event, False if background
-    deadtime : array-like
-        Dead time values
-    bkg : array-like
-        The filtered background event list
-    mask : array-like, optional
-        The mask that filters the input event list and produces the output
-        event list.
-
-    Other Parameters
-    ----------------
-    bkg_ev_list : array-like
-        A background event list that affects dead time
-    dt_sigma : float
-        The standard deviation of a non-constant dead time around deadtime.
-    return_all : bool
-        If True, return the mask that filters the input event list to obtain
-        the output event list.
-
-    """
-    additional_output = _empty()
-    if not isinstance(event_list, EventList):
-        event_list_obj = EventList(event_list)
-    else:
-        event_list_obj = event_list
-
-    ev_list = event_list_obj.time
-    additional_arrays = {}
-    for arr in "pi,energy".split(","):
-        if hasattr(event_list_obj, arr):
-            values = getattr(event_list_obj, arr)
-            if values is not None:
-                additional_arrays[arr] = values
-
-    if deadtime <= 0.0:
-        return copy.copy(event_list)
-
-    # Create the total lightcurve, and a "kind" array that keeps track
-    # of the events classified as "signal" (True) and "background" (False)
-    if bkg_ev_list is not None:
-        Nback = bkg_ev_list.size
-        tot_ev_list = np.append(ev_list, bkg_ev_list)
-        ev_kind = np.append(
-            np.ones(len(ev_list), dtype=bool),
-            np.zeros(len(bkg_ev_list), dtype=bool),
-        )
-        order = np.argsort(tot_ev_list)
-        tot_ev_list = tot_ev_list[order]
-        ev_kind = ev_kind[order]
-        for arr in additional_arrays.keys():
-            newarr = np.append(additional_arrays[arr], np.zeros(Nback))
-            additional_arrays[arr] = newarr[order]
-
-        del order
-    else:
-        tot_ev_list = ev_list
-        ev_kind = np.ones(len(ev_list), dtype=bool)
-
-    nevents = len(tot_ev_list)
-    all_ev_kind = ev_kind.copy()
-
-    if dt_sigma is not None:
-        deadtime_values = ra.normal(deadtime, dt_sigma, nevents)
-    else:
-        deadtime_values = np.zeros(nevents) + deadtime
-
-    initial_len = len(tot_ev_list)
-
-    if paralyzable:
-        tot_ev_list, saved_mask = _paralyzable_dead_time(
-            tot_ev_list, deadtime_values
-        )
-
-    else:
-        tot_ev_list, saved_mask = _non_paralyzable_dead_time(
-            tot_ev_list, deadtime_values
-        )
-
-    ev_kind = ev_kind[saved_mask]
-    deadtime_values = deadtime_values[saved_mask]
-    final_len = len(tot_ev_list)
-    log.info(
-        "filter_for_deadtime: "
-        "{0}/{1} events rejected".format(initial_len - final_len, initial_len)
+def _fill_in_default_information(tbheader):
+    tbheader["OBSERVER"] = "Edwige Bubble"
+    tbheader["COMMENT"] = (
+        "FITS (Flexible Image Transport System) format is"
+        " defined in 'Astronomy and Astrophysics', volume"
+        " 376, page 359; bibcode: 2001A&A...376..359H"
     )
-
-    for arr in additional_arrays.keys():
-        additional_arrays[arr] = additional_arrays[arr][saved_mask][ev_kind]
-
-    retval = EventList(
-        time=tot_ev_list[ev_kind],
-        mjdref=event_list_obj.mjdref,
-        **additional_arrays,
+    tbheader["OBS_ID"] = ("00000000001", "Observation ID")
+    tbheader["TARG_ID"] = (0, "Target ID")
+    tbheader["OBJECT"] = ("Fake X-1", "Name of observed object")
+    tbheader["RA_OBJ"] = (0.0, "[deg] R.A. Object")
+    tbheader["DEC_OBJ"] = (0.0, "[deg] Dec Object")
+    tbheader["RA_NOM"] = (
+        0.0,
+        "Right Ascension used for barycenter corrections",
     )
-
-    if not isinstance(event_list, EventList):
-        retval = retval.time
-
-    if return_all:
-        additional_output.uf_events = tot_ev_list
-        additional_output.is_event = ev_kind
-        additional_output.deadtime = deadtime_values
-        additional_output.mask = saved_mask[all_ev_kind]
-        additional_output.bkg = tot_ev_list[np.logical_not(ev_kind)]
-        retval = [retval, additional_output]
-
-    return retval
+    tbheader["DEC_NOM"] = (0.0, "Declination used for barycenter corrections")
+    tbheader["RA_PNT"] = (0.0, "[deg] RA pointing")
+    tbheader["DEC_PNT"] = (0.0, "[deg] Dec pointing")
+    tbheader["PA_PNT"] = (0.0, "[deg] Position angle (roll)")
+    tbheader["EQUINOX"] = (2.000e03, "Equinox of celestial coord system")
+    tbheader["RADECSYS"] = ("FK5", "Coordinate Reference System")
+    tbheader["TASSIGN"] = ("SATELLITE", "Time assigned by onboard clock")
+    tbheader["TIMESYS"] = ("TDB", "All times in this file are TDB")
+    tbheader["TIMEREF"] = (
+        "SOLARSYSTEM",
+        "Times are pathlength-corrected to barycenter",
+    )
+    tbheader["CLOCKAPP"] = (False, "TRUE if timestamps corrected by gnd sware")
+    tbheader["COMMENT"] = (
+        "MJDREFI+MJDREFF = epoch of Jan 1, 2010, in TT " "time system."
+    )
+    tbheader["TIMEUNIT"] = ("s", "unit for time keywords")
+    return tbheader
 
 
 def generate_fake_fits_observation(
@@ -233,12 +106,23 @@ def generate_fake_fits_observation(
     from astropy.io import fits
     import numpy.random as ra
 
+    inheader = None
     if event_list is None:
         tstart = assign_value_if_none(tstart, 8e7)
         tstop = assign_value_if_none(tstop, tstart + 1025)
         ev_list = sorted(ra.uniform(tstart, tstop, 1000))
+        gti = assign_value_if_none(gti, np.array([[tstart, tstop]]))
     else:
+        if hasattr(event_list, "header") and event_list.header is not None:
+            inheader = Header.fromstring(event_list.header)
         ev_list = event_list.time
+        gti = assign_value_if_none(
+            event_list.gti, np.asarray([[ev_list[0], ev_list[-1]]])
+        )
+        mission = assign_value_if_none(mission, event_list.mission)
+        instr = assign_value_if_none(instr, event_list.instr)
+        tstart = assign_value_if_none(tstart, gti[0, 0])
+        tstop = assign_value_if_none(tstop, gti[-1, 1])
 
     if hasattr(event_list, "pi") and event_list.pi is not None:
         pi = event_list.pi
@@ -250,9 +134,6 @@ def generate_fake_fits_observation(
     else:
         cal_pi = pi / 3
 
-    tstart = assign_value_if_none(tstart, np.floor(ev_list[0]))
-    tstop = assign_value_if_none(tstop, np.ceil(ev_list[-1]))
-    gti = assign_value_if_none(gti, np.array([[tstart, tstop]]))
     filename = assign_value_if_none(filename, "events.evt")
     livetime = assign_value_if_none(livetime, tstop - tstart)
 
@@ -298,32 +179,18 @@ def generate_fake_fits_observation(
     tbhdu.name = "EVENTS"
 
     # ---- Fake lots of information ----
+
     tbheader = tbhdu.header
-    tbheader["OBSERVER"] = "Edwige Bubble"
-    tbheader["COMMENT"] = (
-        "FITS (Flexible Image Transport System) format is"
-        " defined in 'Astronomy and Astrophysics', volume"
-        " 376, page 359; bibcode: 2001A&A...376..359H"
+    tbheader = _fill_in_default_information(tbheader)
+    if inheader is not None:
+        tbheader.update(inheader)
+
+    tbheader["TSTART"] = (
+        tstart,
+        "Elapsed seconds since MJDREF at start of file",
     )
     tbheader["TELESCOP"] = (mission, "Telescope (mission) name")
     tbheader["INSTRUME"] = (instr, "Instrument name")
-    tbheader["OBS_ID"] = ("00000000001", "Observation ID")
-    tbheader["TARG_ID"] = (0, "Target ID")
-    tbheader["OBJECT"] = ("Fake X-1", "Name of observed object")
-    tbheader["RA_OBJ"] = (0.0, "[deg] R.A. Object")
-    tbheader["DEC_OBJ"] = (0.0, "[deg] Dec Object")
-    tbheader["RA_NOM"] = (
-        0.0,
-        "Right Ascension used for barycenter corrections",
-    )
-    tbheader["DEC_NOM"] = (0.0, "Declination used for barycenter corrections")
-    tbheader["RA_PNT"] = (0.0, "[deg] RA pointing")
-    tbheader["DEC_PNT"] = (0.0, "[deg] Dec pointing")
-    tbheader["PA_PNT"] = (0.0, "[deg] Position angle (roll)")
-    tbheader["EQUINOX"] = (2.000e03, "Equinox of celestial coord system")
-    tbheader["RADECSYS"] = ("FK5", "Coordinate Reference System")
-    tbheader["TASSIGN"] = ("SATELLITE", "Time assigned by onboard clock")
-    tbheader["TIMESYS"] = ("TDB", "All times in this file are TDB")
     tbheader["MJDREFI"] = (
         int(mjdref),
         "TDB time reference; Modified Julian Day (int)",
@@ -332,19 +199,6 @@ def generate_fake_fits_observation(
         mjdref - int(mjdref),
         "TDB time reference; Modified Julian Day (frac)",
     )
-    tbheader["TIMEREF"] = (
-        "SOLARSYSTEM",
-        "Times are pathlength-corrected to barycenter",
-    )
-    tbheader["CLOCKAPP"] = (False, "TRUE if timestamps corrected by gnd sware")
-    tbheader["COMMENT"] = (
-        "MJDREFI+MJDREFF = epoch of Jan 1, 2010, in TT " "time system."
-    )
-    tbheader["TIMEUNIT"] = ("s", "unit for time keywords")
-    tbheader["TSTART"] = (
-        tstart,
-        "Elapsed seconds since MJDREF at start of file",
-    )
     tbheader["TSTOP"] = (tstop, "Elapsed seconds since MJDREF at end of file")
     tbheader["LIVETIME"] = (livetime, "On-source time")
     tbheader["TIMEZERO"] = (0.000000e00, "Time Zero")
@@ -380,11 +234,8 @@ def generate_fake_fits_observation(
 
 
 def _read_event_list(filename):
-    if filename is not None:
-        warnings.warn(
-            "Input event lists not yet implemented", AstropyUserWarning
-        )
-    return None, None
+    ev_list = load_events(filename)
+    return ev_list
 
 
 def _read_light_curve(filename):
diff --git a/hendrics/phaseogram.py b/hendrics/phaseogram.py
index bab778b..84b2232 100644
--- a/hendrics/phaseogram.py
+++ b/hendrics/phaseogram.py
@@ -150,7 +150,7 @@ class BasePhaseogram(object):
         self.nph = nph
         if mjdref is None:
             warnings.warn("MJDREF not set. All MJD values will be incorrect.")
-        self.mjdref = assign_value_if_none(mjdref, 0.)
+        self.mjdref = assign_value_if_none(mjdref, 0.0)
         self.ev_times = np.asarray(ev_times)
         self.gti = gti
         if gti is None:


Test patch:
diff --git a/hendrics/tests/test_efsearch.py b/hendrics/tests/test_efsearch.py
index 60940cf..b729f9b 100644
--- a/hendrics/tests/test_efsearch.py
+++ b/hendrics/tests/test_efsearch.py
@@ -478,7 +478,16 @@ class TestEFsearch:
         evfile = self.dum
         with pytest.warns(UserWarning) as record:
             _ = main_zsearch(
-                [evfile, "-f", "9.85", "-F", "9.95", "-n", "2", "--fast",]
+                [
+                    evfile,
+                    "-f",
+                    "9.85",
+                    "-F",
+                    "9.95",
+                    "-n",
+                    "2",
+                    "--fast",
+                ]
             )
         assert np.any(
             [
diff --git a/hendrics/tests/test_fake.py b/hendrics/tests/test_fake.py
index 93e588b..0c79cb5 100644
--- a/hendrics/tests/test_fake.py
+++ b/hendrics/tests/test_fake.py
@@ -15,6 +15,7 @@ from stingray.events import EventList
 from hendrics.tests import _dummy_par
 from hendrics import fake, calibrate, read_events, io
 from hendrics.fake import scramble
+from hendrics.io import load_events
 from hendrics.fold import HAS_PINT
 
 try:
@@ -121,7 +122,7 @@ def test_deadtime_mask_par():
         events, 0.11, bkg_ev_list=bkg_events, paralyzable=True, return_all=True
     )
 
-    assert np.all(filt_events == events[info.mask])
+    assert np.allclose(filt_events, [2.0, 2.2, 3.0])
 
 
 def test_deadtime_conversion():
@@ -239,6 +240,19 @@ class TestFake(object):
         assert os.path.exists(newfile)
         os.remove(newfile)
 
+    def test_fake_fits_input_events_file(self):
+        newfile = "bububu.fit"
+        command = f"-e {self.first_event_file} -o {newfile}"
+        _ = hen.fake.main(command.split())
+        assert os.path.exists(newfile)
+        newfiles = hen.read_events.treat_event_file(newfile)
+        events0 = load_events(newfiles[0])
+        events1 = load_events(self.first_event_file)
+        assert np.allclose(events0.time, events1.time)
+        assert np.allclose(events0.gti, events1.gti)
+
+        os.remove(newfile)
+
     def test_scramble_uncalibrated_events_file_raises(self):
         command = f"{self.first_event_file} -e 3 30"
         with pytest.raises(ValueError):
diff --git a/hendrics/tests/test_fspec.py b/hendrics/tests/test_fspec.py
index 0b4ffa1..4c5ec5a 100644
--- a/hendrics/tests/test_fspec.py
+++ b/hendrics/tests/test_fspec.py
@@ -262,12 +262,14 @@ class TestFullRun(object):
 
     def test_cpds_wrong_norm(self):
         """Test CPDS production."""
-        command = "{0} {1} -f 128 --save-dyn -k CPDS --norm blablabla -o {2}".format(
-            os.path.join(self.datadir, "monol_testA_E3-50_lc")
-            + HEN_FILE_EXTENSION,
-            os.path.join(self.datadir, "monol_testB_E3-50_lc")
-            + HEN_FILE_EXTENSION,
-            os.path.join(self.datadir, "monol_test_E3-50"),
+        command = (
+            "{0} {1} -f 128 --save-dyn -k CPDS --norm blablabla -o {2}".format(
+                os.path.join(self.datadir, "monol_testA_E3-50_lc")
+                + HEN_FILE_EXTENSION,
+                os.path.join(self.datadir, "monol_testB_E3-50_lc")
+                + HEN_FILE_EXTENSION,
+                os.path.join(self.datadir, "monol_test_E3-50"),
+            )
         )
         with pytest.warns(UserWarning) as record:
             hen.fspec.main(command.split())
diff --git a/hendrics/tests/test_phaseogram.py b/hendrics/tests/test_phaseogram.py
index 929127e..669302f 100644
--- a/hendrics/tests/test_phaseogram.py
+++ b/hendrics/tests/test_phaseogram.py
@@ -115,16 +115,16 @@ class TestPhaseogram:
         )
 
     def test_phaseogram_all_defaults(self):
-        times = [0, 3., 4.]
+        times = [0, 3.0, 4.0]
         with pytest.warns(UserWarning) as record:
-            _ = InteractivePhaseogram(times, 1., test=True)
+            _ = InteractivePhaseogram(times, 1.0, test=True)
         assert np.any(["MJDREF not set." in r.message.args[0] for r in record])
 
     def test_binary_all_no_orb(self):
-        times = [0, 3., 4.]
+        times = [0, 3.0, 4.0]
         with pytest.raises(RuntimeError) as excinfo:
             # Missing t0
-            _ = BinaryPhaseogram(times, 1., orbital_period=3., asini=3.)
+            _ = BinaryPhaseogram(times, 1.0, orbital_period=3.0, asini=3.0)
         assert "Please specify all binary parameters" in str(excinfo.value)
 
     def test_phaseogram_input_norm_invalid(self):
diff --git a/hendrics/tests/test_read_events.py b/hendrics/tests/test_read_events.py
index f9839b1..aca691f 100644
--- a/hendrics/tests/test_read_events.py
+++ b/hendrics/tests/test_read_events.py
@@ -251,7 +251,8 @@ class TestReadEvents:
         assert "mjdref" in data
         gtis = data["gti"]
         lengths = np.array([g1 - g0 for (g0, g1) in gtis])
-        assert np.all(lengths <= 100)
+        # add an epsilon for numerical error
+        assert np.all(lengths <= 100 + 1e-7)
 
     def test_split_events(self):
         treat_event_file(self.fits_fileA)