Skip to content

GitLab

Commit c39cae5d authored by Carlo Ferrigno's avatar Carlo Ferrigno
Browse files

1.3.14 

Tested asynchronous call.
Corrected spurious overwriting when LC is normalized over OSA versions.
parent a3368e94
Hide whitespace changes
Inline Side-by-side
Showing with 120 additions and 80 deletions
oda_integral_wrapper/get_osa10_11_factor.py
View file @ c39cae5d
......@@ -3,9 +3,9 @@ import subprocess
import json
import yaml
@logged
def get_existing_dict():
try:
with open('secret.yaml', 'r') as ff:
secret = yaml.full_load(ff)
......@@ -18,9 +18,11 @@ def get_existing_dict():
else:
cmd = 'git clone https://oauth2:%s@gitlab.astro.unige.ch/oda/integral-osa-additional-parameters.git' % (
secret['token'])
cmd += ';cd integral-osa-additional-parameters;git remote set-url origin https://oauth2:%s@gitlab.astro.unige.ch/oda/integral-osa-additional-parameters.git' % (
secret['token'])
cmd += ';cd integral-osa-additional-parameters;' + \
'git remote set-url origin ' + \
'https://oauth2:%s@gitlab.astro.unige.ch/oda/integral-osa-additional-parameters.git' % (
secret['token'])
get_existing_dict._log.info("It will be possible to upload the recomputed factor")
if (subprocess.call(cmd, shell=True)):
......@@ -38,9 +40,10 @@ def get_existing_dict():
return conversion_dict
@logged
def get_osa10_11_conversion_factor(E1_isgri_keV, E2_isgri_keV, radius=10., oda_platform="staging-1-2",
s_max=50, lc_time_bin=3000):
s_max=50, lc_time_bin=3000):
import os.path
conversion_dict = get_existing_dict()
......@@ -54,7 +57,8 @@ def get_osa10_11_conversion_factor(E1_isgri_keV, E2_isgri_keV, radius=10., oda_p
if compute_factor:
if not os.path.isfile('secret.yaml'):
get_osa10_11_conversion_factor._log.warning("As the file secret.yaml does not exist, I do not perform the computation of the value")
get_osa10_11_conversion_factor._log.warning(
"As the file secret.yaml does not exist, I do not perform the computation of the value")
return -1.0
import oda_integral_wrapper.wrapper
......@@ -72,46 +76,45 @@ def get_osa10_11_conversion_factor(E1_isgri_keV, E2_isgri_keV, radius=10., oda_p
simbad = Simbad.query_object(source_name)
coord = SkyCoord(simbad['RA'], simbad['DEC'], unit=[u.hour, u.deg])
coord.fk5
_ = coord.fk5
get_osa10_11_conversion_factor._log.debug("Coordinates for %s are RA=%.4f, Dec=%.4f" % (source_name, coord.ra.deg[0], coord.dec.deg[0]))
get_osa10_11_conversion_factor._log.debug(
"Coordinates for %s are RA=%.4f, Dec=%.4f" % (source_name, coord.ra.deg[0], coord.dec.deg[0]))
ra = coord.ra.deg[0]
dec = coord.dec.deg[0]
revs=[]
revs.append({'coord': coord, 'tstart': tstart, 'tstop': tstop, 'name': source_name, 'label': 'OSA10.2'})
revs.append({'coord': coord, 'tstart': tstart, 'tstop': tstop, 'name': source_name, 'label': 'OSA11.0'})
revs = [{'coord': coord, 'tstart': tstart, 'tstop': tstop, 'name': source_name, 'label': 'OSA10.2'},
{'coord': coord, 'tstart': tstart, 'tstop': tstop, 'name': source_name, 'label': 'OSA11.0'}]
import astroquery.heasarc
Heasarc = astroquery.heasarc.Heasarc()
def _get_scw_list(ra_obj, dec_obj, rr, start_date, end_date):
R = Heasarc.query_region(
position=SkyCoord(ra_obj, dec_obj, unit='deg'),
radius=f"{rr} deg",
mission='intscw',
time=start_date + " .. " + end_date,
good_isgri=">1000",
)
rr = Heasarc.query_region(
position=SkyCoord(ra_obj, dec_obj, unit='deg'),
radius=f"{rr} deg",
mission='intscw',
time=start_date + " .. " + end_date,
good_isgri=">1000",
)
R.sort('SCW_ID')
rr.sort('SCW_ID')
return R
return rr
def get_scw_list(ra_obj, dec_obj, rr, start_date, end_date):
for i in range(10):
for ii in range(10):
try:
return _get_scw_list(ra_obj, dec_obj, rr, start_date, end_date)
except Exception as e:
get_osa10_11_conversion_factor._log.warning(e)
raise RuntimeError
for i, source in enumerate(revs):
for i, source in enumerate(revs):
get_osa10_11_conversion_factor._log.debug(source['coord'].ra.deg, source['tstart'])
r = get_scw_list(ra, dec, radius, tstart, tstop)
scwids = r['SCW_ID']
scwver = r['SCW_VER']
get_osa10_11_conversion_factor._log.debug(source['name'], ' nscw=%d' % (len(scwids)))
......@@ -119,8 +122,7 @@ def get_osa10_11_conversion_factor(E1_isgri_keV, E2_isgri_keV, radius=10., oda_p
revs[i].update(scwver=scwver)
revs[i].update({'RA_SCX': r['RA'], 'DEC_SCX': r['DEC']})
api_cat={
api_cat = {
"cat_frame": "fk5",
"cat_coord_units": "deg",
"cat_column_list": [
......@@ -163,9 +165,10 @@ def get_osa10_11_conversion_factor(E1_isgri_keV, E2_isgri_keV, radius=10., oda_p
all_lc = {}
for source in revs:
name = source['name']+source['label']
scw_list = [ss+'.'+vv.strip() for ss, vv in zip(np.array(source['scwids'], dtype=str),
np.array(source['scwver'], dtype=str)) if ss.endswith('0')]
name = source['name'] + source['label']
scw_list = [ss + '.' + vv.strip() for ss, vv in zip(np.array(source['scwids'], dtype=str),
np.array(source['scwver'], dtype=str)) if
ss.endswith('0')]
out_str = 'We process %d science windows: ' % (len(scw_list)) + ' '.join(scw_list)
get_osa10_11_conversion_factor._log.info(out_str)
......@@ -173,35 +176,35 @@ def get_osa10_11_conversion_factor(E1_isgri_keV, E2_isgri_keV, radius=10., oda_p
if all_lc.get(name, None) is None:
all_lc.update({name: {'scw_list': scw_list}})
combined_data = wrap.long_scw_list_call(scw_list, s_max=s_max,sleep_time=120,
instrument="isgri",
product='isgri_lc',
E1_keV=E1_isgri_keV,
E2_keV=E2_isgri_keV,
query_type='Real',
osa_version=source['label'],
product_type='Real',
time_bin=lc_time_bin,
selected_catalog=json.dumps(api_cat))
combined_data = wrap.long_scw_list_call(scw_list, s_max=s_max, sleep_time=120,
instrument="isgri",
product='isgri_lc',
E1_keV=E1_isgri_keV,
E2_keV=E2_isgri_keV,
query_type='Real',
osa_version=source['label'],
product_type='Real',
time_bin=lc_time_bin,
selected_catalog=json.dumps(api_cat))
all_lc[name].update({'isgri': combined_data})
wrap.write_lc_fits_files(combined_data, source_name, '%d-%d-%s'%(E1_isgri_keV,E2_isgri_keV, source['label']))
wrap.write_lc_fits_files(combined_data, source_name,
'%d-%d-%s' % (E1_isgri_keV, E2_isgri_keV, source['label']))
t1_10, dt_10, r_10, dr_10 = wrap.get_lc(all_lc['CrabOSA10.2']['isgri'],'Crab')
t1_10, dt_10, r_10, dr_10 = wrap.get_lc(all_lc['CrabOSA10.2']['isgri'], 'Crab')
t1_11, dt_11, r_11, dr_11 = wrap.get_lc(all_lc['CrabOSA11.0']['isgri'], 'Crab')
factor = np.mean(r_10)/np.mean(r_11)
factor = np.mean(r_10) / np.mean(r_11)
get_osa10_11_conversion_factor._log.info("Conversion factor is %f" % factor)
conversion_dict.update({'%.0f-%.0f' % (E1_isgri_keV, E2_isgri_keV): float(factor)})
import oda_integral_wrapper
with open('integral-osa-additional-parameters/osa11-10-conversion.json', 'w') as ff:
json.dump(conversion_dict, ff)
subprocess.call('cd integral-osa-additional-parameters;'+
'git config user.email "%s";'%(oda_integral_wrapper.__author_email__) +
'git config --global user.name "%s";' %(oda_integral_wrapper.__author__) +
'git commit osa11-10-conversion.json -m "Update factor %.0f-%.0f";'%(E1_isgri_keV, E2_isgri_keV)+
'git push ', shell=True)
subprocess.call('cd integral-osa-additional-parameters;' +
'git config user.email "%s";' % (oda_integral_wrapper.__author_email__) +
'git config --global user.name "%s";' % (oda_integral_wrapper.__author__) +
'git commit osa11-10-conversion.json -m "Update factor %.0f-%.0f";' % (
E1_isgri_keV, E2_isgri_keV) +
'git push ', shell=True)
return factor
oda_integral_wrapper/wrapper.py
View file @ c39cae5d
import numpy as np
from astropy.io import fits
import oda_api.api
import copy
import matplotlib
from astropy import units as u
from astroquery.simbad import Simbad
import astropy.wcs as wcs
import matplotlib
import matplotlib.pylab as plt
from astropy.coordinates import SkyCoord
from astropy import table
import numpy as np
import oda_api.api
import requests
from astropy import table
from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy.io import fits
from astroquery.simbad import Simbad
from autologging import logged
__author__ = "Carlo Ferrigno"
......@@ -98,12 +99,6 @@ class INTEGRALwrapper(object):
if n_poll == 0:
self.__log.info("At call %d we elaborate %d scw" % (n_call, len(scw_list)))
tot_num += len(scw_list)
# silent = False
# else:
# if n_poll % 10 == 0:
# silent = False
# else:
# silent = True
ys = "%06d" % n_call
......@@ -111,37 +106,44 @@ 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():
self.__log.debug('\n\n\nn_call ', n_call, '\nData ', data, '\nIs failed ', _disp.is_failed)
if data is None and not _disp.is_failed:
revs = get_revs(scw_list)
# Force 0SA10.2 for ISGRI before rev. 1626
self.__log.debug(f"min rev {revs.min()}")
if revs.min() < 1626 and 'isgri' in arguments['product']:
self.__log.debug("Forced OSA10.2 instead of OSA11")
self.__log.info("Forced OSA10.2 instead of OSA11")
local_arguments['osa_version'] = 'OSA10.2'
else:
local_arguments['osa_version'] = arguments['osa_version']
scw_list_str = ",".join([s for s in sorted(set(scw_list))])
self.__log.debug("Is submitted", _disp.is_submitted)
if not _disp.is_submitted:
revs = get_revs(scw_list)
# Force 0SA10.2 for ISGRI before rev. 1626
if revs.min() < 1626 and 'isgri' in arguments['product']:
local_arguments['osa_version'] = 'OSA10.2'
else:
local_arguments['osa_version'] = arguments['osa_version']
scw_list_str = ",".join([s for s in sorted(set(scw_list))])
data = _disp.get_product(scw_list=scw_list_str, **local_arguments)
else:
_disp.poll()
data_by_call[ys] = data
self.__log.debug("Is complete", _disp.is_complete)
if not _disp.is_complete:
continue
else:
data_by_call[ys] = _disp.get_product(scw_list=scw_list_str, **arguments)
data_by_call[ys] = _disp.get_product(scw_list=scw_list_str, **local_arguments)
self.__log.debug('n_poll ', n_poll)
n_complete = len([call for call, _disp in disp_by_call.items() if _disp.is_complete])
self.__log.debug('n_complete ', n_complete)
self.__log.debug(f"complete {n_complete} / {len(disp_by_call)}")
if n_complete == len(disp_by_call):
self.__log.info("done!")
break
self.__log.debug("not done")
n_poll += 1
#TODO test this option
if not wait:
return None
time.sleep(sleep_time)
......@@ -224,8 +226,6 @@ class INTEGRALwrapper(object):
idx_f, idx_s = find_duplicates(stacked_table)
# TODO tests with duplicates
if idx_f is None:
self.__log.info("No duplicates in final catalog")
else:
......@@ -371,7 +371,46 @@ class INTEGRALwrapper(object):
return combined_data
@staticmethod
def normalize_lc_ijd(combined_lc, source_name, Emin, Emax):
new_combined_lc = copy.deepcopy(combined_lc)
t, dt, y, dy = INTEGRALwrapper.get_lc(new_combined_lc, source_name)
ind = t < 5838.086367870370
from oda_integral_wrapper.get_osa10_11_factor import get_osa10_11_conversion_factor
f = get_osa10_11_conversion_factor(Emin, Emax)
y[ind] /= f
dy[ind] /= f
INTEGRALwrapper.put_lc(new_combined_lc, source_name, t, y, dy)
return new_combined_lc
@staticmethod
def put_lc(combined_lc, source_name, x, y, dy):
# In LC name has no "-" nor "+" ??????
patched_source_name = source_name.replace('-', ' ').replace('+', ' ')
hdu = None
j_index = -1
i_index = -1
for j, dd in enumerate(combined_lc._p_list):
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:
hdu = du.to_fits_hdu()
i_index = ii
j_index = j
if hdu is None:
INTEGRALwrapper.__log.info('Source ' + source_name + ' not found in lc')
return None
hdu.data['TIME'] = x
hdu.data['RATE'] = y
hdu.data['ERROR'] = dy
combined_lc._p_list[j_index].data_unit[i_index].from_fits_hdu(hdu)
return
@staticmethod
def get_lc(combined_lc, source_name, systematic_fraction=0):
......@@ -379,7 +418,6 @@ class INTEGRALwrapper(object):
# In LC name has no "-" nor "+" ??????
patched_source_name = source_name.replace('-', ' ').replace('+', ' ')
hdu = None
for j, dd in enumerate(combined_lc._p_list):
INTEGRALwrapper.__log.debug(dd.meta_data['src_name'])
......@@ -389,7 +427,7 @@ class INTEGRALwrapper(object):
hdu = du.to_fits_hdu()
if hdu is None:
INTEGRALwrapper.__log.info('Source ' + source_name + ' not found in lc, not plotting')
INTEGRALwrapper.__log.info('Source ' + source_name + ' not found in lc')
return None
x = hdu.data['TIME']
......@@ -407,7 +445,7 @@ class INTEGRALwrapper(object):
#This could only be valid for ISGRI
try:
dt_lc = hdu.data['XAX_E']
dt_lc = hdu.data['XAX_E']
INTEGRALwrapper.__log.debug('Get time bin directly from light curve')
except:
timedel = hdu.header['TIMEDEL']
......@@ -420,7 +458,6 @@ class INTEGRALwrapper(object):
return x[ind], dt_lc[ind], y[ind], dy[ind]
@staticmethod
def plot_lc(combined_lc, source_name, systematic_fraction=0, ng_sig_limit=3, find_excesses=False):
from scipy import stats
......
setup.py
View file @ c39cae5d
......@@ -15,7 +15,7 @@ include_package_data=True
scripts_list = glob.glob('./bin/*')
setup(name='oda_integral_wrapper',
version="1.3.13",
version="1.3.14",
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