import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

FILENAME = 'events.h5'
station = 501

from datetime import datetime
start = datetime(2015, 5, 1)
end = datetime(2015, 5, 2)

import tables
data = tables.open_file(FILENAME, 'a')

from sapphire import download_data
if '/events' not in data:
    print('downloading events')
    download_data(data, '/', station, start, end)
else:
    print('events already downloaded.')

events = data.root.events
events

t1 = events.col('t1')
t3 = events.col('t3')

dt = t3 - t1
_ = plt.hist(dt, bins = 50, histtype='step')

ph = events.col('pulseheights')
ph1 = ph[:, 0]
ph3 = ph[:, 2]
dt_1_3 = dt.compress((t1 > 0.) & (t3 > 0.) & (ph1 > 200.) & (ph3 > 200.))
_ = plt.hist(dt_1_3, bins=50, histtype='step')

bins=np.arange(-100, 100, 5)
n, _, _ = plt.hist(dt_1_3, bins=bins, histtype='step')

np.mean(dt_1_3)

def gauss(x, N, mu, sigma):
    """Gaussian distribution"""
    return N * np.exp(-0.5 * ((x - mu) / sigma) ** 2)

from scipy.optimize import curve_fit
x = (bins[:-1] + bins[1:]) / 2
popt, pcov = curve_fit(gauss, x, n)

N, mu, sigma = popt
mu, sigma

plt.hist(dt_1_3, bins=bins, histtype='step')
plt.plot(bins, gauss(bins, N, mu, sigma), 'r--')
plt.legend(['fit','dt'])
plt.xlabel('t3 - t1 [ns]')
plt.ylabel('counts')
plt.xlim((-60, 70))

from sapphire import Station
s501 = Station(501)

from sapphire.transformations.clock import datetime_to_gps
ts = datetime_to_gps(datetime(2016,5,1))
ts

s501.detector_timing_offset(ts)

o1, o2, o3, o4 = s501.detector_timing_offset(ts)
print("De detector tijd offset tussen detector 3 en 1 op is: %2.1f ns (t = %d)" % ((o3 - o1), ts))

data.close()