# coding: utf-8
""" Manipulate python data to produce Xspec readable files"""
from astropy.io import fits
import numpy as np
import os
[docs]def read_from_txt(filename):
"""Reads a text file to create xspec readable fits file
Parameters
----------
filename: str
A text files with 2 or 3 columns to read.
Returns
-------
energies, values(, erros): ndarrays
The columns present in the file text
Raises
------
IOError
If the files has less than 2 or more than 3 columns or is read as such
due to formatting.
"""
table = np.loadtxt(filename).T
if table.shape[0] == 2:
energies = table[0]
yvalues = table[1]
return energies, yvalues
elif table.shape[0] == 3:
energies = table[0]
yvalues = table[1]
yerrors = table[2]
return energies, yvalues, yerrors
else:
raise IOError("File has {} columns".format(table.shape[0]))
[docs]def fake_pha(values, energies,
errors=None,
prefix="fake",
respfile="NONE",
telescope="FAKE", instrument="FAKE",
afilter="NONE", backfile="NONE",
corrfile="NONE", ancrfile="NONE",
syserr=0, date="2001-01-01"):
""" Creates a fake pha file to be read fitted with xspec
Parameters
----------
values: ndarray
Values of the fake data
energies: ndarray
Array of bin edges for the x-axis (lenght of values + 1)
errors: ndarray, optional
statistical error of the fake data. Default is zero
prefix: str, optional
A prefix to save file as <prefix>.pha. Defauls is "fake"
Other Parameters
----------------
date: str
File creation date. Default is "2001-01-01"
telescope: str
Name of the mission/telescope
instrument: str
Name of the instrument used
afilter: str
FILTER keyword
respfile: str
Response file.
backfile: str
Background file
corrfile: str
CORRFILE keyword
ancrfile: str
ANCRFILE keyword
syserr: float
systematic error on data
Returns
-------
Creates <prefix>.pha file in the current directory
"""
hdu1 = fits.PrimaryHDU()
hdu1.header["DATE"] = date, "Fake creation date"
hdu1.header["TELESCOP"] = telescope, "Fake mission name"
hdu1.header["INSTRUME"] = instrument, "Fake instrument name"
hdu1.header["CONTENT"] = "SPECTRUM"
hdu1.header["PHAVERSN"] = "1992a"
if errors is None:
errors = np.zeros_like(values)
counts = np.array([val*(energies[i+1]-energies[i]) for i,val in enumerate(values)])
stat_err = np.array([err*(energies[i+1]-energies[i]) for i,err in enumerate(errors)])
channels = np.arange(1, len(energies))
col_channel = fits.Column(name="CHANNEL", format="1I", array=channels)
col_counts = fits.Column(name="COUNTS", format="1D", array=counts)
col_stat_err = fits.Column(name="STAT_ERR", format="1D", array=stat_err)
pha_cols = fits.ColDefs([col_channel, col_counts, col_stat_err])
pha_head = fits.Header()
pha_head["EXTNAME"] = "SPECTRUM"
pha_head["TELESCOP"] = telescope
pha_head["INSTRUME"] = instrument
pha_head["FILTER"] = afilter
pha_head["EXPOSURE"] = 1.0
pha_head["AREASCAL"] = 1.0
pha_head["BACKSCAL"] = 1.0
pha_head["CORRSCAL"] = 1.0
pha_head["BACKFILE"] = backfile
pha_head["CORRFILE"] = corrfile
pha_head["RESPFILE"] = respfile
pha_head["ANCRFILE"] = ancrfile
pha_head["POISSERR"] = False
pha_head["CHANTYPE"] = "PHA"
pha_head["DETCHANS"] = len(channels)
pha_head["SYS_ERR"] = syserr
pha_head["QUALITY"] = 0
pha_head["GROUPING"] = 0
pha_head["HDUCLASS"] = "OGIP"
pha_head["HDUCLAS1"] = "SPECTRUM"
pha_head["HDUVERS"] = "1.1.0"
spechdu = fits.BinTableHDU.from_columns(pha_cols, header=pha_head)
phafile = fits.HDUList([hdu1, spechdu])
phafile.writeto("{}.pha".format(prefix), overwrite=True)
[docs]def fake_rmf(energies, prefix="fake",
telescope="FAKE", instrument="FAKE",
date="2001-01-01", afilter="NONE",
effarea=1.0):
""" Creates a dummy rmf file with a diagonal response
Parameters
----------
energies: ndarray
Array of bin edges for the x-axis (lenght of values + 1)
prefix: str, optional
A prefix to save file as <prefix>.pha. Defauls is "fake"
Other Parameters
----------------
date: str
File creation date. Default is "2001-01-01"
telescope: str
Name of the mission/telescope
instrument: str
Name of the instrument used
afilter: str
FILTER keyword
effarea: float
EFFAREA keyword.
Returns
-------
Creates <prefix>.rmf file in the current directory
"""
prihdu = fits.PrimaryHDU()
prihdu.header["DATE"] = date, "Fake creation date"
prihdu.header["TELESCOP"] = telescope, "Fake mission name"
prihdu.header["INSTRUME"] = instrument, "Fake instrument name"
prihdu.header["CONTENT"] = "SPECTRUM"
channels = np.arange(len(energies) - 1)
energ_lo = np.array(energies[:-1])
energ_hi = np.array(energies[1:])
n_grp = np.ones_like(channels)
f_chan = np.ones_like(channels)
n_chan = np.full_like(channels, len(channels))
matrix = []
identi = np.identity(len(channels))
for row in identi:
matrix.append(np.array(identi))
matrix = np.asarray(matrix)
col_energ_lo = fits.Column(name="ENERG_LO", format="1D",
array=energ_lo, unit="keV")
col_energ_hi = fits.Column(name="ENERG_HI", format="1D",
array=energ_hi, unit="keV")
col_n_grp = fits.Column(name="N_GRP", format="1I", array=n_grp)
col_f_chan = fits.Column(name="F_CHAN", format="1I", array=f_chan)
col_n_chan = fits.Column(name="N_CHAN", format="1I", array=n_chan)
col_matrix = fits.Column(name="MATRIX", format="{}E".format(len(channels)),
array=identi)
matrix_cols = fits.ColDefs([col_energ_lo,
col_energ_hi,
col_n_grp,
col_f_chan,
col_n_chan,
col_matrix])
m_head = fits.Header()
m_head["EXTNAME"] = "MATRIX", "name of this binary table extension"
m_head["TELESCOP"] = telescope, "Fake Mission name"
m_head["INSTRUME"] = instrument, "Fake instrument name"
m_head["FILTER"] = afilter
m_head["CHANTYPE"] = "PHA"
m_head["DETCHANS"] = len(channels)
m_head["LO_THRES"] = 0
m_head["EFFAREA"] = effarea
m_head["RMFVERSN"] = "1992a"
m_head["DATE"] = date, "Fake production date"
matrixhdu = fits.BinTableHDU.from_columns(matrix_cols, header=m_head)
col_emin = fits.Column(name="E_MIN", format="1D", array=energ_lo,
unit="keV")
col_emax = fits.Column(name="E_MAX", format="1D", array=energ_hi,
unit="keV")
col_channel = fits.Column(name="CHANNEL", format="1J", array=channels)
bound_cols = fits.ColDefs([col_emin,
col_emax,
col_channel])
head2 = fits.Header()
head2["EXTNAME"] = "EBOUNDS", "name of this binary table extension"
head2["FILTER"] = afilter
head2["CHANTYPE"] = "PHA"
head2["EFFAREA"] = effarea
head2["RMFVERSN"] = "1992a"
head2["DATE"] = date, "Fake production date"
boundshdu = fits.BinTableHDU.from_columns(bound_cols, header=head2)
rmffile = fits.HDUList([prihdu, matrixhdu, boundshdu])
rmffile.writeto("{}.rmf".format(prefix), overwrite=True)