import astropy.units as u
import numpy as np
from astropy.coordinates import AltAz, EarthLocation, SkyCoord
from astropy.io import fits
from astropy.table import QTable
from astropy.time import Time
[docs]
def create_event_list(data, run_number, source_name):
"""
Create the event_list BinTableHDUs from the given data
Parameters
----------
Data: DL2 data file
'astropy.table.QTable'
Run: Run number
Int
Source_name: Name of the source
Str
Returns
-------
Events HDU: `astropy.io.fits.BinTableHDU`
GTI HDU: `astropy.io.fits.BinTableHDU`
Pointing HDU: `astropy.io.fits.BinTableHDU`
"""
DEFAULT_HEADER = fits.Header()
DEFAULT_HEADER["HDUDOC"] = "https://github.com/open-gamma-ray-astro/gamma-astro-data-formats"
DEFAULT_HEADER["HDUVERS"] = "0.2"
DEFAULT_HEADER["HDUCLASS"] = "GADF"
# Timing parameters
lam = 2800 # Average rate of triggered events, taken by hand for now
t_start = data["trigger_time"].value[0]
t_stop = data["trigger_time"].value[-1]
time = Time(data["trigger_time"], format="unix", scale="utc")
date_obs = time[0].to_value("iso", "date")
obs_time = t_stop - t_start # All corrections excluded
# Position parameters
location = EarthLocation.from_geodetic(-17.89139 * u.deg, 28.76139 * u.deg, 2184 * u.m)
reco_alt = data["reco_alt"]
reco_az = data["reco_az"]
pointing_alt = data["pointing_alt"]
pointing_az = data["pointing_az"]
src_sky_pos = SkyCoord(alt=reco_alt, az=reco_az, frame=AltAz(obstime=time, location=location)).transform_to(
frame="icrs"
)
tel_pnt_sky_pos = SkyCoord(
alt=pointing_alt[0], az=pointing_az[0], frame=AltAz(obstime=time[0], location=location)
).transform_to(frame="icrs")
# ~ try:
# ~ object_radec=SkyCoord.from_name(source_name)
# ~ except:
# ~ log.error('Timeout Error in finding Object in Sesame')
object_radec = SkyCoord(tel_pnt_sky_pos.icrs)
# Observation modes
# ~ source_pointing_diff = object_radec.separation(
# ~ SkyCoord(tel_pnt_sky_pos.ra, tel_pnt_sky_pos.dec)
# ~ ).deg
# Assuming wobble offset is fixed to 0.4
# ~ if round(source_pointing_diff, 1) == 0.4:
mode = "WOBBLE"
# ~ elif round(source_pointing_diff, 1) > 1:
# ~ mode = 'OFF'
# ~ elif round(source_pointing_diff, 1) == 0.:
# ~ mode = 'ON'
# ~ else:
# ~ mode = 'ON-MISPOINTING' # Either this or modify the method of getting ON mode
# ~ log.error(f'Source pointing difference with camera pointing is {source_pointing_diff:.3f} deg' )
# Event table
event_table = QTable(
{
"EVENT_ID": u.Quantity(data["event_id"]),
"TIME": u.Quantity(data["trigger_time"]),
"RA": u.Quantity(src_sky_pos.ra.to(u.deg)),
"DEC": u.Quantity(src_sky_pos.dec.to(u.deg)),
"ENERGY": u.Quantity(data["reco_energy"]),
}
)
# GTI table
gti_table = QTable({"START": u.Quantity(t_start, ndmin=1), "STOP": u.Quantity(t_stop, ndmin=1)})
# Adding the meta data
# Event table metadata
ev_header = DEFAULT_HEADER.copy()
ev_header["HDUCLAS1"] = "EVENTS"
ev_header["OBS_ID"] = run_number
ev_header["DATE_OBS"] = date_obs
ev_header["TSTART"] = t_start
ev_header["TSTOP"] = t_stop
ev_header["MJDREFI"] = "40587" # mjd format
ev_header["MJDREFF"] = "0"
ev_header["TIMEUNIT"] = "s"
ev_header["TIMESYS"] = "UTC"
ev_header["OBJECT"] = source_name
ev_header["OBS_MODE"] = mode
ev_header["N_TELS"] = data["tel_id"][0]
ev_header["TELLIST"] = f'LST-{data["tel_id"][0]}'
ev_header["RA_PNT"] = tel_pnt_sky_pos.ra.value
ev_header["DEC_PNT"] = tel_pnt_sky_pos.dec.value
ev_header["ALT_PNT"] = round(np.rad2deg(data["pointing_alt"].value.mean()), 6)
ev_header["AZ_PNT"] = round(np.rad2deg(data["pointing_az"].value[0]), 6)
ev_header["RA_OBJ"] = object_radec.ra.value
ev_header["DEC_OBJ"] = object_radec.dec.value
ev_header["FOVALIGN"] = "ALTAZ"
ev_header["ONTIME"] = obs_time
# Dead time for DRS4 chip is 26 u_sec
ev_header["DEADC"] = 1 / (1 + 2.6e-5 * lam) # 1/(1 + dead_time*lambda)
ev_header["LIVETIME"] = ev_header["DEADC"] * ev_header["ONTIME"]
# GTI table metadata
gti_header = DEFAULT_HEADER.copy()
gti_header["HDUCLAS1"] = "GTI"
gti_header["OBS_ID"] = run_number
gti_header["MJDREFI"] = ev_header["MJDREFI"]
gti_header["MJDREFF"] = ev_header["MJDREFF"]
gti_header["TIMESYS"] = ev_header["TIMESYS"]
gti_header["TIMEUNIT"] = ev_header["TIMEUNIT"]
# Pointing table metadata
pnt_header = DEFAULT_HEADER.copy()
pnt_header["HDUCLAS1"] = "POINTING"
pnt_header["OBS_ID"] = run_number
pnt_header["RA_PNT"] = tel_pnt_sky_pos.ra.value
pnt_header["DEC_PNT"] = tel_pnt_sky_pos.dec.value
pnt_header["ALT_PNT"] = ev_header["ALT_PNT"]
pnt_header["AZ_PNT"] = ev_header["AZ_PNT"]
pnt_header["TIME"] = t_start
# Create HDUs
pnt_table = QTable()
event = fits.BinTableHDU(event_table, header=ev_header, name="EVENTS")
gti = fits.BinTableHDU(gti_table, header=gti_header, name="GTI")
pointing = fits.BinTableHDU(pnt_table, header=pnt_header, name="POINTING")
return event, gti, pointing