Commit 24c5a8f3 authored by Carlo Ferrigno's avatar Carlo Ferrigno
Browse files

v 1.3.6 : ploting of lightcurves improved. Patched a bug of ddosa (possibly)

that uses names of sources in light curves with spaces at the place of + and -
parent 1b67f395
......@@ -55,7 +55,7 @@ class INTEGRALwrapper(object):
except:
raise ConnectionError
def long_scw_list_call(self, in_total_scw_list, s_max=50, wait=True, **arguments):
def long_scw_list_call(self, in_total_scw_list, s_max=50, wait=True, sleep_time=10, **arguments):
import time
total_scw_list = sorted(in_total_scw_list)
......@@ -93,7 +93,7 @@ class INTEGRALwrapper(object):
disp_by_call[ys] = oda_api.api.DispatcherAPI(url=self.disp.url, wait=False)
_disp = disp_by_call[ys]
#print(type(_disp))
data = data_by_call.get(ys, None)
if data is None or not _disp.is_failed():
......@@ -120,7 +120,7 @@ class INTEGRALwrapper(object):
n_poll += 1
if not wait:
return None
time.sleep(5)
time.sleep(sleep_time)
loc_keys = data_by_call.keys()
......@@ -330,18 +330,21 @@ class INTEGRALwrapper(object):
@staticmethod
def plot_lc(combined_lc, source_name, systematic_fraction=0, ng_sig_limit=3):
# In LC name has no "-" nor "+" ??????
patched_source_name = source_name.replace('-', ' ').replace('+', ' ')
from scipy import stats
hdu = None
for j, dd in enumerate(combined_lc._p_list):
if dd.meta_data['src_name'] == source_name:
IND_src_combined = j
INTEGRALwrapper.__log.debug(dd.meta_data['src_name'])
if dd.meta_data['src_name'] == source_name or dd.meta_data['src_name'] == patched_source_name:
for ii, du in enumerate(dd.data_unit):
if 'LC' in du.name:
IND_lc_combined = ii
hdu = du.to_fits_hdu()
for dd in combined_lc._p_list:
if dd.meta_data['src_name'] == source_name:
INTEGRALwrapper.__log.info('Source ' + source_name)
hdu = combined_lc._p_list[IND_src_combined].data_unit[IND_lc_combined].to_fits_hdu()
if hdu is None:
INTEGRALwrapper.__log.info('Source ' + source_name + ' not found in lc, not plotting')
return None
x = hdu.data['TIME']
y = hdu.data['RATE']
......@@ -357,8 +360,10 @@ class INTEGRALwrapper(object):
meany = np.sum(y / dy ** 2) / np.sum(1. / dy ** 2)
err_mean = np.sum(1 / dy ** 2)
std_dev = np.std(y)
fig = plt.figure()
_ = plt.errorbar(x, y, yerr=dy, marker='o', capsize=0)
_ = plt.errorbar(x, y, yerr=dy, marker='o', capsize=0, linestyle='', label='Lightcurve')
_ = plt.axhline(meany, color='green', linewidth=3)
_ = plt.xlabel('Time [IJD]')
_ = plt.ylabel('Rate')
......@@ -367,9 +372,35 @@ class INTEGRALwrapper(object):
prob_limit = stats.norm().sf(ng_sig_limit)
chi2_limit = stats.chi2(ndof).isf(prob_limit)
band_width = np.sqrt(chi2_limit / err_mean)
INTEGRALwrapper.__log.debug('%g %g %g %g %g %g %g' %(
meany, err_mean, std_dev, prob_limit, chi2_limit, band_width, ng_sig_limit))
_ = plt.axhspan(meany - band_width, meany + band_width, color='green', alpha=0.3,
label=f'{ng_sig_limit} $\sigma$, {100 * systematic_fraction}% syst')
label=f'{ng_sig_limit} $\sigma_m$, {100 * systematic_fraction}% syst')
_ = plt.axhspan(meany - std_dev*ng_sig_limit, meany + std_dev*ng_sig_limit,
color='cyan', alpha=0.3,
label=f'{ng_sig_limit} $\sigma_d$, {100 * systematic_fraction}% syst')
_ = plt.legend()
ind = (y - band_width)/dy > ng_sig_limit
if np.sum(ind) > 0:
_ = plt.plot(x[ind], y[ind], marker='x', color='red', linestyle='', markersize=10)
INTEGRALwrapper.__log.info('We found positive excesses on the lightcurve at times')
good_ind = np.where(ind)
#print(good_ind[0][0:-1], good_ind[0][1:])
old_time = -1
for i,j in zip(good_ind[0][0:-1], good_ind[0][1:]):
#print(i,j)
if j-i > 2:
if x[i] != old_time :
INTEGRALwrapper.__log.info('%f' % x[i])
INTEGRALwrapper.__log.info('%f' % (x[j]))
# else:
# INTEGRALwrapper.__log.info('%f' % ((x[i]+x[j])/2))
old_time = x[j]
plot_title = source_name
_ = plt.title(plot_title)
......@@ -611,7 +642,10 @@ class INTEGRALwrapper(object):
@staticmethod
def write_lc_fits_files(lc, source_name, subcases_pattern, output_dir='.'):
lcprod = [l for l in lc._p_list if l.meta_data['src_name'] == source_name]
# In LC name has no "-" nor "+" ??????
patched_source_name = source_name.replace('-', ' ').replace('+', ' ')
lcprod = [l for l in lc._p_list if l.meta_data['src_name'] == source_name or \
l.meta_data['src_name'] == patched_source_name]
if (len(lcprod) < 1):
INTEGRALwrapper.__log.warning("source %s not found in light curve products" % source_name)
return "none", 0, 0, 0
......@@ -621,16 +655,20 @@ class INTEGRALwrapper(object):
"source %s is found more than once light curve products, writing only the first one" % source_name)
instrument = lcprod[0].data_unit[1].header['INSTRUME']
if instrument == 'IBIS':
ind_extension = 1
else:
ind_extension = 2
lc_fn = output_dir + "/%s_lc_%s_%s.fits" % (instrument, source_name.replace(' ', '_'), subcases_pattern)
lcprod[0].write_fits_file(lc_fn)
ff = fits.open(lc_fn) # , mode='update'
mjdref = float(ff[1].header['MJDREF'])
tstart = float(ff[1].header['TSTART']) + mjdref
tstop = float(ff[1].header['TSTOP']) + mjdref
mjdref = float(ff[ind_extension].header['MJDREF'])
tstart = float(ff[ind_extension].header['TSTART']) + mjdref
tstop = float(ff[ind_extension].header['TSTOP']) + mjdref
try:
exposure = ff[1].header['EXPOSURE']
exposure = float(ff[ind_extension].header['EXPOSURE'])
except:
exposure = -1
......@@ -721,36 +759,6 @@ class INTEGRALwrapper(object):
return spec_fn, tstart, tstop, exposure
@staticmethod
def write_lc_fits_files(lc, source_name, subcases_pattern, output_dir='.'):
lcprod = [l for l in lc._p_list if l.meta_data['src_name'] == source_name]
if (len(lcprod) < 1):
INTEGRALwrapper.__log.warning("source %s not found in light curve products" % source_name)
return "none", 0, 0, 0
if (len(lcprod) > 1):
INTEGRALwrapper.__log.warning(
"source %s is found more than once light curve products, writing only the first one" % source_name)
instrument = lcprod[0].data_unit[1].header['INSTRUME']
lc_fn = output_dir + "/%s_lc_%s_%s.fits" % (instrument, source_name.replace(' ', '_'), subcases_pattern)
lcprod[0].write_fits_file(lc_fn)
ff = fits.open(lc_fn) # , mode='update'
mjdref = ff[1].header['MJDREF']
tstart = float(ff[1].header['TSTART']) + mjdref
tstop = float(ff[1].header['TSTOP']) + mjdref
try:
exposure = ff[1].header['EXPOSURE']
except:
exposure = -1
ff.close()
return lc_fn, tstart, tstop, exposure
@staticmethod
def show_spectral_products(summed_data):
for dd, nn in zip(summed_data._p_list, summed_data._n_list):
......
......@@ -19,7 +19,7 @@ include_package_data=True
scripts_list = glob.glob('./bin/*')
setup(name='oda_integral_wrapper',
version="1.3.5",
version="1.3.6",
description='wrapper for INTEGRAL analysis using the API plugin for CDCI online data analysis',
author='Carlo Ferrigno',
author_email='carlo.ferrigno@unige.ch',
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment