import numpy as np
import matplotlib.pyplot as plt
from Chandra.Time import DateTime
from astropy.table import Table, vstack
from mica.archive import aca_l0
%matplotlib inline
# from Ska.Matplotlib import plot_cxctime

from Ska.DBI import DBI

from Ska.Matplotlib import plot_cxctime

dbh = DBI(dbi='sybase', server='sybase', user='aca_read')

tstart= DateTime('2013:001').secs
tstop = DateTime('2014:180').secs

stars = dbh.fetchall('select obsid, slot, kalman_tstart, kalman_tstop, aoacmag_median, aoacmag_5th, aoacmag_95th, '
                     'n_samples, not_tracking_samples'
                     ' from trak_stats_data'
                     ' where obsid > 50000 and'
                     ' aoacmag_median > 10.2 and'
                     ' not_tracking_samples < 300 and'
                     ' n_samples > 3000 and'
                     ' kalman_tstart > {} AND kalman_tstart < {}'.format(tstart, tstop))

len(stars)

stars = Table(stars)

stars[:5]

def get_img0(star):
    files = aca_l0.get_files(start=star['kalman_tstart'], stop=star['kalman_tstop'], slots=[star['slot']], imgsize=[8])
    if not files:
        raise ValueError('No 8x8')
    tables = []
    for fn in files:
        tables.append(Table.read(fn))
    
    if len(tables) > 1:
        t = vstack(tables, metadata_conflicts='silent')
    else:
        t = tables[0]
    ok = (t['TIME'] > star['kalman_tstart']) & (t['TIME'] < star['kalman_tstop']) & (t['QUALITY'] == 0)
    t = t[ok]
    return t

def sum_imgraw(imr):
    for r0 in (0, 7):
        for c0 in (0, 1, 6, 7):
            imr[:, r0, c0] = 0

    return np.sum(imr.reshape(-1, 64), axis=1)

stars['p1'] = -1
stars['p5'] = -1
stars['p16'] = -1
stars['p50'] = -1
stars['p84'] = -1
stars['p95'] = -1
stars['p99'] = -1
stars['sigma'] = -1

for star in stars:
    obsid = star['obsid']
    print 'Processing obsid {}'.format(obsid),
    try:
        img0 = get_img0(star)
    except Exception as err:
        print('Error {}'.format(err))
        continue
    flux = sum_imgraw(img0['IMGRAW'])
    p1, p5, p16, p50, p84, p95, p99 = np.percentile(flux, [1, 5, 16, 50, 84, 95, 99])
    star['p1'] = p1
    star['p5'] = p5
    star['p16'] = p16
    star['p50'] = p50
    star['p84'] = p84
    star['p95'] = p95
    star['p99'] = p99
    star['sigma'] = (p84 - p16) * 5 / 1.7  # 5 e-/DN / (1.7 sec/readout)
    print(' sigma={:.1f} e-/sec'.format(star['sigma']))


stars[:10]

ok = stars['sigma'] > 0
sigma = stars['sigma'][ok]

plt.hist(sigma, bins=20)
plt.title('1-$\sigma$ noise in integrated counts)')
plt.xlabel('Noise (e-/sec)');

np.mean(sigma)