Fix timestamps in DL3 files, fixes #986

lstchain/high_level/hdu_table.py

@@ -34,6 +34,16 @@
 
 wobble_offset = 0.4 * u.deg
 
+LST_EPOCH = Time('2018-10-01T00:00:00', scale='utc')
+
+
+def time_to_fits(time, epoch=LST_EPOCH):
+    '''Convert time to time since epoch
+
+    This is the real elapsed time, i.e. including leap seconds for UTC.
+    '''
+    return (time - epoch).to(u.s)
+
 
 def create_obs_index_hdu(filename_list, fits_dir, obs_index_file, overwrite):
     """
@@ -230,7 +240,7 @@ def create_hdu_index_hdu(
     hdu_index_list.writeto(hdu_index_file, overwrite=overwrite)
 
 
-def get_timing_params(data):
+def get_timing_params(data, epoch=LST_EPOCH):
     """
     Get event lists and retrieve some timing parameters for the DL3 event list
     as a dict
@@ -240,16 +250,19 @@ def get_timing_params(data):
     t_start_iso = time_utc[0].to_value("iso", "date_hms")
     t_stop_iso = time_utc[-1].to_value("iso", "date_hms")
 
+    mjdreff, mjdrefi = np.modf(epoch.mjd)
+
     time_pars = {
-        "t_start": data["dragon_time"].value[0],
-        "t_stop": data["dragon_time"].value[-1],
+        "t_start": time_to_fits(time_utc[0], epoch=epoch),
+        "t_stop": time_to_fits(time_utc[-1], epoch=epoch),
         "t_start_iso": t_start_iso,
         "t_stop_iso": t_stop_iso,
         "date_obs": t_start_iso[:10],
         "time_obs": t_start_iso[11:],
         "date_end": t_stop_iso[:10],
         "time_end": t_stop_iso[11:],
-        "MJDREF": Time("1970-01-01T00:00", scale="utc")
+        "MJDREFI": int(mjdrefi),
+        "MJDREFF": mjdreff,
     }
     return time_pars
 
@@ -336,7 +349,8 @@ def set_expected_pos_to_reco_altaz(data):
     data["reco_az"]  = expected_src_altaz.az
 
 def create_event_list(
-        data, run_number, source_name, source_pos, effective_time, elapsed_time
+    data, run_number, source_name, source_pos, effective_time, elapsed_time,
+    epoch=LST_EPOCH,
 ):
     """
     Create the event_list BinTableHDUs from the given data
@@ -375,7 +389,7 @@ def create_event_list(
     event_table = QTable(
         {
             "EVENT_ID": data["event_id"],
-            "TIME": data["dragon_time"],
+            "TIME": time_to_fits(Time(data["dragon_time"], format='unix')),
             "RA": data["RA"].to(u.deg),
             "DEC": data["Dec"].to(u.deg),
             "ENERGY": data["reco_energy"],
@@ -418,15 +432,18 @@ def create_event_list(
     ev_header["TIME-OBS"] = time_params["time_obs"]
     ev_header["DATE-END"] = time_params["date_end"]
     ev_header["TIME-END"] = time_params["time_end"]
-    ev_header["TSTART"] = time_params["t_start"]
-    ev_header["TSTOP"] = time_params["t_stop"]
-    ev_header["MJDREFI"] = int(time_params["MJDREF"].mjd)
-    ev_header["MJDREFF"] = time_params["MJDREF"].mjd - int(time_params["MJDREF"].mjd)
+    ev_header["TSTART"] = (time_params["t_start"].value, time_params["t_start"].unit)
+    ev_header["TSTOP"] = (time_params["t_stop"].value, time_params["t_stop"].unit)
+    ev_header["MJDREFI"] = time_params["MJDREFI"]
+    ev_header["MJDREFF"] = time_params["MJDREFF"]
     ev_header["TIMEUNIT"] = "s"
     ev_header["TIMESYS"] = "UTC"
     ev_header["TIMEREF"] = "TOPOCENTER"
     ev_header["ONTIME"] = elapsed_time
-    ev_header["TELAPSE"] = time_params["t_stop"] - time_params["t_start"]
+    ev_header["TELAPSE"] = (
+        (time_params["t_stop"] - time_params["t_start"]).to_value(u.s),
+        u.s
+    )
     ev_header["DEADC"] = effective_time / elapsed_time
     ev_header["LIVETIME"] = effective_time
 

lstchain/high_level/tests/test_hdus.py

@@ -1,5 +1,8 @@
 import pytest
+import astropy.units as u
 from astropy.table import Table
+from astropy.time import Time
+import numpy as np
 
 
 @pytest.fixture(scope='session')
@@ -72,3 +75,26 @@ def test_create_obs_hdu_index(tmp_path, dl3_file):
 
     assert "HDU_CLASS" in Table.read(hdu_index).columns
     assert "OBJECT" in Table.read(obs_index).columns
+
+
+def test_get_timing_params():
+    from lstchain.high_level.hdu_table import get_timing_params, LST_EPOCH
+    t = Time(['2020-01-01T00:00:00', '2020-01-01T00:00:01', '2020-01-01T00:00:10'])
+
+    data = Table({
+        'event_id': np.arange(len(t)),
+        'dragon_time': t.unix,
+    })
+
+    params = get_timing_params(data)
+
+    # Converting to a single mjd number loses precision but should be better than 0.1 us
+    epoch = Time(params["MJDREFI"], params["MJDREFF"], format="mjd")
+    assert (epoch + params['t_start']).isclose(t[0], atol=0.1 * u.us)
+    assert (epoch + params['t_stop']).isclose(t[-1], atol=0.1 * u.us)
+
+    assert params["date_obs"] == "2020-01-01"
+    assert params["date_end"] == "2020-01-01"
+    assert params["time_obs"] == "00:00:00.000"
+    assert params["time_end"] == "00:00:10.000"
+    assert epoch == LST_EPOCH