#!/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)