#!/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 = 'srcTimeSeriesEnergyBands'
VERSION = '1.4';
author = 'Jose Vicente Perea';
date = '2018-06-07';
ChangeLog
=========
Version 1.4 - 2018-06-07 (JVP)
------------
+ Include Uncorrected Light Curves only if the corresponding Corrected ones exist.
We could have the uncorrected LCs per E band for a Corrected LC not selected for the final products (ex.: source spectrum with not enough number of counts).
Version 1.3 - 2018-05-16 (JVP)
------------
+ Important change to ensure good merging of tables. Tables join option : join_type='right'
Version 1.2 - 2018-05-08 (JVP)
------------
+ Final FITS file : write RATE table and modified header in different 'open file' instances to avoid writing conflicts
Version 1.1 - 2018-05-04 (JVP)
------------
+ Light Curves per Energy Band and per Exposure a intermediate/'srcTimeSeriesEnergyBands-correctedLightCurve-' + + '-plus_Rates_per_Energy_band.ftz'
Some verbose messages included
Module searching order removed : sys.path.insert()
Version 1.0 - 2018-04-24 (JVP)
------------
+ Initial version
Task for Light Curves merging per Energy band
"""
__version__ = "1.4"
__author__ = "Jose V. Perea"
__copyright__ = "ESA (C) 2000-2018"
__license__ = "GPL"
#print( __doc__ )
#------------------------------------------------------------------------------------------------
def usage():
#------------------------------------------------------------------------------------------------
'''
Help.
'''
print("""
usage: srcTimeSeriesEnergyBands.py --allEnergyRangeCorrLC=<...> --SrcLCsperEnergyBand=<...> --BkgLCsperEnergyBand=<...> [options]
output: ...
INTRODUCTION
The srcTimeSeriesEnergyBands task reads all the EPIC whole Energy band corrected (epiclccorr output) Light Curve FITS files (all exposures),
then all uncorrected Light Curves per Energy band. Uncorrected LCs are corrected using the 'BKGRATIO' value from the corrected ones.
All the resulting corrected LCs (whole E range LC + LCs per E band) are merged in the same FITS file table.
A mereged EPIC LCs plot is also created.
OPTIONS:
Input:
--allEnergyRangeCorrLC= Input EPIC corrected Light Curves
Output from 'epiclccorr'
--SrcLCsperEnergyBand Uncorrected Source Light Curves.
All Energy bands
--BkgLCsperEnergyBand Uncorrected Background Light Curves.
All Energy bands
Output:
-o, --outfits= Output FITS file.
All Light Curves included ( whole E range LC + LCs per E band )
--plotting= Plot all the merged Light Curves (Yes/No, Y/N). Default: Yes
--plotfile= Plot file name. Default: figure.png
--plotformat={png,pdf} Output figure format. Default: png
--interactplot Output plot in an interactive graphics window (Y) or in a PNG file (N)
(not implemented yet)
--verbose= Increase verbosity (Yes/No, Y/N)
-v, --version Show version message
Help:
-h, --help Print this help.
DETAILS
...
LICENCE
...
DOCUMENTATION
...
...
===================================================================================================
""")
#------------------------------------------------------------------------------------------------
def error(string):
'''
prints the string in red.
'''
#import sys
#sys.stdout.write("\033[39;31m" + string + "\033[39;29m") # Red color string
sys.stdout.write( string ) # No color string
def version(version_num):
'''
prints version number.
'''
#import sys
#versionString = sys.argv[0].split('/')[-1] + ' ' + version_num
versionString = '{0} ({0}-{1}) [-]'.format(sys.argv[0].split('/')[-1], version_num)
print(versionString)
def main(arguments):
# ----------------------------------------------------------------------------------------
'''
This is the main function doing the full loop.
'''
import getopt
### sys.path.insert(0,"/lhome/jperea/ppsdev/pcms/pipeline/python/")
import ModuleUtil
version_num = __version__ # Version number to be updated
verbose = False
allEnergyRangeCorrLC = ''
SrcLCsperEnergyBand = ''
BkgLCsperEnergyBand = ''
outfits = 'output.fits'
pn_timing = False
plotting = True
points_density = 512
plotfile = "figure.png"
plotformat = 'png'
interactplot = False
arguments = arguments.split()
progName = arguments[0]
argsList = arguments[1:]
try:
opts, args = getopt.getopt(argsList, "vo:h", ["allEnergyRangeCorrLC=", "SrcLCsperEnergyBand=", "BkgLCsperEnergyBand=", "outfits=", \
"plotting=", "plotfile=", "plotformat=", "interactplot=", \
"version", "verbose=", "help"])
except getopt.GetoptError as err:
usage()
error("ERROR: Could not understand the options!\n\n")
print(err) # will print something like "option -a not recognized"
print("================================================================================================")
return
for o, a in opts:
if o in ("-v", "--version"):
version(version_num)
return 0
if o in ("-V", "--verbose"):
verbose = ModuleUtil.varCheck(client = progName, inVar = a).checkBool()
if o in ("-h", "--help"):
usage()
return 0
if o in ("--allEnergyRangeCorrLC"):
allEnergyRangeCorrLC = a
if o in ("--SrcLCsperEnergyBand"):
SrcLCsperEnergyBand = a
if o in ("--BkgLCsperEnergyBand"):
BkgLCsperEnergyBand = a
if o in ("-o", "--outfits"):
outfits = a
if o in ("--plotting"):
plotting = ModuleUtil.varCheck(client = progName, inVar = a).checkBool()
if o in ("--plotfile"):
plotfile = a
if o in ("--plotformat"):
plotformat = a
# If plotfile name already has a valid extension (png or pdf): keep it.
if plotfile.lower().endswith(('.png', '.pdf')):
plotfile = plotfile
else:
plotfile = plotfile + '.' + plotformat
# Convert input string in lists
allEnergyRangeCorrLC = allEnergyRangeCorrLC.split(',')
SrcLCsperEnergyBand = SrcLCsperEnergyBand.split(',')
BkgLCsperEnergyBand = BkgLCsperEnergyBand.split(',')
if verbose:
import datetime
time1 = datetime.datetime.utcnow()
print('\nImporting main modules ...\t', time1)
import warnings
import numpy as np
from astropy.table import Table, join, vstack
from astropy.io import fits
from astropy.time import Time
from astropy.coordinates import SkyCoord
from astropy import units as u
if verbose:
time2 = datetime.datetime.utcnow()
print('......Done\t\t\t', time2, '\nImport modules time = ',time2 - time1, '\n')
# ##################################################################################################################################################
uncorrectedSrcRateSet_all_set = SrcLCsperEnergyBand
uncorrectedBkgRateSet_all_set = BkgLCsperEnergyBand
correctedAllEnergyRange_all_set = allEnergyRangeCorrLC
# ##################################################################################################################################################
#
# Inputs for the Python task:
#
# - correctedTemp_all_set[0]
# - uncorrectedSrcRateSet_all_set
# - uncorrectedBkgRateSet_all_set
#
# Previous has to be done from the Perl module: LCMerge.pm :
#
# Input for each SRC and EXPOSURE ?:
#
# WARNING: SRC_NUM = 0 ==> Timing
#
# Number of Exposures
#
#allInputLCs = uncorrectedSrcRateSet_all_set + uncorrectedBkgRateSet_all_set + correctedAllEnergyRange_all_set
exposures = []
correctedLCfiles = {}
uncorrectedSrcRateSets = {}
uncorrectedBkgRateSets = {}
for file in correctedAllEnergyRange_all_set:
exp = fits.getval(file, 'EXPIDSTR', 0)
correctedLCfiles[exp] = file
if exp not in exposures:
exposures.append(exp)
for exp in exposures:
uncorrectedSrcRateSets[exp] = []
uncorrectedBkgRateSets[exp] = []
skipped_exps = []
for file in uncorrectedSrcRateSet_all_set: # Dict per exposures containing all the uncorrected LCs, bamds: 1,2,3,4,5. Source
exp = fits.getval(file, 'EXPIDSTR', 0)
if exp in exposures:
uncorrectedSrcRateSets[exp].append(file)
else:
skipped_exps.append(exp)
print('There is not corrected LC for exposures', set(skipped_exps),'. Skipping this exposure.' )
for file in uncorrectedBkgRateSet_all_set: # Dict per exposures containing all the uncorrected LCs, bamds: 1,2,3,4,5. Background
exp = fits.getval(file, 'EXPIDSTR', 0)
if exp in exposures:
uncorrectedBkgRateSets[exp].append(file)
# Corrected Time Series ( from 'epiclccorr' )
#
print('--------- Corrected Light Curves. All Energy range -----------------------------------------------')
bkgratio = {}
TIMEDEL = {}
exposureTimes = []
correctTIME = {}
correctRATE = {}
correctERROR = {}
FRACEXP = {}
correctBACKV = {}
correctBACKE = {}
for exp, correctedLCfile in correctedLCfiles.items():
print(exp, '\t', correctedLCfile)
with fits.open( correctedLCfile ) as correctedRate:
bkgratio[exp] = correctedRate['RATE'].header['BKGRATIO']
TIMEDEL[exp] = correctedRate['RATE'].header['TIMEDEL']
exposureTimes.append(correctedRate['RATE'].header['EXPOSURE'])
correctTIME[exp] = correctedRate['RATE'].data['TIME']
correctRATE[exp] = correctedRate['RATE'].data['RATE']
correctERROR[exp] = correctedRate['RATE'].data['ERROR']
FRACEXP[exp] = correctedRate['RATE'].data['FRACEXP']
correctBACKV[exp] = correctedRate['RATE'].data['BACKV']
correctBACKE[exp] = correctedRate['RATE'].data['BACKE']
print('\tBKGRATIO from Corrected TS = ', bkgratio[exp], '\n')
totalExposureTime = sum( exposureTimes )
print('\tTotal Exposure times from all Corrected TSs = ', totalExposureTime)
print('--------- Corrected Light Curves. All Energy range -----------------------------------------------')
#-------------------------------------------------------------------
# Return a Dict where :
# keys = Exposures
# values = Data tables ( 'TIME', 'RATE', 'RATE_ERR' )
#
def dict_of_tables_from_files( LCs_files ):
output = {}
energy_band_ranges = {}
band_name = 99
for file in LCs_files:
with fits.open( file ) as RateSet:
try:
pilo = RateSet['RATE'].header['CHANMIN']
pihi = RateSet['RATE'].header['CHANMAX']
except:
print('CHANMIN or CHANMAX keywords missing in file.')
print('Mandatory keywords to know what is the energy range of each input Time Series. Exiting.')
sys.exit(1)
try:
id_band = RateSet['RATE'].header['ID_BAND']
except KeyError:
band_name += 1
id_band = band_name
print( 'ID_BAND keyword missing in file. An ID_BAND created: ', id_band )
#print('\t', file, '\t\tID_BAND, pilo, pihi = ', id_band, pilo, pihi )
energy_band_ranges[id_band] = (pilo, pihi) # ID_BAND = (pilo, pihi) definition
t = RateSet['RATE'].data['TIME']
r = RateSet['RATE'].data['RATE']
Dr = RateSet['RATE'].data['ERROR']
table = Table([t, r, Dr], names=['TIME', 'RATE', 'RATE_ERR'])
output[ id_band ] = table
return (output, energy_band_ranges) # output = { '1': table, '2': table, ..., '5': table }, where table = ['TIME', 'RATE', 'RATE_ERR'] columns
# energy_band_ranges = { 'id_band_1': (pilo, pihi), 'id_band_2': (pilo, pihi), ...}
#-------------------------------------------------------------------
# Find the Min and Max values of RATE* columns
# (for plotting)
#
def getMinMax(hdu):
minmax = {}
allValues = []
for exp, hdul in hdu.items():
table = hdul.data
minmax[exp] = []
for i in ['', '1','2','3','4','5']:
minmax[exp].append( min( table['RATE'+i] ) )
minmax[exp].append( max( table['RATE'+i] ) )
allValues += minmax[exp]
#print('\n', allValues, '\n', len(allValues) )
ylimits = [ min(allValues), max(allValues) ]
return ylimits
#-------------------------------------------------------------------
# Plot all Time Series in different panels
#
def plot_all_LCs( hdu, energy_ranges = '' ):
import matplotlib
matplotlib.use('Agg') # To avoid problems with DISPLAY variable
# raise RuntimeError('Invalid DISPLAY variable')
# RuntimeError: Invalid DISPLAY variable
import matplotlib.pyplot as plt
for exp, bands in energy_ranges.items():
n_bands = len(bands)
# Plot Time Series only if the number of Energy Bands is < 5
if n_bands > 5:
print('Exposure '+ exp + ' has more than 5 energy bands. Plot not produced.')
sys.exit(0)
#fig, ax = plt.subplots(3, 2)
fig, ax = plt.subplots(6, sharex=True )
fig.suptitle('TimeSeries per Energy Band', fontsize=12)
fig.set_size_inches(8.0,10.0)
#pos = [(x, y) for x in [0,1,2] for y in [0,1]] # Panel Positions array
#pos = [0,1,2,3,4]
pos = [1,2,3,4,5]
colors = {0: 'tab:red',
1: 'tab:orange',
2: 'tab:green',
3: 'tab:blue',
4: 'tab:purple'}
###### ylimits = getMinMax( hdu )
new_time_bin = (TSTOP - TSTART) / ( points_density - 1 ) # (Total Exposure time. All Exposures included) / ( # points in plot )
t_ref = TSTART
#=====================================================================================================================
# Reference Time in UTC:
# Reference values from the first Exposure EPIC LC
# timesys, timeref, mjdref = TT LOCAL 50814.0
if 'mjdref' in locals() and 'timesys' in locals():
t_ref_utc = Time(mjdref + t_ref/86400,format='mjd',scale=(timesys).lower(), precision=0).utc.iso # Conversion to UTC
if verbose: print('Minimum TSTART from all Timeseries: Reference Time (UTC) = ', t_ref_utc)
else:
t_ref_utc = t_ref # No UTC conversion
if verbose: print('Minimum TSTART from all Timeseries: Reference Time could not be converted to UTC = ', t_ref)
#=====================================================================================================================
for exp, hdul in hdu.items():
if verbose: print('Plotting data from Exposure : ', exp)
table = hdul.data
time_bin = TIMEDEL[exp]
t_bin_frac = int(np.round( new_time_bin / time_bin ))
if t_bin_frac <= 1:
rebin = False # ==> No re-binning
print(' No re-binning needed for this Light Curve. Input / Output time bins ratio : ', t_bin_frac)
else:
rebin = True
# All Energy range (Band 8) :
tsr = Table([ table['TIME'], table['RATE'], table['ERROR'], table['FRACEXP'] ], names=['TIME', 'RATE', 'ERROR', 'FRACEXP'])
if rebin:
if verbose: print('\tStart re-binning for ALL Bands (band = 8) ...')
ntsr = ModuleUtil.TSrebin(tsr, time_bin, new_time_bin) # import module 'ModuleUtil' required
else:
ntsr = tsr
ntsr['TIME'] = ntsr['TIME'] - t_ref
ranges = energy_ranges[exp]
sortedEnergyRanges = sorted( ranges, key=ranges.get ) # Energy ID_BAND's sorted by Energy
#ax[5].set_title('Band 8: {}-{} eV'.format( ranges[sortedEnergyRanges[0]][0] , ranges[sortedEnergyRanges[-1]][1]), fontsize=10, fontweight='bold')
ax[0].annotate ('Band 8: {}-{} eV'.format( ranges[sortedEnergyRanges[0]][0] , ranges[sortedEnergyRanges[-1]][1]), xy=(0.01, 0.03), xycoords='axes fraction', fontsize=10)
ax[0].errorbar( ntsr['TIME'], ntsr['RATE'], ntsr['ERROR'], fmt='.', c='k', ms=3, ecolor='k')
#ax[2,1].set_ylim( ylimits )
#xticks = ax[5].get_xticks()
idx = -1
for id_band in sortedEnergyRanges:
idx += 1
#if verbose: print('idx: ', idx, '\tpos: ',pos[idx], '\tid_band: ',id_band, '\tranges: ',ranges[id_band])
# Re-Binning for Plot ---------------
tsr = Table([ table['TIME'], table['RATE'+str(id_band)], table['RATE'+str(id_band)+'_ERR'], table['FRACEXP'] ], names=['TIME', 'RATE', 'ERROR', 'FRACEXP'])
if rebin:
if verbose: print('\tStart re-binning for ID_BAND: ', id_band)
ntsr = ModuleUtil.TSrebin(tsr, time_bin, new_time_bin) # import module 'ModuleUtil' required
else:
ntsr = tsr
ntsr['TIME'] = ntsr['TIME'] - t_ref
if verbose: print('\t\tRe-Binning, Exp: ' + exp + '\tID_BAND: '+ str(id_band) + '\ttime_bin = ', time_bin, '\tnew_time_bin = ', new_time_bin)
#------------------------------------
#ax[ pos[idx] ].set_title('Band '+str(id_band)+': {}-{} eV'.format(ranges[id_band][0], ranges[id_band][1]), fontsize=10)
ax[ pos[idx] ].annotate('Band '+str(id_band)+': {}-{} eV'.format(ranges[id_band][0], ranges[id_band][1]), xy=(0.01, 0.03), xycoords='axes fraction', fontsize=10)
ax[ pos[idx] ].errorbar( ntsr['TIME'], ntsr['RATE'], ntsr['ERROR'], fmt='.', c=colors[idx], ms=3, ecolor='k')
#ax[ pos[idx] ].scatter( table['TIME'], table['RATE'+str(id_band)], c=colors[idx], s=3)
#ax[ pos[idx] ].plot ( table['TIME'], table['RATE'+str(id_band)], c=colors[idx], lw=1)
#ax[ pos[idx] ].scatter( table['TIME'], table['BACK'+str(id_band)+'V'], c='k', s=2)
#ax[ pos[idx] ].set_xticks(xticks)
fig.subplots_adjust(hspace=0.0) # Horizontal shift between plots
plt.annotate('OBS_ID: {0}'.format(obs_id), xy=(0.13, 0.93), xycoords='figure fraction')
plt.annotate('Instrument: {0}'.format(instrument), xy=(0.13, 0.91), xycoords='figure fraction')
if 'srcnum' in locals():
plt.annotate('EPIC SRCNUM: {0}'.format(srcnum), xy=(0.58, 0.93), xycoords='figure fraction')
plt.annotate('RA,Dec: {0}'.format(coords_hmsdms), xy=(0.58, 0.91), xycoords='figure fraction')
if pn_timing:
plt.annotate('EPN Timing', xy=(0.58, 0.93), xycoords='figure fraction')
plt.annotate("rate (count/s)", xy=(0.05, 0.6), xycoords='figure fraction', rotation=90)
plt.xlabel('time (s) after {0} (UTC)'.format(t_ref_utc))
#plt.show()
print('\nSaving figure ', plotfile, ' ...')
#plt.savefig(figure, orientation='portrait', papertype='a4')
plt.savefig(plotfile)
#-------------------------------------------------------------------
print('\n--------- Un-Corrected Light Curves per Energy band -----------------------------------------------')
SrcUncorrected_Band = {}
BkgUncorrected_Band = {}
src_energy_ranges = {}
bkg_energy_ranges = {}
print('Reading files for Source ...')
for exp, uncorrectedSrcRates_per_Band in uncorrectedSrcRateSets.items():
(SrcUncorrected_Band[exp], src_energy_ranges[exp]) = dict_of_tables_from_files( uncorrectedSrcRates_per_Band ) # SrcUncorrected_Band['EXPOSURE'] = [ Ban 1 table, Band 2 table, ... ]
print('Reading files for Background ...')
for exp, uncorrectedBkgRates_per_Band in uncorrectedBkgRateSets.items():
(BkgUncorrected_Band[exp], bkg_energy_ranges[exp]) = dict_of_tables_from_files( uncorrectedBkgRates_per_Band ) # BkgUncorrected_Band['EXPOSURE'] = [ Ban 1 table, Band 2 table, ... ]
#
# Check Energy bands of Src and Bkg files
# If Energy bands in Src an Bkg are not the same we can not get the Corrected Time Series per Energy band
#
print('Checking Src and Bkg Energy Bands ...')
for exp, bands in src_energy_ranges.items():
if src_energy_ranges[exp] != bkg_energy_ranges[exp]:
print('Energy bands for Source and Background differ in exposure', exp, '. Exiting.')
print('\t\t{band: (pilo, pihi), band: (pilo, pihi), ...}')
print('Src bands:\t', src_energy_ranges[exp])
print('Bkg bands:\t', bkg_energy_ranges[exp])
sys.exit(1)
for id_band, ranges in bands.items():
print(exp, '\tID_BAND: ',id_band, '\tEnergy range: ',ranges, 'eV')
print('-----------------------------------------------------------')
np.seterr(divide='ignore', invalid='ignore') # Set Treatment for "division by zero" and "invalid floating-point operation" to 'ignore'
correctedTables = {}
for exp in SrcUncorrected_Band.keys():
if verbose: print('\t\t EXPOSURE = ', exp)
if verbose: print('\t\t BKGRATIO = ', bkgratio[exp])
if verbose: print('\t\t FRACEXP = ', FRACEXP[exp])
correctedTables[exp] = {}
bands = src_energy_ranges[exp].keys()
for band in bands:
if verbose: print('\t\t\t Band = ', band)
#---------------------------------------------------------------------
# Src and Bkg TIME, RATE, and ERROR's
#
t = SrcUncorrected_Band[exp][band]['TIME']
r = SrcUncorrected_Band[exp][band]['RATE']
Dr = SrcUncorrected_Band[exp][band]['RATE_ERR']
b = BkgUncorrected_Band[exp][band]['RATE']
Db = BkgUncorrected_Band[exp][band]['RATE_ERR']
srcRATE = (r - (b/bkgratio[exp])) / (FRACEXP[exp])
srcRATE_ERROR = ( 1/FRACEXP[exp] ) * ( (Dr)**2 + (Db/bkgratio[exp])**2 )**(1/2)
bkgRATE = b/( bkgratio[exp]*FRACEXP[exp] )
bkgRATE_ERROR = Db/( bkgratio[exp]*FRACEXP[exp] )
# Corrrected Time Series per Energy Band:
bandstr = str(band)
correctedTables[exp][band] = Table([t, srcRATE, srcRATE_ERROR, bkgRATE, bkgRATE_ERROR], names=['TIME', 'RATE'+bandstr, 'RATE'+bandstr+'_ERR', 'BACK'+bandstr+'V', 'BACK'+bandstr+'E'])
if verbose: print('\t\t\t\t Table size = ', len( correctedTables[exp][band] ))
#
#---------------------------------------------------------------------
# Write down 'finalTableJoin' into the original Time Series file : correctedLCfile
#-----------------------------------------------------------------------------------------------------------------
#
# 1) Copy the file:
#
from shutil import copy
intermediate = 'intermediate/' # To be modified in a final SAS task
# Loop running over the EXPOSUREs
# correctedLCfiles (product/)
# 1) Copy to intermediate/
# 2) Add RATE1, RATE2, etc to Corrected from each EXPOSURE
# 3) Save it as temporal file. One file per EXPOSURE
#
hdu = {}
orig_cols = {}
new_cols = {}
orig_header = {}
print('Creating intermediate Light Curve files per exposure')
for exp, correctedLCfile in correctedLCfiles.items(): # Loop running over the EXPOSUREs and the proper files
# Copy
filename = correctedLCfile.split('/')[-1]
file = intermediate + 'srcTimeSeriesEnergyBands-correctedLightCurve-' + filename.split('.')[0] + '-plus_Rates_per_Energy_band.ftz'
if verbose: print(exp, '\tCopy ', correctedLCfile, 'to', file)
copy( correctedLCfile, file )
print('Populating intermediate file ', file , 'with Corrected Lihgt Curves per Energy band')
# i) Open/Read the All Energy Range Time Series File (Band 8): *SRCTSR* -> hdu table
# ii) Join the Corrected Time Series per Energy band (correctedTables[band]) to the previous hdu table
# iii) Paste the new column values to the input FITS file
with fits.open( file, mode='update' ) as hdul:
orig_header[exp] = hdul['RATE'].header
hdu[exp] = hdul['RATE'].data
orig_cols[exp] = hdu[exp].columns
# Columns Definitions:
fitsColumnDefinitions = []
for band in correctedTables[exp].keys():
hdu[exp] = join(hdu[exp] , correctedTables[exp][band], keys='TIME', join_type='right')
bandstr = str(band)
fitsColumnDefinitions.append( fits.Column( name='RATE'+bandstr , format='E', unit='count/s', array=hdu[exp]['RATE'+bandstr]) )
fitsColumnDefinitions.append( fits.Column( name='RATE'+bandstr+'_ERR', format='E', unit='count/s', array=hdu[exp]['RATE'+bandstr+'_ERR']) )
fitsColumnDefinitions.append( fits.Column( name='BACK'+bandstr+'V' , format='E', unit='count/s', array=hdu[exp]['BACK'+bandstr+'V']) )
fitsColumnDefinitions.append( fits.Column( name='BACK'+bandstr+'E' , format='E', unit='count/s', array=hdu[exp]['BACK'+bandstr+'E']) )
fitsColumnDefinitions.append( fits.Column( name='EXP_ID' , format='A4' , array= [exp] * len(hdu[exp]) ) )
fitsColumnDefinitions.append( fits.Column( name='TIMEDEL' , format='D', unit='s' , array= [ TIMEDEL[exp] ] * len(hdu[exp]) ) )
new_cols[exp] = fits.ColDefs( fitsColumnDefinitions )
# Update the RATE extension => hdu = hdul['RATE'].data
hdu[exp] = fits.BinTableHDU.from_columns( orig_cols[exp] + new_cols[exp], header=orig_header[exp] )
hdul['RATE'] = hdu[exp]
hdul.flush() # Flush to the Complete HDUList == hdul
print()
#-----------------------------------------------------------------------------------------------------------------
################################################################################################################
#
# Tables from table = Table.read( file, hdu='RATE' )
#
print('Merging all intermediate Corrected Light Curves into the final product')
table = {}
allTables = []
times = {}
for exp, correctedLCfile in correctedLCfiles.items():
#file = intermediate + correctedLCfile.split('/')[-1]
filename = correctedLCfile.split('/')[-1]
file = intermediate + 'srcTimeSeriesEnergyBands-correctedLightCurve-' + filename.split('.')[0] + '-plus_Rates_per_Energy_band.ftz'
if verbose: print( 'Reading ...', exp, ' ... ', file )
with warnings.catch_warnings():
warnings.simplefilter("ignore", category=Warning) # To ignore: WARNING: UnitsWarning: 'cts' did not parse as fits unit
table[exp] = Table.read( file, hdu='RATE' )
times[exp] = ( table[exp]['TIME'][0], table[exp]['TIME'][-1] ) # First and Last time values
# Vertical stack of all Exposure tables
# 1) Sort the exposure names ( dict keys ) by Start Time ( times dict values )
sortedExposures = sorted( times, key=times.get )
# 2) Vertical stacking of teh tables ( table[exp] )
for exp in sortedExposures:
allTables.append( table[exp] )
TSTART = allTables[0]['TIME'][0] # First Time value from the 1st LC
TSTOP = allTables[-1]['TIME'][-1] # Last Time value from the last LC
allTable = vstack( allTables, metadata_conflicts='silent' )
# Convert 'Table object' to 'BinTableHDU object'
allTableHDU = fits.table_to_hdu( allTable ) # BinTableHDU object which is complient with writing the FITS Extension and Flush :
# Lost: Header keywords comments
# Write the Final Time Series file:
#
# Get the 1st Exposure file:
exp_1st = sortedExposures[0] # 1st exposure == 1st exposure in time
#file_1st = intermediate + correctedLCfiles[exp_1st].split('/')[-1]
file_1st = correctedLCfiles[exp_1st]
print('Copy (1st exposure Light Curve) ', file_1st, ' to ', outfits)
# Copy
copy( file_1st, outfits )
print('\nWriting all into the Final FITS file : ', outfits, '\n')
# Write down the final stacked table into the 1st Exposure FITS file:
with fits.open( outfits, mode='update' ) as hdul:
# --- Get some important info from the original LC file ---
if file_1st.find('PN') != -1 and file_1st.find('SRCTSR0000.F') != -1:
pn_timing = True
instrument = hdul['RATE'].header['INSTRUME']
obs_id = hdul['RATE'].header['OBS_ID']
objectID = hdul['RATE'].header['OBJECT']
#print('===========> ', instrument, obs_id, objectID, pn_timing)
# Get keywords fo Source coordinates
if instrument == 'EPN' or instrument == 'EMOS1' or instrument == 'EMOS2':
try:
ra_deg = hdul['RATE'].header['SRC_RA']
dec_deg = hdul['RATE'].header['SRC_DEC']
srcnum = hdul['RATE'].header['SRCNUM']
coord = SkyCoord(ra=ra_deg*u.degree, dec=dec_deg*u.degree)
coords_hmsdms = coord.to_string('hmsdms')
#print('===========> ', srcnum, ra_deg, dec_deg)
except KeyError:
pass
# Get keywords for Time conversion
try:
timesys = hdul['RATE'].header['TIMESYS']
timeref = hdul['RATE'].header['TIMEREF']
mjdref = hdul['RATE'].header['MJDREF']
#print('===========> ', timesys, timeref, mjdref)
except KeyError:
pass
# Write the Final Table into the FITS HDU extension 'RATE' = hdul['RATE']
hdul['RATE'] = allTableHDU
hdul.flush()
# - - - -
# - - - - FITS Header modifications
with fits.open( outfits, mode='update' ) as hdul:
header = hdul['RATE'].header
header['EXP_ID'] = ' '.join(exposures)
header['EXPIDSTR'] = ' '.join(exposures)
header['TSTART'] = TSTART
header['TSTOP'] = TSTOP
header['EXPOSURE'] = totalExposureTime
header.remove('SLCTEXPR')
header.remove('TIMEDEL')
header.remove('BKGRATIO')
try:
header.remove('FLCUTTHR')
except:
pass
header.remove('AVRATE')
header.remove('CHISQUAR')
header.remove('CHI2PROB')
header.remove('N_POINTS')
header.remove('FVAR')
header.remove('FVARERR')
i=0
for exp in correctedLCfiles.keys():
i +=1
header['BKGRAT'+str(i)] = ( bkgratio[exp], 'Ratio of BKG to SRC collection areas for ' + exp)
for id_band, ranges in src_energy_ranges[exp].items():
chanmin_label = 'CHANMIN' + str(id_band)
chanmax_label = 'CHANMAX' + str(id_band)
header['CHAMI'+str(id_band)] = ( ranges[0], 'Minimum of the energy column ID_BAND '+ str(id_band) +', ' + exp)
header['CHAMA'+str(id_band)] = ( ranges[1], 'Maximum of the energy column ID_BAND '+ str(id_band) +', ' + exp)
# Write the Final modified Header into the FITS HDU extension 'RATE' = hdul['RATE']
hdul['RATE'].header = header
hdul.flush()
# - - - -
if plotting:
print ('Plotting Time Series ...')
#dataFromFile = Table.read( outfits, hdu='RATE' )
#plot_all_LCs( dataFromFile )
plot_all_LCs( hdu, energy_ranges = src_energy_ranges )
sys.exit()
##################################################################################################################
##################################################################################################################
# ------------------------------------------------------------------------------------------------------------------------
if __name__=='__main__':
import sys
main(" ".join(sys.argv))
# ------------------------------------------------------------------------------------------------------------------------