# 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 .
#
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
`_.
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 ():
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