#!/usr/bin/env python
version_num = "1.1" # Version number to be updated
#
# VERSION = '1.1';
# author = 'Jose Vicente Perea';
# date = '2017-04-04';
#
#
# ChangeLog
# =========
#
# Version 1.1 - 2017-04-04 (JVP)
# ------------
#
# + Include matplotlib.use('Agg') to avoid issue:
# 'This application failed to start because it could not find or load the Qt platform plugin "xcb"'
#
# Version 1.0 - 2017-03-28 (JVP)
# ------------
#
# + Initial version
# Python script to get the Flaring Background filtering threshold cut
# from the Source Time Series. That cut value will be used to produce
# the Source Spectrum.
# + In addition a diagnostic plot is created. Default: No
#
# ------------------------------------------------------------------------------
#
# usage: srcBkgGti.py [-h] [--version] --srctslist SourceTimeSeries.ds [--outfile OUTFILE] --outformat {png,pdf}
#
# Get the Flaring Background filtering threshold cut from the Source Time Series.
#
# optional arguments:
# -h, --help show this help message and exit
# --version, -V show program's version number and exit
#
# --srctslist INFILE Input Source Time Series FITS file. Default: SourceTimeSeries.ds
# --outfile OUTFILE Output figure file. Default: figure
# --outformat {png,pdf} Output figure format. Default: png
#
# ------------------------------------------------------------------------------
#
import sys,os
import warnings
import argparse
###################
# Inut Options
#
parser = argparse.ArgumentParser(description='Get the Flaring Background filtering threshold cut from the Source Time Series', prog='srcBkgGti.py')
#parser.add_argument('--version', '-v', action='version', version='%(prog)s (%(prog)s-1.0) [-]')
parser.add_argument('--version', '-v', action='version', version='%(prog)s (%(prog)s-{}) [-]'.format(version_num))
parser.add_argument('--srctslist', help='Input Source Time Series FITS file. Default: SourceTimeSeries.ds', default='SourceTimeSeries.ds')
parser.add_argument('--diagplot' , help='Create a diagnostic plot. Default: No', default='No')
parser.add_argument('--outfile' , help='Output figure file. Default: figure. If the output filename already has a valid extension, it will be kept', default='figure')
parser.add_argument('--outformat', help='Output figure format. Default: png', choices=['png', 'pdf'], default='png')
args = parser.parse_args()
#
###################
from astropy.io import fits
from astropy.table import Table
from astropy.coordinates import SkyCoord
from astropy import units as u
import numpy as np
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
srcTSList = args.srctslist
ext = args.outformat
diagPlot = args.diagplot
diagPlotTest = 0
yes = ('yes', 'y')
no = ('no', 'n')
if diagPlot.lower() in yes:
diagPlotTest = 1
elif diagPlot.lower() in no:
print('No plot produced.')
else:
print('Wrong option for --diagplot. No plot produced.')
def isNotNaN(num):
return num == num
# If input file name already has a valid extension (png or pdf): keep it
if args.outfile.lower().endswith(('.png', '.pdf')):
#print('\n ----> outfile name already has a valid format')
figure = args.outfile
else:
figure = args.outfile + '.' + ext
#print(' ----> figure = ', figure, '\n')
if os.path.isfile(srcTSList) and os.access(srcTSList, os.W_OK):
#print("\nInput file " + srcTSList + " exists and has write persmissions.")
pass
else:
print("\nEither file " + srcTSList + " is missing or has not write persmissions.\n")
sys.exit(1)
###################
#
# Read input FITS file
#
with warnings.catch_warnings():
warnings.simplefilter("ignore", category=Warning) # Ignore warnings about non-standard FITS units
tsr = Table.read(srcTSList, hdu=1) # Read FITS file table and converted to Table
tsr_not_sorted = Table.read(srcTSList, hdu=1)
# Writing Keywords:
#
# tsr.meta['MY_KEYWD'] = 'my value'
# tsr.meta['COMMENT'] = ['First comment', 'Second comment', 'etc']
# tsr.write('my_table.fits', overwrite=True)
#
# Reading FITS keywords:
obs_id = tsr.meta['OBS_ID']
instrument = tsr.meta['INSTRUME']
if instrument == 'EPN' or instrument == 'EMOS1' or instrument == 'EMOS2':
try:
ra_deg = tsr.meta['SRC_RA']
dec_deg = tsr.meta['SRC_DEC']
srcnum = tsr.meta['SRCNUM']
coord = SkyCoord(ra=ra_deg*u.degree, dec=dec_deg*u.degree)
coords_hmsdms = coord.to_string('hmsdms')
except KeyError:
pass
# Reading Table data and sort by 'BACKV'
tsr.sort('BACKV') # Sort Table by 'BACKV'
rate = np.nan_to_num(tsr['RATE']) * np.nan_to_num(tsr['FRACEXP'] ) # rate = RATE * FRACEXP == rate weighted by the ExposureFraction
crate = (np.fmax(rate, rate*0)).cumsum() # cummulative rate
erate = np.sqrt( ((np.nan_to_num(tsr['ERROR'])**2)*tsr['FRACEXP']).cumsum() ) # cummulative error
crate = np.ma.masked_array(crate, mask=(erate == 0)) # mask with 0 any NaN
erate = np.ma.masked_array(erate, mask=(erate == 0)) # mask with 0 any NaN
s2nc = np.nan_to_num(crate/erate) # cummulative S/N
imax = s2nc.argmax() # indices of the maximum values
s2nmax = s2nc[imax] #
backvmax = tsr['BACKV'][imax] # Value of 'BACKV' for the maximum S/N
covered_time_fraction = 100.*(imax+1.)/len(s2nc)
s2nlast = s2nc[-1]
diff_max_last = 100*(abs(s2nmax - s2nlast)/s2nmax)
#print('====> ', s2nmax, s2nlast, diff_max_last)
# If diff max-last < 20 % OR max(s2nc) < 10 ==> Include all the the range == no-threshold
if diff_max_last < 20 or s2nmax < 10:
backvmax = tsr['BACKV'][ isNotNaN(tsr['BACKV']) ][-1]
#backvmax = tsr['BACKV'][-1] # Value of 'BACKV' for the last S/N point
covered_time_fraction = 100. # Coverage = 100 %
# Remove entries where 'BACKV' = NaN
# Print result to stdout
#
#print ("Max. rate for maximum S/N =", backvmax," ; fraction of time covered =",100.*(imax+1.)/len(s2nc))
print('Bkg rate threshold cut = {0:.4f}; fraction of time covered = {1:.4f} %'.format(backvmax, covered_time_fraction))
# Write Keyword (FLCUTTHR = backvmax) in the original FITS
#
fits.setval(srcTSList, 'FLCUTTHR', value = backvmax, ext=1)
#-------------------------------------------------------------
# Plot
#
if diagPlotTest:
f = plt.figure(1) # Default figure size = 640x480 = 6.4 x 4.8 inch
#print(f.get_size_inches()) # Return the figure size in inches
#f.set_size_inches(8.0,6.0) # ratio 8/6 = 640/480
f.set_size_inches(8.0,10.0)
#print(f.get_size_inches()) #
#print(f)
plt.subplot(3, 1, 1)
#plt.title('Flaring background filtering\nfor maximum S/N diagnostic plot')
plt.plot(tsr['BACKV'], s2nc, ms=4, ls='-')
plt.plot( backvmax, s2nmax, 'rs', fillstyle='none', ms=8)
plt.axvline( x=backvmax, color='r', ls='--')
plt.xlabel('{0}'.format('Background rate (BACKV)'))
plt.ylabel('{0}'.format('Cummulative S/N'))
plt.annotate('Max. Cumm. S/N: {0:.2f}'.format(s2nmax), xy=(0.4, 0.2), xycoords='axes fraction')
plt.annotate('Diff. Max. Cumm. S/N vs last value: {0:.2f} %'.format(diff_max_last), xy=(0.4, 0.1), xycoords='axes fraction')
time = tsr_not_sorted['TIME'] - tsr_not_sorted['TIME'][0] # Ref. Time 0 at TimeSeries start time
plt.subplot(3, 1, 2)
plt.plot( time ,tsr_not_sorted['RATE'], ms=0.5, ls='-')
plt.ylabel('{0}'.format('Source rate'))
plt.subplot(3, 1, 3)
plt.plot( time ,tsr_not_sorted['BACKV'], ms=0.5, ls='-')
plt.axhline( y=backvmax, color='r', ls='--')
plt.xlabel('{0}'.format('TIME (sec)'))
plt.ylabel('{0}'.format('Background rate (BACKV)'))
f.subplots_adjust(hspace=0.4) # hspace between plots
plt.annotate('Flaring background filtering\nfor maximum S/N diagnostic plot', xy=(0.5, 0.95), xycoords='figure fraction', horizontalalignment='center', size=12)
plt.annotate('OBS_ID: {0}'.format(obs_id), xy=(0.13, 0.91), xycoords='figure fraction')
plt.annotate('EPIC SRCNUM: {0}'.format(srcnum), xy=(0.50, 0.91), xycoords='figure fraction')
plt.annotate('RA,Dec: {0}'.format(coords_hmsdms), xy=(0.50, 0.89), xycoords='figure fraction')
#ext='png'
#figure='figure.' + ext
plt.savefig( figure, format=ext)
#plt.show()
#-------------------------------------------------------------