1087 lines
47 KiB
Python
1087 lines
47 KiB
Python
# drm.py: Detector response matrix class
|
|
#
|
|
# Authors: William Cleveland (USRA),
|
|
# Adam Goldstein (USRA) and
|
|
# Daniel Kocevski (NASA)
|
|
#
|
|
# Portions of the code are Copyright 2020 William Cleveland and
|
|
# Adam Goldstein, Universities Space Research Association
|
|
# All rights reserved.
|
|
#
|
|
# Written for the Fermi Gamma-ray Burst Monitor (Fermi-GBM)
|
|
#
|
|
# 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, see <https://www.gnu.org/licenses/>.
|
|
#
|
|
import os
|
|
import astropy.io.fits as fits
|
|
import numpy as np
|
|
from collections import OrderedDict
|
|
from scipy.interpolate import interp1d
|
|
from scipy.integrate import trapz
|
|
|
|
from gbm.time import Met
|
|
from gbm.data.primitives import Bins
|
|
from .data import DataFile
|
|
from . import headers as hdr
|
|
|
|
|
|
class RSP(DataFile):
|
|
"""Class for single- and multi-DRM response files
|
|
|
|
Attributes:
|
|
channel_centroids (np.array): Geometric centroids of the energy channel bins
|
|
channel_widths (np.array): Widths of the energy channel bins
|
|
datatype (str): The datatype of the file
|
|
detector (str): The GBM detector the file is associated with
|
|
directory (str): The directory the file is located in
|
|
ebounds (:obj:`np.recarray`): The energy channel bounds
|
|
filename (str): The filename
|
|
full_path (str): The full path+filename
|
|
headers (dict): The headers of the RSP file
|
|
id (str): The GBM file ID
|
|
is_gbm_file (bool): True if the file is a valid GBM standard file,
|
|
False if it is not.
|
|
is_trigger (bool): True if the file is a GBM trigger file, False if not
|
|
numchans (int): Number of energy channels
|
|
numdrms (int): Number of DRMs
|
|
numebins (int): Number of photon bins
|
|
photon_bins (np.array, np.array): The photon bin edges
|
|
photon_bin_centroids (np.array): Geometric centroids of the photon bins
|
|
photon_bin_widths (np.array): Widths of the photon bins
|
|
trigtime (float): The trigger time, if applicable
|
|
tcent (np.array): The center times of the intervals over which the DRMs are valid
|
|
tstart (np.array): The start times of the intervals over which the DRMs are valid
|
|
tstop (np.array): The end times of the intervals over which the DRMs are valid
|
|
"""
|
|
|
|
def __init__(self):
|
|
super(RSP, self).__init__()
|
|
self._headers = OrderedDict()
|
|
self._ebounds = None
|
|
self._pbounds = None
|
|
self._drm_list = []
|
|
self._ngrp_list = []
|
|
self._fchan_list = []
|
|
self._nchan_list = []
|
|
|
|
@property
|
|
def numchans(self):
|
|
return self._drm_list[0].shape[1]
|
|
|
|
@property
|
|
def numebins(self):
|
|
return self._drm_list[0].shape[0]
|
|
|
|
@property
|
|
def numdrms(self):
|
|
return len(self._drm_list)
|
|
|
|
@property
|
|
def trigtime(self):
|
|
if 'TRIGTIME' in self._headers['PRIMARY'].keys():
|
|
trigtime = self._headers['PRIMARY']['TRIGTIME']
|
|
if trigtime == '':
|
|
trigtime = None
|
|
return trigtime
|
|
|
|
@property
|
|
def tstart(self):
|
|
return self._tstart
|
|
|
|
@property
|
|
def tstop(self):
|
|
return self._tstop
|
|
|
|
@property
|
|
def tcent(self):
|
|
return (self.tstart + self.tstop) / 2.0
|
|
|
|
@property
|
|
def headers(self):
|
|
return self._headers
|
|
|
|
@property
|
|
def ebounds(self):
|
|
return self._ebounds
|
|
|
|
@property
|
|
def photon_bins(self):
|
|
return self._pbounds[0], self._pbounds[1]
|
|
|
|
@property
|
|
def photon_bin_centroids(self):
|
|
return np.sqrt(self._pbounds[0] * self._pbounds[1])
|
|
|
|
@property
|
|
def photon_bin_widths(self):
|
|
return self._pbounds[1] - self._pbounds[0]
|
|
|
|
@property
|
|
def channel_centroids(self):
|
|
return np.sqrt(self._ebounds['E_MIN'] * self._ebounds['E_MAX'])
|
|
|
|
@property
|
|
def channel_widths(self):
|
|
return self._ebounds['E_MAX'] - self._ebounds['E_MIN']
|
|
|
|
@classmethod
|
|
def open(cls, filename):
|
|
"""Read a RSP or RSP2 file from disk
|
|
|
|
Args:
|
|
filename (str): The filename
|
|
|
|
Returns:
|
|
:class:`RSP` The RSP object
|
|
"""
|
|
obj = cls()
|
|
obj._file_properties(filename)
|
|
|
|
# open FITS file
|
|
with fits.open(filename) as hdulist:
|
|
# read the headers and valid times (can be any number of drms)
|
|
obj._tstart = []
|
|
obj._tstop = []
|
|
ispec = 1
|
|
for hdu in hdulist:
|
|
if hdu.name == 'SPECRESP MATRIX':
|
|
name = hdu.name + str(ispec)
|
|
ispec += 1
|
|
obj._tstart.append(hdu.header['TSTART'])
|
|
obj._tstop.append(hdu.header['TSTOP'])
|
|
else:
|
|
name = hdu.name
|
|
obj._headers.update({name: hdu.header})
|
|
|
|
# set the time ranges of each DRM
|
|
obj._tstart = np.array(obj._tstart)
|
|
obj._tstop = np.array(obj._tstop)
|
|
if obj.trigtime is not None:
|
|
obj._tstart -= obj.trigtime
|
|
obj._tstop -= obj.trigtime
|
|
|
|
# now pull out the actual data
|
|
obj._ebounds = hdulist['EBOUNDS'].data
|
|
num_drms = obj._headers['PRIMARY']['DRM_NUM']
|
|
for idrm in range(1, num_drms + 1):
|
|
drm_data = hdulist['SPECRESP MATRIX', idrm].data
|
|
obj._ngrp_list.append(drm_data['N_GRP'])
|
|
obj._fchan_list.append(drm_data['F_CHAN'])
|
|
obj._nchan_list.append(drm_data['N_CHAN'])
|
|
obj._pbounds = (drm_data['ENERG_LO'], drm_data['ENERG_HI'])
|
|
decompressed = obj._decompress_drm(drm_data, idrm)
|
|
obj._drm_list.append(decompressed)
|
|
|
|
return obj
|
|
|
|
@classmethod
|
|
def from_arrays(cls, chan_lo, chan_hi, energies_lo, energies_hi, eff_area,
|
|
trigtime=None, tstart=None, tstop=None, detnam=None,
|
|
object=None, ra_obj=None, dec_obj=None, err_rad=None,
|
|
infile=None, ch2e_ver=None, gain=None, mat_type=None,
|
|
src_az=None, src_el=None, geo_az=None, geo_el=None,
|
|
det_ang=None, geo_ang=None, atscat=None, datatype=None,
|
|
filetype='DRM', filename=None):
|
|
"""Create a single-DRM RSP object from arrays of channel energies,
|
|
incident photon energies and effective area.
|
|
|
|
Args:
|
|
chan_lo (np.array): The low edges of the energy channels
|
|
chan_hi (np.array): The high edges of the energy channels
|
|
energies_lo (np.array): The low edges of the photon bins
|
|
energies_hi (np.array): The high edges of the photon bins
|
|
eff_area (np.array): The effective area of each photon bin mapped
|
|
to each energy channel
|
|
trigtime (float, optional): The trigger time, if applicable.
|
|
tstart (float, optional): The start time of the observation for
|
|
which the response is valid
|
|
tstop (float, optional): The stop time of the observation for
|
|
which the response is valid
|
|
detnam (str, optional): The detector name the corresponds to the response
|
|
object (str, optional): The object being observed
|
|
ra_obj (float, optional): The RA of the object
|
|
dec_obj (float, optional): The Dec of the object
|
|
err_rad (float, optional): The localization error radius of the object
|
|
infile (str, optional): The file from which the response was made
|
|
ch2e_ver (str, optional): The version of the channel to energy calibration
|
|
gain (float, optional): The gain, if applicable
|
|
mat_type (int, optional): If set: 0=instrument scattering only,
|
|
1=atmospheric scattering only,
|
|
2=instrument+atmospheric scattering
|
|
src_az (float, optional): The azimuth in spacecraft coordinates for
|
|
which the response is valid
|
|
src_el (float, optional): The elevation in spacecraft coordinates for
|
|
which the response is valid
|
|
geo_az (float, optional): The location of the geocenter in spacecraft azimuth
|
|
geo_el (float, optional): The location of the geocenter in spacecraft elevation
|
|
det_ang (float, optional): The angle between detector normal and
|
|
position for which the response is valid
|
|
geo_ang (float, optional): The angle between the geocenter and the
|
|
position for which the response is valid
|
|
atscat (str, optional): The atmospheric scattering database used
|
|
datatype (str, optional): The datatype to which the response corresponds
|
|
filetype (str, optional): If not set, default is 'DRM'
|
|
filename (str, optional): If not set, a default name will be
|
|
constructed from other inputs
|
|
|
|
Returns:
|
|
:class:`RSP`: The newly created RSP object
|
|
"""
|
|
obj = cls()
|
|
|
|
# set up the photon bounds and ebounds
|
|
obj._pbounds = (np.unique(energies_lo), np.unique(energies_hi))
|
|
nbounds = obj._pbounds[0].shape[0]
|
|
|
|
chan_lo_unique = np.unique(chan_lo)
|
|
chan_hi_unique = np.unique(chan_hi)
|
|
nchannels = chan_lo_unique.size
|
|
obj._ebounds = np.recarray(nchannels,
|
|
dtype=[('CHANNEL', '>i2'), ('E_MIN', '>f4'),
|
|
('E_MAX', '>f4')])
|
|
obj._ebounds['CHANNEL'] = np.arange(nchannels)
|
|
obj._ebounds['E_MIN'] = chan_lo_unique
|
|
obj._ebounds['E_MAX'] = chan_hi_unique
|
|
|
|
# put the matrix in the correct format
|
|
obj._drm_list = [eff_area.reshape(nbounds, -1)]
|
|
|
|
# Set dummy compression information, equal to no compression at all
|
|
obj._ngrp_list = [np.ones(nbounds)]
|
|
obj._fchan_list = [[np.ones(1)] * nbounds]
|
|
obj._nchan_list = [[nchannels * np.ones(1)] * nbounds]
|
|
|
|
# set the metadata and headers
|
|
obj._tstart = np.array([tstart])
|
|
obj._tstop = np.array([tstop])
|
|
if trigtime is not None:
|
|
obj._tstart -= trigtime
|
|
obj._tstop -= trigtime
|
|
|
|
pri_hdr = hdr.primary(detnam=detnam, filetype=filetype, tstart=tstart,
|
|
tstop=tstop, filename=filename,
|
|
trigtime=trigtime, object=object, ra_obj=ra_obj,
|
|
dec_obj=dec_obj, err_rad=err_rad, infile=infile)
|
|
pri_hdr.append(('DRM_NUM', 1, 'Number of DRMs stored in this file'))
|
|
|
|
eb_hdr = hdr.ebounds(detnam=detnam, tstart=tstart, tstop=tstop,
|
|
trigtime=trigtime, object=object, ra_obj=ra_obj,
|
|
dec_obj=dec_obj, err_rad=err_rad,
|
|
detchans=obj.numchans, ch2e_ver=ch2e_ver,
|
|
gain=gain,
|
|
infile=infile)
|
|
|
|
spec_hdr = hdr.specresp(detnam=detnam, tstart=tstart, tstop=tstop,
|
|
trigtime=trigtime, object=object,
|
|
ra_obj=ra_obj,
|
|
dec_obj=dec_obj, mat_type=mat_type, rsp_num=1,
|
|
src_az=src_az, src_el=src_el, geo_az=geo_az,
|
|
geo_el=geo_el, det_ang=det_ang,
|
|
geo_ang=geo_ang,
|
|
numebins=obj.numebins, detchans=obj.numchans,
|
|
atscat=atscat)
|
|
|
|
obj._headers = {'PRIMARY': pri_hdr, 'EBOUNDS': eb_hdr,
|
|
'SPECRESP MATRIX': spec_hdr}
|
|
|
|
# set file properties
|
|
if (trigtime is None) and (tstart is None):
|
|
tstart = 0.0
|
|
obj.set_properties(detector=detnam, trigtime=obj.trigtime,
|
|
tstart=tstart, datatype=str(datatype).lower(),
|
|
extension='rsp')
|
|
|
|
return obj
|
|
|
|
def effective_area(self, energy, index=None, atime=None,
|
|
interp_kind='linear'):
|
|
"""Calculate the effective area at a given energy or an array of
|
|
energies. This function interpolates the DRM (in log10(cm^2)) to
|
|
provide the effective area for any energy within the photon energy
|
|
bounds of the DRM.
|
|
|
|
Args:
|
|
energy (float or np.array): The photon energy or energies at which
|
|
to calculate the effective area.
|
|
index (int, optional): The DRM index. Default is 0
|
|
atime (float, optional): The time corresponding to a DRM.
|
|
The nearest DRM to the time will be chosen
|
|
interp_kind (str, optional): The kind of interpolation to be
|
|
passed to scipy.interp1d. Default is
|
|
'linear'.
|
|
|
|
Returns:
|
|
(np.array): The effective area
|
|
"""
|
|
if not isinstance(energy, np.ndarray):
|
|
try:
|
|
energy = np.array([energy], dtype=float).squeeze()
|
|
except:
|
|
raise TypeError('energy must be a positive float')
|
|
if (energy <= 0.0).sum() > 0:
|
|
raise ValueError('energy must be positive')
|
|
|
|
if (atime is None) and (index is None):
|
|
index = 0
|
|
|
|
effarea = self.photon_effective_area(index=index, atime=atime)
|
|
diff_effarea = effarea.counts
|
|
|
|
|
|
# create the DRM interpolator
|
|
x = np.append(self.photon_bins[0], self.photon_bins[1][-1])
|
|
y = np.append(diff_effarea, diff_effarea[-1])
|
|
nz_mask = (y > 0)
|
|
effarea_interpolator = interp1d(x[nz_mask], np.log10(y[nz_mask]),
|
|
kind=interp_kind,
|
|
fill_value='extrapolate')
|
|
|
|
# do the interpolation
|
|
effarea = 10.0**effarea_interpolator(energy)
|
|
return effarea
|
|
|
|
def extract_drm(self, index=None, atime=None):
|
|
"""Extract a single DRM from a multi-DRM RSP2 object and return
|
|
a new single-DRM RSP object.
|
|
|
|
Args:
|
|
index (int, optional): The DRM index to retrieve
|
|
atime (float, optional): The time corresponding to a DRM.
|
|
The nearest DRM to the time will be chosen
|
|
|
|
Returns:
|
|
:class:`RSP` : The single-DRM RSP object
|
|
"""
|
|
if (atime is None) and (index is None):
|
|
index = 0
|
|
if atime is not None:
|
|
index = self.drm_index((atime, atime))[0]
|
|
drm = self._drm_list[index]
|
|
|
|
drm_ext = 'SPECRESP MATRIX{}'.format(index + 1)
|
|
atscat = None
|
|
if 'INFILE04' in self.headers[drm_ext]:
|
|
atscat = self.headers[drm_ext]['INFILE04']
|
|
obj = RSP.from_arrays(self._ebounds['E_MIN'], self._ebounds['E_MAX'],
|
|
*self._pbounds, drm, trigtime=self.trigtime,
|
|
tstart=self.tstart[index] + self.trigtime,
|
|
tstop=self.tstop[index] + self.trigtime,
|
|
detnam=self.detector,
|
|
object=self.headers['PRIMARY']['OBJECT'],
|
|
ra_obj=self.headers['PRIMARY']['RA_OBJ'],
|
|
dec_obj=self.headers['PRIMARY']['DEC_OBJ'],
|
|
infile=self.filename,
|
|
ch2e_ver=self.headers['EBOUNDS']['CH2E_VER'],
|
|
gain=self.headers['EBOUNDS']['GAIN_COR'],
|
|
mat_type=self.headers[drm_ext]['MAT_TYPE'],
|
|
src_az=self.headers[drm_ext]['SRC_AZ'],
|
|
src_el=self.headers[drm_ext]['SRC_EL'],
|
|
geo_az=self.headers[drm_ext]['GEO_AZ'],
|
|
geo_el=self.headers[drm_ext]['GEO_EL'],
|
|
det_ang=self.headers[drm_ext]['DET_ANG'],
|
|
geo_ang=self.headers[drm_ext]['GEO_ANG'],
|
|
atscat=atscat, datatype=self.datatype)
|
|
|
|
obj._ngrp_list.append(self._ngrp_list[index])
|
|
obj._fchan_list.append(self._fchan_list[index])
|
|
obj._nchan_list.append(self._nchan_list[index])
|
|
|
|
return obj
|
|
|
|
def extract_weighted_drm(self, weights, bin_centers):
|
|
raise NotImplementedError('Not implemented yet.')
|
|
|
|
def interpolate(self, atime, **kwargs):
|
|
"""Interpolate over a multi-DRM object for a given time and return
|
|
a new single-DRM RSP object. This function uses `scipy.interpolate.interp1d
|
|
<https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html>`_.
|
|
|
|
If the given time is before (after) the times covered by the DRMs
|
|
within the RSP object, then the first (last) DRM will be returned.
|
|
|
|
Args:
|
|
atime (float): The time corresponding to a DRM. The nearest DRM to
|
|
the time will be chosen
|
|
**kwargs: Keyword arguments to be passed to scipy.interpolate.interp1d
|
|
|
|
Returns:
|
|
:class:`RSP`: The interpolated single-DRM RSP object
|
|
"""
|
|
if self.numdrms == 1:
|
|
raise ValueError('Single DRM response. Cannot interpolate')
|
|
|
|
# create the interpolate and interpolate
|
|
tcent = (self.tstart + self.tstop) / 2.0
|
|
drm_interp = interp1d(tcent, np.array(self._drm_list), axis=0,
|
|
fill_value=(self._drm_list[0], self._drm_list[-1]),
|
|
**kwargs, bounds_error=False)
|
|
new_drm = drm_interp(atime)
|
|
|
|
# create the object
|
|
index = self.drm_index((atime, atime))[0]
|
|
drm_ext = 'SPECRESP MATRIX{}'.format(index + 1)
|
|
atscat = None
|
|
if 'INFILE04' in self.headers[drm_ext]:
|
|
atscat = self.headers[drm_ext]['INFILE04']
|
|
obj = RSP.from_arrays(self._ebounds['E_MIN'], self._ebounds['E_MAX'],
|
|
*self._pbounds, new_drm, trigtime=self.trigtime,
|
|
tstart=atime + self.trigtime,
|
|
tstop=atime + self.trigtime,
|
|
detnam=self.detector,
|
|
object=self.headers['PRIMARY']['OBJECT'],
|
|
ra_obj=self.headers['PRIMARY']['RA_OBJ'],
|
|
dec_obj=self.headers['PRIMARY']['DEC_OBJ'],
|
|
infile=self.filename,
|
|
ch2e_ver=self.headers['EBOUNDS']['CH2E_VER'],
|
|
gain=self.headers['EBOUNDS']['GAIN_COR'],
|
|
mat_type=self.headers[drm_ext]['MAT_TYPE'],
|
|
src_az=self.headers[drm_ext]['SRC_AZ'],
|
|
src_el=self.headers[drm_ext]['SRC_EL'],
|
|
geo_az=self.headers[drm_ext]['GEO_AZ'],
|
|
geo_el=self.headers[drm_ext]['GEO_EL'],
|
|
det_ang=self.headers[drm_ext]['DET_ANG'],
|
|
geo_ang=self.headers[drm_ext]['GEO_ANG'],
|
|
atscat=atscat, datatype=self.datatype)
|
|
|
|
obj._ngrp_list.append(self._ngrp_list[index])
|
|
obj._fchan_list.append(self._fchan_list[index])
|
|
obj._nchan_list.append(self._nchan_list[index])
|
|
|
|
return obj
|
|
|
|
def fold_spectrum(self, function, params, index=0, atime=None,
|
|
channel_mask=None):
|
|
"""Fold a photon spectrum through a DRM to get a count spectrum
|
|
|
|
Args:
|
|
function (<function>):
|
|
A photon spectrum function. The function must accept a list of
|
|
function parameters as its first argument and an array of photon
|
|
energies as its second argument. The function must return the
|
|
evaluation of the photon model in units of ph/s-cm^2-keV.
|
|
params (list of float): A list of parameter values to be passed to
|
|
the photon spectrum function
|
|
index (int, optional): The index of the DRM to be used from the list
|
|
of DRMs in the response. Default is 0.
|
|
atime (float, optional):
|
|
A time associated with a DRM in the list of available DRMs.
|
|
The DRM covering the time closest to atime will be used.
|
|
If set, this overrides the index argument.
|
|
channel_mask (np.array, optional):
|
|
A boolean mask where True indicates the channel is to be used
|
|
for folding and False indicates the channel is to not be used
|
|
for folding. If omitted, all channels are used.
|
|
|
|
Returns:
|
|
np.array: The count spectrum
|
|
"""
|
|
if atime is not None:
|
|
drm = self.nearest_drm(atime)
|
|
else:
|
|
drm = self._drm_list[index]
|
|
if channel_mask is not None:
|
|
drm = drm[:, channel_mask]
|
|
|
|
# evaluate photon model
|
|
photon_model = function(params, self.photon_bin_centroids)
|
|
|
|
# fold photon model through DRM
|
|
counts = np.dot(drm.T, photon_model * self.photon_bin_widths)
|
|
|
|
return counts
|
|
|
|
def drm_index(self, time_range):
|
|
"""Return the indices of the DRMs that cover the requested time range
|
|
|
|
If the time range is before (after) the times covered by the DRMs in
|
|
the RSP object, then the first (last) index will be returned.
|
|
|
|
Args:
|
|
time_range (float, float): The time range
|
|
|
|
Returns:
|
|
np.array: An array of the indices of the DRMs covering the time range
|
|
"""
|
|
start, stop = self._assert_range(time_range)
|
|
mask = (self.tstop > start) & (self.tstart < stop)
|
|
|
|
if mask.sum() == 0:
|
|
if stop < self.tstart[0]:
|
|
return np.array([0], dtype=int)
|
|
else:
|
|
return np.array([self.numdrms-1], dtype=int)
|
|
|
|
idx = np.arange(self.numdrms, dtype=int)
|
|
return idx[mask]
|
|
|
|
def nearest_drm(self, atime):
|
|
"""Return the nearest DRM to the requested time
|
|
|
|
Args:
|
|
atime (float): The requested time
|
|
|
|
Returns:
|
|
np.array: The nearest DRM
|
|
"""
|
|
# only one DRM in the file; return it
|
|
if self.numdrms == 1:
|
|
return self._drm_list[0]
|
|
|
|
idx = self.drm_index((atime, atime))[0]
|
|
drm = self._drm_list[idx]
|
|
return drm
|
|
|
|
def drm(self, index):
|
|
"""Return the DRM based on the index into the DRM list
|
|
|
|
Args:
|
|
index (int): The DRM index
|
|
|
|
Returns:
|
|
np.array: The requested DRM
|
|
"""
|
|
return self._drm_list[index]
|
|
|
|
def photon_effective_area(self, index=None, atime=None):
|
|
"""Returns the effective area as a function of incident photon energy
|
|
(integrated over recorded energy channels).
|
|
|
|
Args:
|
|
index (int, optional): The DRM index
|
|
atime (float, optional): The time corresponding to a DRM
|
|
|
|
Returns:
|
|
:class:`~.primitives.Bins`: A Bins object containing the effective area histogram
|
|
"""
|
|
if index is not None:
|
|
drm = self.drm(index)
|
|
elif atime is not None:
|
|
drm = self.nearest_drm(atime)
|
|
else:
|
|
raise ValueError("Must set either 'index' or 'atime'")
|
|
|
|
bins = Bins(drm.sum(axis=1), *self.photon_bins)
|
|
return bins
|
|
|
|
def channel_effective_area(self, index=None, atime=None):
|
|
"""Returns the effective area as a function of recorded channel energy
|
|
(integrated over incident photon bins).
|
|
|
|
Args:
|
|
index (int, optional): The DRM index
|
|
atime (float, optional): The time corresponding to a DRM
|
|
|
|
Returns:
|
|
:class:`~.primitives.Bins`: A Bins object containing the effective area histogram
|
|
"""
|
|
if index is not None:
|
|
drm = self.drm(index)
|
|
elif atime is not None:
|
|
drm = self.nearest_drm(atime)
|
|
else:
|
|
raise ValueError("Must set either 'index' or 'atime'")
|
|
|
|
bins = Bins(drm.sum(axis=0), self.ebounds['E_MIN'],
|
|
self.ebounds['E_MAX'])
|
|
return bins
|
|
|
|
def weighted(self, weights, bin_centers):
|
|
"""Return a counts-weighted DRM.
|
|
For long spectral accumulations, the time series of responses needs
|
|
to be weighted by the counts history
|
|
|
|
Args:
|
|
weights (np.array): The weights of the counts history. Must sum to 1.0
|
|
bin_centers (np.array): The times at the center of the time bins
|
|
|
|
Returns:
|
|
np.array: The weighted DRM
|
|
"""
|
|
# only one DRM in the file; return it
|
|
if self.numdrms == 1:
|
|
return self._drm_list[0]
|
|
|
|
num_weights = len(weights)
|
|
if len(bin_centers) != num_weights:
|
|
raise ValueError("weights and bin_centers are not equal length")
|
|
# weight the matrices by counts
|
|
weightV = []
|
|
weightS = 0.0
|
|
weightD = []
|
|
currDRM = -1
|
|
for i in range(num_weights):
|
|
idx = np.argmin(np.abs(self.tstart - bin_centers[i]))
|
|
if (currDRM == -1) or (currDRM == idx):
|
|
weightS += weights[i]
|
|
else:
|
|
weightV.append(weightS)
|
|
weightS = weights[i]
|
|
weightD.append(idx)
|
|
currDRM = idx
|
|
|
|
weightV.append(weightS)
|
|
response = np.zeros((self.numebins, self.numchans))
|
|
for i, weight in zip(weightD, weightV):
|
|
response += self._drm_list[i] * weight
|
|
|
|
return response
|
|
|
|
def write_weighted(self, tstart, tstop, weights, directory, filename=None):
|
|
"""Writes a counts-weighted DRM to a FITS file.
|
|
|
|
Args:
|
|
tstart (np.array): The start times of the count history
|
|
tstop (np.array): The end times of the count history
|
|
weights (np.array): The normalized weights
|
|
directory (str): The directory to write the file
|
|
filename (str, optional): If set, will override the standardized name
|
|
"""
|
|
# set filename
|
|
if (self.filename is None) and (filename is None):
|
|
raise NameError('Filename not set')
|
|
if filename is None:
|
|
filename = self.filename
|
|
full_path = os.path.join(directory, filename)
|
|
|
|
# set dates/times
|
|
trigtime = self.trigtime if self.trigtime is not None else 0.0
|
|
startStr = Met(tstart[0] + trigtime).iso()
|
|
stopStr = Met(tstop[-1] + trigtime).iso()
|
|
todaysDate = Met.now().iso()
|
|
|
|
# adjust primary header
|
|
prihdr = self._headers['PRIMARY']
|
|
prihdr['FILENAME'] = filename
|
|
prihdr.set("DRM_NUM", 1)
|
|
prihdr['TSTART'] = tstart[0] + trigtime
|
|
prihdr['TSTOP'] = tstop[-1] + trigtime
|
|
prihdr['DATE-OBS'] = startStr
|
|
prihdr['DATE-END'] = stopStr
|
|
prihdr['DATE'] = todaysDate
|
|
prihdr.set("HISTORY",
|
|
"Created by writeResponse() in data.py from GSPEC")
|
|
prihdr.set("HISTORY", "from original file: %s" % self.filename)
|
|
prihdr.set("COMMENT",
|
|
"The matrix has been weighted by the count history")
|
|
|
|
# make ebounds table
|
|
ebounds = self.ebounds
|
|
|
|
# make weighted matrix
|
|
bin_centers = (tstart + tstop) / 2.0
|
|
matrix = self.weighted(weights, bin_centers)
|
|
|
|
drm = np.recarray(self.numebins, dtype=[('ENERG_LO', '>f4'),
|
|
('ENERG_HI', '>f4'),
|
|
('N_GRP', '>i2'),
|
|
('F_CHAN', '>i4'),
|
|
('N_CHAN', '>i4'),
|
|
('MATRIX', '>i4',
|
|
(self.numchans,))])
|
|
drm['ENERG_LO'], drm['ENERG_HI'] = self.photon_bins
|
|
drm['N_GRP'] = self._ngrp_list[0]
|
|
drm['F_CHAN'] = self._fchan_list[0]
|
|
drm['N_CHAN'] = self._nchan_list[0]
|
|
drm['MATRIX'] = matrix
|
|
|
|
drmhdr = self._headers['SPECRESP MATRIX1']
|
|
drmhdr['TSTART'] = tstart[0] + trigtime
|
|
drmhdr['TSTOP'] = tstop[-1] + trigtime
|
|
drmhdr['DATE-OBS'] = startStr
|
|
drmhdr['DATE-END'] = stopStr
|
|
drmhdr['DATE'] = todaysDate
|
|
|
|
# create FITS
|
|
hdulist = fits.HDUList()
|
|
primary_hdu = fits.PrimaryHDU(header=prihdr)
|
|
hdulist.append(primary_hdu)
|
|
|
|
ebounds_hdu = fits.BinTableHDU(data=ebounds, name='EBOUNDS',
|
|
header=self._headers['EBOUNDS'])
|
|
hdulist.append(ebounds_hdu)
|
|
|
|
drm_hdu = fits.BinTableHDU(data=drm, name='SPECRESP MATRIX',
|
|
header=drmhdr)
|
|
hdulist.append(drm_hdu)
|
|
|
|
hdulist.writeto(full_path, checksum=True)
|
|
|
|
def write(self, directory, filename=None, **kwargs):
|
|
"""Writes a single-DRM RSP object to a FITS file.
|
|
|
|
Args:
|
|
directory (str): Where to write the file
|
|
filename (str, optional): The filename to write to
|
|
"""
|
|
# set filename
|
|
if (self.filename is None) and (filename is None):
|
|
raise NameError('Filename not set')
|
|
if filename is None:
|
|
filename = self.filename
|
|
self._full_path = os.path.join(directory, filename)
|
|
self._headers['PRIMARY']['FILENAME'] = filename
|
|
|
|
hdulist = fits.HDUList()
|
|
primary_hdu = fits.PrimaryHDU(header=self.headers['PRIMARY'])
|
|
hdulist.append(primary_hdu)
|
|
|
|
ebounds_hdu = fits.BinTableHDU(data=self._ebounds, name='EBOUNDS',
|
|
header=self.headers['EBOUNDS'])
|
|
hdulist.append(ebounds_hdu)
|
|
|
|
compressed_drm = self._compress_drm(0)
|
|
|
|
col0 = fits.Column(name='ENERG_LO', format='F8',
|
|
array=self.photon_bins[0])
|
|
col1 = fits.Column(name='ENERG_HI', format='F8',
|
|
array=self.photon_bins[1])
|
|
col2 = fits.Column(name='N_GRP', format='I2', array=self._ngrp_list[0])
|
|
col3 = fits.Column(name='F_CHAN', format='PI()',
|
|
array=self._fchan_list[0])
|
|
col4 = fits.Column(name='N_CHAN', format='PI()',
|
|
array=self._nchan_list[0])
|
|
col5 = fits.Column(name='MATRIX', format='PE()', array=compressed_drm)
|
|
|
|
drm_hdu = fits.BinTableHDU.from_columns(
|
|
[col0, col1, col2, col3, col4, col5],
|
|
name='SPECRESP MATRIX',
|
|
header=self.headers['SPECRESP MATRIX'])
|
|
hdulist.append(drm_hdu)
|
|
|
|
hdulist.writeto(self.full_path, checksum=True, **kwargs)
|
|
|
|
def resample(self, index=0, num_photon_bins=None, photon_bin_edges=None,
|
|
num_interp_points=20, interp_kind='linear'):
|
|
"""Resamples the incident photon axis of a DRM and returns a new
|
|
response object. Resampling can be used to downgrade the photon energy
|
|
resolution, upgrade the resolution, and/or redefine the edges of the
|
|
incident photon bins. By definition, the resampling can only be
|
|
performed within the photon energy range of the current object.
|
|
|
|
The process for resampling is to create a high-resolution grid for each
|
|
new photon bin, interpolate the differential effective area onto that
|
|
grid (interpolation is performed on log10(effective area)), and then
|
|
integrate the differential effective area over the grid points in each
|
|
bin. The final effective area is then scaled by the ratio of the
|
|
initial number of photon bins to the requested number of photon bins.
|
|
|
|
Args:
|
|
index (int, optional): The DRM index. Default is 0
|
|
num_photon_bins (int, optional): The number of photon bins in the
|
|
new DRM. The bin edges will be
|
|
generated logarithmically. Only set
|
|
this or `photon_bin_edges`.
|
|
photon_bin_edges (np.array, optional): The array of photon bin edges.
|
|
Only set this or
|
|
`num_photon_bins`
|
|
num_interp_points (int, optional): The number of interpolation points
|
|
used to integrate over for each
|
|
new photon bin. Default is 20.
|
|
interp_kind (str, optional): The kind of interpolation to be
|
|
passed to scipy.interp1d. Default is
|
|
'linear'.
|
|
|
|
Returns:
|
|
(:class:`RSP`): A new resampled single-DRM response object
|
|
|
|
"""
|
|
if (num_photon_bins is None) and (photon_bin_edges is None):
|
|
raise ValueError('Either num_photon_bins or photon_bin_edges must' \
|
|
' be set')
|
|
elif num_photon_bins is not None:
|
|
try:
|
|
num_photon_bins = int(num_photon_bins)
|
|
except:
|
|
raise TypeError('num_photon_bins must be a positive integer')
|
|
if num_photon_bins < 1:
|
|
raise ValueError('num_photon_bins must be a positive integer')
|
|
orig_edges = self.photon_bins
|
|
photon_bin_edges = np.logspace(*np.log10((orig_edges[0][0],
|
|
orig_edges[1][-1])),
|
|
num_photon_bins+1)
|
|
elif photon_bin_edges is not None:
|
|
if not isinstance(photon_bin_edges, np.ndarray):
|
|
raise TypeError('photon_bin_edges must be a numpy array')
|
|
orig_edges = self.photon_bins
|
|
badmask = (photon_bin_edges < orig_edges[0][0]) | \
|
|
(photon_bin_edges > orig_edges[1][-1])
|
|
if badmask.sum() > 0:
|
|
raise ValueError('photon_bin_edges are beyond valid range of '\
|
|
'{0:.2f}-{1:.2f}'.format(orig_edges[0][0],
|
|
orig_edges[1][-1]))
|
|
photon_bin_edges = np.sort(photon_bin_edges)
|
|
num_photon_bins = photon_bin_edges.size-1
|
|
|
|
# differential effective area
|
|
diff_effarea = self.drm(index)/self.photon_bin_widths[:,np.newaxis]
|
|
|
|
# create an interpolation array over each new bin
|
|
photon_bounds = np.vstack((photon_bin_edges[:-1], photon_bin_edges[1:]))
|
|
interp_arr = np.logspace(*np.log10(photon_bounds), num_interp_points+2,
|
|
axis=0)
|
|
|
|
# create the DRM interpolator
|
|
x = np.append(self.photon_bins[0], self.photon_bins[1][-1])
|
|
y = np.vstack((diff_effarea, diff_effarea[-1,:]))
|
|
y[y == 0] = 1e-10
|
|
effarea_interpolator = interp1d(x, np.log10(y), axis=0,
|
|
kind=interp_kind,
|
|
fill_value='extrapolate')
|
|
|
|
# interpolate the DRM over the photon interpolation array
|
|
diff_effarea_interp = 10.**effarea_interpolator(interp_arr)
|
|
badmask = np.isnan(diff_effarea_interp)
|
|
diff_effarea_interp[badmask] = 0.0
|
|
|
|
# integrate the DRM over the photon interpolation array
|
|
diff_effarea_new = trapz(diff_effarea_interp,
|
|
interp_arr[:,:,np.newaxis], axis=0)
|
|
|
|
# scale by ratio of photon bins
|
|
drm = diff_effarea_new / (self.numebins/num_photon_bins)
|
|
|
|
# create response object
|
|
drm_ext = 'SPECRESP MATRIX'
|
|
if drm_ext not in self.headers.keys():
|
|
drm_ext = 'SPECRESP MATRIX{}'.format(index + 1)
|
|
atscat = None
|
|
if 'INFILE04' in self.headers[drm_ext]:
|
|
atscat = self.headers[drm_ext]['INFILE04']
|
|
obj = RSP.from_arrays(self._ebounds['E_MIN'], self._ebounds['E_MAX'],
|
|
*photon_bounds, drm, trigtime=self.trigtime,
|
|
tstart=self.tstart[index] + self.trigtime,
|
|
tstop=self.tstop[index] + self.trigtime,
|
|
detnam=self.detector,
|
|
object=self.headers['PRIMARY']['OBJECT'],
|
|
ra_obj=self.headers['PRIMARY']['RA_OBJ'],
|
|
dec_obj=self.headers['PRIMARY']['DEC_OBJ'],
|
|
infile=self.filename,
|
|
ch2e_ver=self.headers['EBOUNDS']['CH2E_VER'],
|
|
gain=self.headers['EBOUNDS']['GAIN_COR'],
|
|
mat_type=self.headers[drm_ext]['MAT_TYPE'],
|
|
src_az=self.headers[drm_ext]['SRC_AZ'],
|
|
src_el=self.headers[drm_ext]['SRC_EL'],
|
|
geo_az=self.headers[drm_ext]['GEO_AZ'],
|
|
geo_el=self.headers[drm_ext]['GEO_EL'],
|
|
det_ang=self.headers[drm_ext]['DET_ANG'],
|
|
geo_ang=self.headers[drm_ext]['GEO_ANG'],
|
|
atscat=atscat, datatype=self.datatype)
|
|
|
|
obj._ngrp_list.append(self._ngrp_list[index])
|
|
obj._fchan_list.append(self._fchan_list[index])
|
|
obj._nchan_list.append(self._nchan_list[index])
|
|
return obj
|
|
|
|
def rebin(self, index=0, factor=None, edge_indices=None):
|
|
"""Rebins the channel energy axis of a DRM and returns a new response
|
|
object. Rebinning can only be used to downgrade the channel resolution
|
|
and is constrained to the channel edges of the current DRM.
|
|
|
|
Rebinning can be performed by either defining an integral factor of the
|
|
number of current energy channels to combine (e.g. factor=4 for 128
|
|
energy channels would result in 32 energy channels), or by defining an
|
|
index array into the current channel edges that will define the new set
|
|
of edges.
|
|
|
|
Args:
|
|
index (int, optional): The DRM index. Default is 0
|
|
factor (int, optional): The rebinning factor. Must set either this
|
|
or `edge_indices`
|
|
edge_indices (np.array, optional): The index array that represents
|
|
which energy edges should remain
|
|
in the rebinned DRM.
|
|
|
|
Returns:
|
|
(:class:`RSP`): A new rebinned single-DRM response object
|
|
"""
|
|
if (factor is None) and (edge_indices is None):
|
|
raise ValueError('Either factor or edge_indices must be set')
|
|
elif factor is not None:
|
|
try:
|
|
factor = int(factor)
|
|
except:
|
|
raise TypeError('factor must be a positive integer')
|
|
if factor < 1:
|
|
raise ValueError('factor must be a positive integer')
|
|
if (self.numchans % factor) > 0:
|
|
raise ValueError('factor must be factor of numchans: '\
|
|
'{}'.format(self.numchans))
|
|
emin = self.ebounds['E_MIN'].reshape(-1, factor)[:,0]
|
|
emax = self.ebounds['E_MAX'].reshape(-1, factor)[:,-1]
|
|
elif edge_indices is not None:
|
|
if isinstance(edge_indices, (list, tuple)):
|
|
edge_indices = np.array(edge_indices)
|
|
if not isinstance(edge_indices, np.ndarray):
|
|
raise TypeError('new_edges must be a list, tuple or array')
|
|
edge_indices = np.unique(edge_indices.astype(int))
|
|
if (edge_indices[0] < 0) or (edge_indices[-1] > self.numchans-1):
|
|
raise ValueError('edge_indices outside valid range of ' \
|
|
'0-{}'.format(self.numchans-1))
|
|
|
|
|
|
old_edges = np.append(self.ebounds['E_MIN'],
|
|
self.ebounds['E_MAX'][-1])
|
|
new_edges = old_edges[edge_indices]
|
|
emin = new_edges[:-1]
|
|
emax = new_edges[1:]
|
|
|
|
|
|
effarea = self.drm(index)
|
|
# the new effective area matrix will just be summed over the channel
|
|
# axis for the bins we are combining
|
|
new_effarea = np.zeros((self.numebins, emin.size))
|
|
old_emin, old_emax = self.ebounds['E_MIN'], self.ebounds['E_MAX']
|
|
sidx = (old_emin[:,np.newaxis] == emin[np.newaxis,:]).nonzero()[0]
|
|
eidx = (old_emax[:,np.newaxis] == emax[np.newaxis,:]).nonzero()[0]
|
|
for i in range(emin.size):
|
|
new_effarea[:,i] = effarea[:,sidx[i]:eidx[i]+1].sum(axis=1)
|
|
|
|
# create the new object
|
|
drm_ext = 'SPECRESP MATRIX'
|
|
if drm_ext not in self.headers.keys():
|
|
drm_ext = 'SPECRESP MATRIX{}'.format(index + 1)
|
|
atscat = None
|
|
if 'INFILE04' in self.headers[drm_ext]:
|
|
atscat = self.headers[drm_ext]['INFILE04']
|
|
obj = RSP.from_arrays(emin, emax, *self._pbounds, new_effarea,
|
|
trigtime=self.trigtime,
|
|
tstart=self.tstart[index] + self.trigtime,
|
|
tstop=self.tstop[index] + self.trigtime,
|
|
detnam=self.detector,
|
|
object=self.headers['PRIMARY']['OBJECT'],
|
|
ra_obj=self.headers['PRIMARY']['RA_OBJ'],
|
|
dec_obj=self.headers['PRIMARY']['DEC_OBJ'],
|
|
infile=self.filename,
|
|
ch2e_ver=self.headers['EBOUNDS']['CH2E_VER'],
|
|
gain=self.headers['EBOUNDS']['GAIN_COR'],
|
|
mat_type=self.headers[drm_ext]['MAT_TYPE'],
|
|
src_az=self.headers[drm_ext]['SRC_AZ'],
|
|
src_el=self.headers[drm_ext]['SRC_EL'],
|
|
geo_az=self.headers[drm_ext]['GEO_AZ'],
|
|
geo_el=self.headers[drm_ext]['GEO_EL'],
|
|
det_ang=self.headers[drm_ext]['DET_ANG'],
|
|
geo_ang=self.headers[drm_ext]['GEO_ANG'],
|
|
atscat=atscat, datatype=self.datatype)
|
|
|
|
obj._ngrp_list.append(self._ngrp_list[index])
|
|
obj._fchan_list.append(self._fchan_list[index])
|
|
obj._nchan_list.append(self._nchan_list[index])
|
|
return obj
|
|
|
|
def _decompress_drm(self, drm_data, spec_index):
|
|
"""The GBM DRM is not stored as 2D matrix and may be compressed, so this
|
|
constructs and decompresses the DRM.
|
|
|
|
Args:
|
|
drm_data (np.recarray): The DRM data
|
|
spec_index (int): The index of the DRM in the event of multiple DRMs
|
|
|
|
Returns:
|
|
np.array: 2D matrix containing the unpacked DRM
|
|
"""
|
|
num_photon_bins = self._headers['SPECRESP MATRIX' + str(spec_index)][
|
|
'NAXIS2']
|
|
first_chans = drm_data['F_CHAN']
|
|
num_chans = drm_data['N_CHAN']
|
|
|
|
matrix = drm_data['MATRIX']
|
|
energy_channel_max = self._ebounds['E_MAX']
|
|
num_channels = energy_channel_max.size
|
|
|
|
# The format of the compress matrix is a series of groups, for each
|
|
# energy bin, of channels with non-zero values.
|
|
# fchan stands for the first channel of each of these groups
|
|
# and nchan for the number of channels in the group group.
|
|
# Each row in the matrix is a 1D list consisting on the contatenated
|
|
# values of all groups for a given energy bin
|
|
# Note that in FITS the first index is 1
|
|
drm = np.zeros((num_photon_bins, num_channels))
|
|
for fchans,nchans,effective_area,drm_row \
|
|
in zip(first_chans,num_chans,matrix,drm):
|
|
|
|
channel_offset = 0
|
|
|
|
for fchan,nchan in zip(fchans,nchans):
|
|
|
|
start_idx = fchan - 1
|
|
end_idx = start_idx + nchan
|
|
|
|
drm_row[start_idx:end_idx] = \
|
|
effective_area[channel_offset:channel_offset+nchan]
|
|
|
|
channel_offset += nchan
|
|
|
|
return drm
|
|
|
|
def _compress_drm(self, spec_index):
|
|
"""Compress a DRM by removing all the zeros. This can result in a
|
|
file that is ~50% the size of an uncompressed DRM because the DRM is
|
|
largely a triangle matrix.
|
|
|
|
Also modifies the helper FITS columns that are used to decompress the
|
|
matrix
|
|
|
|
Args:
|
|
spec_index (int): The index of the DRM in the event of multiple DRMs
|
|
|
|
Returns:
|
|
np.array: An array of variable length arrays containing the compressed matrix
|
|
"""
|
|
|
|
# function to split channels into contiguous groups
|
|
def group_consecutive_indices(indices):
|
|
# change points where consecutive indices > 1
|
|
diff = indices[1:] - indices[0:-1]
|
|
diff_idx = np.where(diff > 1)[0] + 1
|
|
# add first, last indices
|
|
diff_idx = np.concatenate(([0], diff_idx, [idx.size]))
|
|
# group into consecutive indices
|
|
groups = [idx[diff_idx[i]:diff_idx[i + 1]] \
|
|
for i in range(diff_idx.size - 1)]
|
|
return groups
|
|
|
|
drm = self._drm_list[spec_index]
|
|
|
|
matrix = np.zeros((self.numebins,), dtype=np.object_)
|
|
first_chans = []
|
|
num_chans = []
|
|
num_groups = []
|
|
for ibin in range(self.numebins):
|
|
idx = np.where(drm[ibin, :] > 0.0)[0]
|
|
# if all zeros
|
|
if idx.size == 0:
|
|
num_groups.append(1)
|
|
first_chans.append(np.array([self.numchans]))
|
|
num_chans.append(np.array([1]))
|
|
matrix[ibin] = np.array([0.0])
|
|
else:
|
|
groups = group_consecutive_indices(idx)
|
|
num_groups.append(len(groups))
|
|
first_chans.append(
|
|
np.array([group[0] + 1 for group in groups]))
|
|
num_chans.append(np.array([len(group) for group in groups]))
|
|
matrix[ibin] = np.concatenate(
|
|
[drm[ibin, group] for group in groups])
|
|
|
|
self._ngrp_list[spec_index] = np.array(num_groups)
|
|
self._fchan_list[spec_index] = first_chans
|
|
self._nchan_list[spec_index] = num_chans
|
|
|
|
return matrix
|
|
|
|
def _assert_range(self, valrange):
|
|
assert valrange[0] <= valrange[1], \
|
|
'Range must be in increasing order: (lo, hi)'
|
|
return valrange
|