Source code for rta_selection.dl3_writer

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