#!/usr/bin/env python __author__ = "Roman Krivonos" __copyright__ = "Space Research Institute (IKI)" """ Работает только в архиве """ import re import sys from glob import glob from astropy.io import fits from astropy.io.fits import getheader from astropy.coordinates import SkyCoord # High-level coordinates from astropy.coordinates import ICRS, Galactic, FK4, FK5 # Low-level frames from astropy.coordinates import Angle, Latitude, Longitude # Angles import astropy.units as u def file_grep(pattern=None, filename=None, filt=None): fl = open(filename, "r") for line in fl: if re.search(pattern, line): if(filt): if(filt in line): return(line.replace("\n", "")) else: return(line.replace("\n", "")) return(None) bands={ 'E01':'25.00 60.00', 'E02':'17.00 22.10', 'E03':'22.10 28.73', 'E04':'28.73 37.35', 'E05':'37.35 48.55', 'E06':'48.55 63.12', 'E07':'63.12 82.06', 'E08':'82.06 106.67', 'E09':'106.67 138.67', 'E10':'138.67 180.28', 'E11':'180.28 234.36', 'E12':'234.36 304.67', } crab_crd = SkyCoord(83.6287, 22.0147, frame=FK5(), unit="deg") key = sys.argv[1] do_mjd = False revs = sorted(glob("????", recursive = False)) erange = bands[key] fout="detcnts.{}.fits".format(key) obs_id=[] ra=[] dec=[] rota=[] lon=[] lat=[] t_start=[] t_stop=[] exposure=[] ph1=[] ph2=[] ph=[] raw=[] src=[] clean=[] npix=[] rms=[] orbit=[] tstart=[] tstop=[] mjd=[] crab_sep=[] for rev in revs: logs = sorted(glob(rev+"/????????????.log", recursive = False)) for log in logs: res = file_grep(pattern="rev_phase1",filename=log) phase1 = float(res.split()[1]) if(res) else 0.0 res = file_grep(pattern="rev_phase2",filename=log) phase2 = float(res.split()[1]) if(res) else 0.0 res = file_grep(pattern="Final FSD RMS",filename=log) rms0 = float(res.split()[-1]) if(res) else None res = file_grep(pattern="OBTMISDC",filename=log) obtmisdc = float(res.split()[-1]) if(res) else None res = file_grep(pattern="DETCNTS2",filename=log,filt=erange) d = res.split() if(res) else None if not (d): continue # ['DETCNTS2', '011000130010', '00001', '266.570', '-26.994', '89.359', '2639.38', '28148860.0', '28152310.0', '3.55128E-03', '8.53612E-04', '2.69767E-03', '1.27125E+04', '28.73', '37.35'] obsid=d[1] ra0=float(d[3]) dec0=float(d[4]) rota0=float(d[5]) texp=float(d[6]) t1=float(d[7]) t2=float(d[8]) raw0=float(d[9]) src0=float(d[10]) clean0=float(d[11]) npix0=float(d[12]) e1=float(d[13]) e2=float(d[14]) if not (raw0 > 0.0): continue obs_id.append(obsid) orbit.append(int(obsid[:4])) ra.append(ra0) dec.append(dec0) rota.append(rota0) sc = SkyCoord(ra0, dec0, frame=FK5(), unit="deg") lon.append(sc.galactic.l.value if(sc.galactic.l.value<180) else sc.galactic.l.value-360) lat.append(sc.galactic.b.value) crab_sep.append(crab_crd.separation(sc).deg) t_start.append(t1) t_stop.append(t2) exposure.append(texp) raw.append(raw0) src.append(src0) clean.append(clean0) npix.append(npix0) ph1.append(phase1) ph2.append(phase2) ph.append((phase1+phase2)/2.0) rms.append(rms0) if(do_mjd): # ~/results1/osa11/scw/0101/010100040010.001/isgri_cor_events.fits.gz corfile="/afs/mpa/project/integral/data1/scw/{}/{}.001/swg.fits".format(obsid[:4],obsid) d = getheader(corfile, 1) mjdref=float(d['MJDREF']) # tjd2/86400+obtmisdc+mjdref-40000 tj1=float(d['TSTART']) tj2=float(d['TSTOP']) tstart.append(tj1) tstop.append(tj2) mjd.append((tj1+tj2)/2.0+mjdref) #print(obtmisdc,mjdref) #print("{} {} {} {}".format(log,phase1,phase2,obsid)) else: mjd.append(0.0) indices = sorted( range(len(t_start)), key=lambda index: t_start[index] ) coldefs = fits.ColDefs([ fits.Column(name='REV', format='J', array=[orbit[index] for index in indices]), fits.Column(name='OBSID', format='12A', array=[obs_id[index] for index in indices]), fits.Column(name='RA', format='D', unit='deg', array=[ra[index] for index in indices]), fits.Column(name='DEC', format='D', unit='deg', array=[dec[index] for index in indices]), fits.Column(name='ROTA', format='D', unit='deg', array=[rota[index] for index in indices]), fits.Column(name='LON', format='D', unit='deg', array=[lon[index] for index in indices]), fits.Column(name='LAT', format='D', unit='deg', array=[lat[index] for index in indices]), #fits.Column(name='SC_TSTART', format='D', unit='s', array=[t_start[index] for index in indices]), #fits.Column(name='SC_TSTOP', format='D', unit='s', array=[t_stop[index] for index in indices]), #fits.Column(name='TSTART', format='D', unit='day', array=[tstart[index] for index in indices]), #fits.Column(name='TSTOP', format='D', unit='day', array=[tstop[index] for index in indices]), fits.Column(name='MJD', format='D', unit='day', array=[mjd[index] for index in indices]), fits.Column(name='EXPOSURE', format='D', unit='s', array=[exposure[index] for index in indices]), #fits.Column(name='PHASE1', format='D', unit='', array=[ph1[index] for index in indices]), #fits.Column(name='PHASE2', format='D', unit='', array=[ph2[index] for index in indices]), fits.Column(name='PHASE', format='D', unit='', array=[ph[index] for index in indices]), #fits.Column(name='NROWS', format='J', unit='', array=[nrows[index] for index in indices]), #fits.Column(name='LIVETIME', format='D', unit='s', array=[livetime[index] for index in indices]), fits.Column(name='RAW', format='D', unit='cts/s', array=[raw[index] for index in indices]), fits.Column(name='SRC', format='D', unit='', array=[src[index] for index in indices]), fits.Column(name='CLEAN', format='D', unit='', array=[clean[index] for index in indices]), fits.Column(name='RMS', format='D', unit='', array=[rms[index] for index in indices]), fits.Column(name='CRAB_SEP', format='D', unit='deg', array=[crab_sep[index] for index in indices]), ]) hdu = fits.BinTableHDU.from_columns(coldefs, name='RATE') hdu.header['MISSION'] = ('INTEGRAL', '') hdu.header['ERANGE'] = (key, '') hdu.header['E1'] = (erange.split()[0], '') hdu.header['E2'] = (erange.split()[1], '') hdu.header['TELESCOP'] = (key, '') hdu.header['INSTITUT'] = ('IKI', 'Affiliation') hdu.header['AUTHOR'] = ('Roman Krivonos', 'Responsible person') hdu.header['EMAIL'] = ('krivonos@cosmos.ru', 'E-mail') #hdu.add_checksum() print(hdu.columns) hdu.writeto(fout, overwrite=True) with fits.open(fout, mode='update') as hdus: hdus[1].add_checksum()