# Copyright (C) 2022, the ixpeobssim team.
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
"""Basic data format definitions for the event lists.
"""
from __future__ import print_function, division
from ixpeobssim.core.fitsio import xPrimaryHDU, xBinTableHDUBase
from ixpeobssim.instrument.du import du_physical_name, du_logical_name
from ixpeobssim.utils.astro import xy_columns_kwargs, set_xy_header_limits, build_wcs
from ixpeobssim.utils.logging_ import logger
from ixpeobssim.utils.time_ import MISSION_START_MJD, MISSION_START_MJDREFF
from ixpeobssim.utils.time_ import met_to_string
from ixpeobssim.utils.units_ import arcsec_to_degrees, degrees_to_arcmin
# pylint: disable=invalid-name, too-many-ancestors, too-many-arguments, no-member
# Header keywords related to the telescope.
_TELESCOPE_HEADER_KEYWORDS = [
('TELESCOP', 'IXPE', 'Telescope name'),
('INSTRUME', 'GPD', 'Instrument name'),
('DETNAM' , 'N/A', 'Name of the logical detector unit'),
('DET_ID' , 'N/A', 'Name of the physical detector unit')
]
# Header keywords related to time.
_TIME_HEADER_KEYWORDS = [
('TSTART' , -1, 'Observation start time in MET'),
('TSTOP' , -1, 'Observation end time in MET'),
('DATE-OBS', 'N/A', 'Observation start datetime'),
('DATE-END', 'N/A', 'Observation end datetime'),
('TELAPSE' , -1, 'TSTOP-TSTART'),
('TIMESYS' , 'TT', 'Time system'),
('TIMEUNIT', 'seconds', 'Time units'),
('TIMEREF' , 'LOCAL', 'Time reference'),
('MJDREFI' , MISSION_START_MJD, 'MJD ref day 01 Jan 2017 00:00:00 UTC'),
('MJDREFF' , MISSION_START_MJDREFF, 'Frac part of MJD ref (32.184secs+37leapsecs)'),
('TIMEZERO', 0., 'Zero time'),
('ONTIME' , -1, 'On source time'),
('LIVETIME', -1, 'On source time corrected for dead time'),
('DEADC' , -1, 'Dead time correction'),
('DEADAPP' , 'F', 'Has DEADC been applied to data')
]
#Header keywords related to the object being observed.
_OBJECT_HEADER_KEYWORDS = [
('RA_OBJ' , 'N/A', '[deg] R.A. Object'),
('DEC_OBJ', 'N/A', '[deg] Dec Object'),
('RA_PNT' , 'N/A', '[deg] RA pointing'),
('DEC_PNT', 'N/A', '[deg] Dec pointing'),
('OBJECT' , 'N/A', 'Name of observed object')
]
COMMON_HEADER_KEYWORDS = \
_TELESCOPE_HEADER_KEYWORDS + _TIME_HEADER_KEYWORDS + _OBJECT_HEADER_KEYWORDS
#Header keywords related to the WCS.
_WCS_HEADER_KEYWORDS = [
('EQUINOX' , 2000, 'equinox of celestial coord system'),
('RADECSYS', 'FK5', 'celestial coord system')
]
#Header keywords related to file.
_VERSION_HEADER_KEYWORDS = [
('FILE_LVL', 'LV2', 'File level'),
('LV2_VER' , -1 , 'Version of the LV2 data format')
]
[docs]
def set_version_keywords(hdu, version):
"""Set the version-related header keywords for a given HDU.
"""
logger.info('Updating version header keywords for the %s HDU...', hdu.name)
header = hdu.header
header.set('LV2_VER', version)
# Basic parameters for the output sky coordinates.
_SKYCOORD_NUM_SIDE_PIXELS = 600
_SKYCOORD_PIXEL_SIZE = arcsec_to_degrees(2.6)
[docs]
def standard_xy_columns_kwargs(ra0, dec0):
"""Return the appropriate keywords for the X and Y columns in event files.
See https://bitbucket.org/ixpesw/ixpeobssim/issues/552
"""
return xy_columns_kwargs(ra0, dec0, _SKYCOORD_NUM_SIDE_PIXELS, _SKYCOORD_PIXEL_SIZE)
[docs]
def set_standard_xy_header_limits(hdu):
"""Set the TLMIN and TLMAX keywords for the X and Y columns for a given hdu.
"""
set_xy_header_limits(hdu, _SKYCOORD_NUM_SIDE_PIXELS)
[docs]
def build_standard_wcs(ra0, dec0):
"""Build the standard WCS object for event files programmatically.
See https://bitbucket.org/ixpesw/ixpeobssim/issues/552
"""
return build_wcs(ra0, dec0, _SKYCOORD_NUM_SIDE_PIXELS, _SKYCOORD_PIXEL_SIZE)
# Definition of the WCS origin.
WCS_ORIGIN = 1
[docs]
def standard_radec_to_xy(ra, dec, ra0, dec0):
"""Convert sky coordinates to X and Y digital coordinates in the sky frame.
This is builfding on the fly a standard IXPE WCS object centered at the given
sky position, and using it for the conversion.
Args
----
ra : array_like
The array of input Ra coordinates.
dec : array_like
The array of input Dec coordinates.
ra0 : float
The Ra coordinate of the center of the field in the sky.
dec0 : float
The Dec coordinate of the center of the field in the sky.
"""
return build_standard_wcs(ra0, dec0).wcs_world2pix(ra, dec, WCS_ORIGIN)
[docs]
def standard_xy_to_radec(x, y, ra0, dec0):
"""Convert X and Y digital coordinates to RA and DEC.
This is builfding on the fly a standard IXPE WCS object centered at the given
sky position, and using it for the conversion.
Args
----
x : array_like
The array of input X coordinates.
y : array_like
The array of input Y coordinates.
ra0 : float
The Ra coordinate of the center of the field in the sky.
dec0 : float
The Dec coordinate of the center of the field in the sky.
"""
return build_standard_wcs(ra0, dec0).wcs_pix2world(x, y, WCS_ORIGIN)
[docs]
class xLvl2PrimaryHDU(xPrimaryHDU):
"""Level 2 primary header definition.
"""
HEADER_KEYWORDS = [
('OBS_ID' , -1, 'Observation ID'),
('CONTNUMB', -1, 'Contact number'),
('OBS_MODE', 'OBSERVATION SIMULATION', 'Observation mode'),
('SRC_CONF', 'ASTRO', 'Source configuration'),
('ORIGIN' , 'IXPE team', 'Organization responsible for the data'),
('CALDBVER', -1, 'CALDB version'),
('CLOCKCOR', 'UNKNOWN', 'Whether the time has been corrected')
] + COMMON_HEADER_KEYWORDS + _VERSION_HEADER_KEYWORDS
[docs]
class xBinTableHDUEvents(xBinTableHDUBase):
"""Binary table description for the EVENTS extension of the observation
output files.
"""
NAME = 'EVENTS'
HEADER_KEYWORDS = COMMON_HEADER_KEYWORDS + _WCS_HEADER_KEYWORDS + _VERSION_HEADER_KEYWORDS
DATA_SPECS = [
('TRG_ID' , 'J', None , 'Trigger identifier'),
('SEC' , 'J', 's' , 'Integral part of event time (MET)'),
('MICROSEC', 'J', 'us' , 'Fractional part of event time (MET)'),
('TIME' , 'D', 's' , 'Event time in seconds (MET)'),
('LIVETIME', 'J', 'us' , 'Live time since the previous event in microseconds'),
('PHA' , 'J', None , 'Event pulse height'),
('PI' , 'E', None , 'Event pulse invariant'),
('ENERGY' , 'E', 'keV' , 'Event energy in keV'),
('NUM_CLU' , 'I', '' , 'Number of clusters in the event'),
('DETX' , 'E', 'mm' , 'Reconstructed absorption point X (GPD frame)'),
('DETY' , 'E', 'mm' , 'Reconstructed absorption point Y (GPD frame)'),
('RA' , 'E', 'deg' , 'Event right ascension'),
('DEC' , 'E', 'deg' , 'Event declination'),
('X' , 'E', 'pixel', 'Event X position (SKY frame)'),
('Y' , 'E', 'pixel', 'Event Y position (SKY frame)'),
('DETPHI' , 'E', 'rad' , 'Photolectron emission angle (GPD frame)'),
('PHI' , 'E', 'rad' , 'Photoelectron emission angle (SKY frame)'),
('Q' , 'E', None , 'Corrected event q Stokes parameter'),
('U' , 'E', None , 'Corrected event u Stokes parameter'),
('W_MOM' , 'E', None , 'Event weight')
]
def __init__(self, ra0, dec0, data=None, keywords=None, comments=None):
"""Overloaded constructor.
We need this in order to update the DATA_KWARGS class member, so that
when the constructor of the base class is called, the proper keyword
arguments are passed to the X and Y FITS column creation.
"""
xkwargs, ykwargs = standard_xy_columns_kwargs(ra0, dec0)
self.DATA_KWARGS = {'X': xkwargs, 'Y': ykwargs}
xBinTableHDUBase.__init__(self, data, keywords, comments)
# Note that we also have to set the limits for the X and Y columns for the
# thing to be properly displayed in ds9.
set_standard_xy_header_limits(self)
[docs]
class xBinTableHDUMonteCarlo(xBinTableHDUBase):
"""Binary table description for the MONTE_CARLO extension of the observation
output files, including the additional Monte Carlo fields.
"""
NAME = 'MONTE_CARLO'
HEADER_KEYWORDS = [
('IRFNAME' , 'N/A', 'name of the IRFs used for the simulation')
] + COMMON_HEADER_KEYWORDS
DATA_SPECS = [
('SRC_ID' , 'I', None , 'Monte Carlo source identifier'),
('MC_ENERGY', 'E', 'keV' , 'Monte Carlo event energy'),
('MC_PHA' , 'J', None , 'Monte Carlo pulse height'),
('MC_PI' , 'E', None , 'Monte Carlo pulse invariant'),
('MC_RA' , 'E', 'degrees', 'Monte Carlo right ascension'),
('MC_DEC' , 'E', 'degrees', 'Monte Carlo declination'),
('MC_X' , 'I', 'degrees', 'Monte Carlo event X position (SKY frame)'),
('MC_Y' , 'I', 'degrees', 'Monte Carlo event Y position (SKY frame)'),
('MC_GAIN' , 'E', None , 'Relative GEM gain used for the event')
]
[docs]
def set_irf_name(self, irf_name):
"""Set the IRFNAME keyword.
"""
self.set_keyword('IRFNAME', irf_name)
[docs]
class xBinTableHDUGTI(xBinTableHDUBase):
"""Binary table for the good time intervals (GTI).
"""
NAME = 'GTI'
HEADER_KEYWORDS = COMMON_HEADER_KEYWORDS
DATA_SPECS = [
('START', 'D', 's', 'GTI start time'),
('STOP' , 'D', 's', 'GTI stop time')
]
[docs]
class xBinTableHDURoiTable(xBinTableHDUBase):
"""Binary table for the good time intervals (GTI).
"""
NAME = 'ROITABLE'
HEADER_KEYWORDS = [
('ROIRA' , -1, 'right ascension of the ROI center'),
('ROIDEC' , -1, 'declination of the ROI center'),
('EQUINOX' , 2000., 'equinox for RA and DEC')
] + COMMON_HEADER_KEYWORDS
DATA_SPECS = [
('SRCID' , 'I' , None, 'source identifier'),
('SRCNAME', 'A20', None, 'source name')
]
[docs]
def set_center(self, ra0, dec0):
"""Set the keywords for the ROI center.
"""
self.set_keyword('ROIRA', ra0)
self.set_keyword('ROIDEC', dec0)
[docs]
class xBinTableHDUSpacecraftData(xBinTableHDUBase):
"""Binary table for the spacecraft data.
"""
NAME = 'SC_DATA'
HEADER_KEYWORDS = [
('ROLL' , -1, 'spacecraft roll angle'),
] + COMMON_HEADER_KEYWORDS
DATA_SPECS = [
('MET' , 'D', 's' , 'Mission elapsed time'),
('RA_PNT' , 'E', 'degrees', 'Pointing RA'),
('DEC_PNT' , 'E', 'degrees', 'Pointing DEC'),
('LAT_GEO' , 'E', 'degrees', 'Spacecraft latitude'),
('LON_GEO' , 'E', 'degrees', 'Spacecraft longitude'),
('ALT_GEO' , 'E', 'km' , 'Spacecraft elevation'),
('SUN_ANGLE' , 'E', 'degrees', 'Angle between the Sun and the target'),
('IN_SAA' , 'I', None , 'SAA flag'),
('TARGET_OCCULT', 'I', None , 'Earth occultation flag')
]
[docs]
def set_roll_angle(self, roll_angle):
"""Set the roll angle.
"""
self.set_keyword('ROLL', roll_angle)
[docs]
class xBinTableHDUTimeline(xBinTableHDUBase):
"""Binary table for the timeline data.
"""
NAME = 'TIMELINE'
HEADER_KEYWORDS = COMMON_HEADER_KEYWORDS
DATA_SPECS = [
('START' , 'D', 's' , 'Epoch start time'),
('STOP' , 'D', 's' , 'Epoch stop time'),
('IN_SAA' , 'I', None , 'SAA flag'),
('TARGET_OCCULT', 'I', None , 'Earth occultation flag')
]
[docs]
class xBinTableHDUOCTI(xBinTableHDUBase):
"""Binary table for the on-orbit calibration time intervals (OCTIs).
"""
NAME = 'OCTI'
HEADER_KEYWORDS = [
('CALRUNS', -1, 'Number of on-orbit calibration runs'),
('CALTIME', -1, 'Total time for on-orbit calibration in s')
] + COMMON_HEADER_KEYWORDS
DATA_SPECS = [
('START', 'D', 's', 'OCTI start time'),
('STOP' , 'D', 's', 'OCTI stop time')
]
[docs]
def set_cal_stats(self, num_runs, total_time):
"""Set the calibration statistics keywords.
"""
self.set_keyword('CALRUNS', num_runs)
self.set_keyword('CALTIME', total_time)