#!/usr/bin/env python
#
"""
This file is part of ESA's XMM-Newton Scientific Analysis System (SAS).
SAS 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.
SAS 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 SAS. If not, see .
Name = 'spec2typeii'
VERSION = '1.5';
author = 'Jose Vicente Perea';
date = '2021-02-19';
ChangeLog
=========
Version 1.5 - 2021-02-19 (JVP)
------------
+ Version output string fixed
Version 1.4 - 2020-12-21 (JVP)
------------
+ Column value to NaN if no EXPOSURE keyword in spec file.
evselect produces empty table and no EXPOSURE keyword FITS file
if all the events are filtered out.
Version 1.3 - 2020-10-16 (JVP)
------------
+ Minor change about FITS columns dimensions
Some cleaning of commented lines
Version 1.2 - 2020-10-13 (JVP)
------------
+ One more spectrum included in the merged combined spectra file (OGIP TYPE II FITS file):
Source extraction region: Circle, but only Doubles (PATTERN == [1:4])
(Input patternD doubles spectrum file)
Version 1.1 - 2020-09-14 (JVP)
------------
+ BACKFILE, RESPFILE and ANCRFILE good filenames included
in the proper columns of the final multispectra file
Version 1.0 - 2020-08-13 (JVP)
------------
+ Initial version
Task to merge a set of spectre
and produce a final OGIP TypeII output
"""
__version__ = "1.5"
__author__ = "Jose V. Perea"
__copyright__ = "ESA (C) 2000-2020"
__license__ = "GPL"
version_num = __version__
# -------------------
## For draft version
#files = [
#'spec_pattern0.fits',
#'spec_ann_1.fits',
#'spec_ann_2.fits',
#'spec_ann_3.fits',
#'spec_ann_4.fits',
#'spec_ann_5.fits',
#'spec_ann_6.fits',
#'spec_ann_7.fits',
#'spec_ann_8.fits'
#]
# -------------------
########################################################################
#
# outspec,outbkgd,outarf
#
#import os
#pwd = os.getcwd()
#if pwd == '/home/xmm/jperea/Utils/multispec_annulus':
# lab = 'pere'
#elif pwd == '/home/xmm/jperea/Utils/multispec_annulus/from_PR':
# lab = 'PR'
#
## outspec = 'spec_tests_' + lab + '.fits'
## outbkgd = 'back_tests_' + lab + '.fits'
## outarf = 'arfs_tests_' + lab + '.fits'
########################################################################
###############################
# Input files
#
import sys, os, argparse
parser = argparse.ArgumentParser(description='Merge spectra to a single OGIP TypeII FITS output',
formatter_class=argparse.ArgumentDefaultsHelpFormatter)
parser.add_argument('--version', '-v', action='version', version=f'%(prog)s (%(prog)s-{version_num}) [-]')
parser.add_argument('--specPattern0', help='Input pattern0 singles spectrum file', default='spec_pattern0.fits')
parser.add_argument('--specPatternD', help='Input patternD doubles spectrum file', default='spec_patternD.fits')
parser.add_argument('--spectraAnnulus', help='Input spec from Annulus files', default='spec_annulus_001.fits,spec_annulus_002.fits')
parser.add_argument('--outspec', help='Output source combined spectrum file', default='outspec.fits')
parser.add_argument('--outbkgd', help='Output background spectrum file', default='outbkgd.fits')
parser.add_argument('--outarf', help='Output arf file', default='outarf.fits')
args = parser.parse_args()
# files = [ args.specPattern0 ]
files = [ args.specPattern0, args.specPatternD ]
for i in sorted(args.spectraAnnulus.split(',')):
files.append(i)
for f in files:
if not os.path.isfile(f):
print(f'Some of the input files does not exist: {f}. Exiting')
sys.exit(1)
outspec = args.outspec
outbkgd = args.outbkgd
outarf = args.outarf
###############################
from astropy.io import fits
from astropy.table import Table
nf = len(files)
print(files)
print(nf)
with fits.open(files[0]) as hdul:
hin = hdul['SPECTRUM'].header
nchan = hin['NAXIS2']
nebin = fits.getval( hin['ANCRFILE'], 'NAXIS2', 'SPECRESP' )
#print(f'SPECTRUM NAXIS2 = {nchan}')
#print(f'ANCRFILE NAXIS2 = {nebin}')
regs = []
out = []
outb = []
outa = []
for i in range(0, nf):
print(f'\nOpen file : {files[i]}')
with fits.open( files[i] ) as fhdul:
reg_ext = [ s.name for s in fhdul if s.name.startswith('REG') ][0]
regs.append( {'data': fhdul[reg_ext].data , 'hdr': fhdul[reg_ext].header} )
a = fhdul['SPECTRUM'].data
ha = fhdul['SPECTRUM'].header
#print( len(a) )
bfile = ha['BACKFILE']
ancrfile = ha['ANCRFILE']
respfile = ha['RESPFILE']
# == PENDING == INCLUDE fileexists check
print(f'\tOpen backfile : {bfile}')
with fits.open( bfile ) as bhdul:
b = bhdul['SPECTRUM'].data
hb = bhdul['SPECTRUM'].header
# == PENDING == INCLUDE fileexists check
print(f'\tOpen ancrfile : {ancrfile}')
with fits.open( ancrfile ) as ahdul:
arf = ahdul['SPECRESP'].data
harf = ahdul['SPECRESP'].header
if i == 0:
hout = ha
hout['HDUCLAS4'] = 'TYPE:II'
hout['HDUVERS1'] = '1.2.1'
hout['HDUVERS'] = '1.2.1'
hout.pop('BACKFILE') # remove keywords from final FITS file.
hout.pop('RESPFILE') # BACKFILE, RESPFILE, ANCRFILE will be column names
hout.pop('ANCRFILE')
houtb = hb
houtb['HDUCLAS4'] = 'TYPE:II'
houtb['HDUVERS1'] = '1.2.1'
houtb['HDUVERS'] = '1.2.1'
houtarf = harf
houtarf['HDUCLAS4'] = 'TYPE:II'
outarfFileName = os.path.basename(outarf).replace('FIT', 'FTZ')
backfileFileName = os.path.basename(bfile).replace('FIT', 'FTZ')
respfileFileName = os.path.basename(respfile).replace('FIT', 'FTZ')
if 'EXPOSURE' in ha:
ha_exposure = ha['EXPOSURE']
else:
import numpy as np
ha_exposure = np.nan
if 'EXPOSURE' in hb:
hb_exposure = hb['EXPOSURE']
else:
import numpy as np
hb_exposure = np.nan
out.append( {'spec_num': i+1,
'channel' : a.CHANNEL,
'counts' : a.COUNTS,
'quality' : a.QUALITY,
'grouping': a.GROUPING,
'rowid' : str(i+1),
'exposure': ha_exposure,
'backfile': backfileFileName,
'backsc' : ha['BACKSCAL'],
'respfile': respfileFileName,
'ancrfile': f'{outarfFileName}{{{i+1}}}' # ha['ANCRFILE'] or outarf ???
} )
outb.append( {'spec_num': i+1,
'channel' : b.CHANNEL,
'counts' : b.COUNTS,
'rowid' : str(i+1),
'exposure': hb_exposure,
'backsc' : hb['BACKSCAL']
} )
outa.append( {'arf_num' : i+1,
'energ_lo' : arf.ENERG_LO,
'energ_hi' : arf.ENERG_HI,
'specresp' : arf.SPECRESP,
'rowid' : str(i+1)
} )
#########################################
hdul = []
h0 = fits.getheader( files[0], 0 )
hdu0 = fits.PrimaryHDU( header=h0 )
hdul.append(hdu0)
# SPECTRUM --------------------------------------------------------------------------------
#
formats = {
'spec_num' : 'I'
, 'channel' : str(len(a)) + 'J'
, 'counts' : str(len(a)) + 'J'
, 'quality' : str(len(a)) + 'I'
, 'grouping': str(len(a)) + 'I'
, 'rowid' : '8A'
, 'exposure': 'E'
, 'backfile': '32A'
, 'respfile': '32A'
, 'ancrfile': '36A'
, 'backsc' : 'E'
}
formatsb = {
'spec_num' : 'I'
, 'channel' : str(len(b)) + 'J'
, 'counts' : str(len(b)) + 'J'
, 'rowid' : '8A'
, 'exposure': 'E'
, 'backsc' : 'E'
}
formatsarf = {
'arf_num' : 'I'
, 'energ_lo' : str(len(arf)) + 'E'
, 'energ_hi' : str(len(arf)) + 'E'
, 'specresp' : str(len(arf)) + 'E'
, 'rowid' : '8A'
}
fitsColumnDefinitions = []
for col in out[0].keys():
#
arr = [ s[col] for s in out ]
#print (f'COL : {col.upper()} {arr}')
#multi = ['channel', 'counts', 'quality', 'grouping']
#inArr = []
#if col not in multi:
# inArr = arr
#else:
# inArr = [ k for k in arr ]
inArr = [ k for k in arr ]
unit = None
if col == 'counts': unit='count'
fitsColumnDefinitions.append( fits.Column( name=col.upper(), format=formats[col], unit=unit, array= inArr ))
cols = fits.ColDefs( fitsColumnDefinitions )
hduspec = fits.BinTableHDU.from_columns( cols, header = hout )
hdul.append(hduspec)
# REGIONS --------------------------------------------------------------------------------
#
for i in range(0, len(regs)):
hdr = regs[i]['hdr']
hdr['EXTNAME'] = f'REG{i+1:05d}'
hdul.append( fits.BinTableHDU( data=regs[i]['data'], header=hdr) )
# ALL TOGETHER --------------------------------------------------------------------------------
#
#hdu1 = fits.BinTableHDU(header=hout)
#hdul = fits.HDUList([hdu0, hduspec])
hdul = fits.HDUList( hdul )
hdul.writeto(outspec, overwrite = True )
def writeFITS( h0=None, h1=None, out=None, formats=None, outFile=None ):
hdul = []
hdul.append( fits.PrimaryHDU(header=h0) )
fitsColumnDefinitions = []
for col in out[0].keys():
arr = [ s[col] for s in out ]
inArr = [ k for k in arr ]
unit = None
if col == 'counts': unit='count'
if col == 'energ_lo' or col == 'energ_hi': unit='keV'
if col == 'specresp': unit='cm2'
fitsColumnDefinitions.append( fits.Column( name=col.upper(), format=formats[col], unit=unit, array= inArr ))
cols = fits.ColDefs( fitsColumnDefinitions )
hdul.append( fits.BinTableHDU.from_columns( cols, header = h1 ) )
hdul = fits.HDUList( hdul )
hdul.writeto(outFile, overwrite = True )
# -----------------------------------------------------
# BACKFILE
#
hb0 = fits.getheader(bfile, 0)
writeFITS( h0=hb0, h1=houtb, out=outb, formats=formatsb, outFile=outbkgd )
# -----------------------------------------------------
# ANCRFILE
#
h_0 = fits.getheader(ancrfile, 0)
writeFITS( h0=h_0, h1=houtarf, out=outa, formats=formatsarf, outFile=outarf )
# --------------------------------------------------------------------------------------------------------------
sys.exit(0)