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