STATION = 15  # gvdveen

import sapphire
import matplotlib.pyplot as plt
import numpy as np
import datetime as dt
%matplotlib inline

data = sapphire.quick_download(STATION)

# Download data van een andere dag
# start = dt.datetime(2019, 3, 20)
# end = dt.datetime(2019, 3, 21)
# sapphire.download_data(data, '/s%d' % STATION, start, end)

print(data)

node_naam = '/s%d/events' % STATION
node_naam

event_tabel = data.get_node(node_naam)
event_tabel

event_tabel[0]

event_tabel[1]

first_event = event_tabel[0]
first_event['timestamp']

event_tabel[0][6]  # 7de kolom van 1ste rij

first_event = event_tabel[0]
first_event['n1']

print(first_event['n1'])
print(event_tabel[0][6])

first_event['pulseheights']

print("pulshoogte detector 1: %d ADC (eerste event)" % first_event['pulseheights'][0])

events = event_tabel.read()

ts = events['timestamp']
ts[:30]

ns = events['nanoseconds']
ns[:30]

plt.figure(figsize=(10,4))
plt.hist(ns, histtype='step')
plt.ylabel('aantal')
plt.xlabel('nanoseconds deel van de timestamp')
plt.title('ESD Events van een enkele dag van station %d' % STATION)

plt.figure(figsize=(10,4))
plt.hist(ts, bins=24, histtype='step')
plt.ylabel('aantal')
plt.xlabel('timestamp [s]')
plt.title('ESD Events van een enkele dag van station %d' % STATION)

eerste_ts = ts[0]
eerste_ts

# linker en rechter grenzen van bins van 1 uur (3600 s) breed vanaf de eerste timestamp (1 dag)
bins = [eerste_ts + 3600 * h for h in range(25)]

plt.figure(figsize=(10,4))
plt.hist(ts, bins=bins, histtype='step')
plt.ylabel('aantal')
plt.xlabel('tijd vanaf middernacht [h]')
plt.xticks(bins, range(25))
plt.title('ESD Events van een enkele dag van station %d' % STATION)

n1 = event_tabel.col('n1')
plt.figure(figsize=(10,4))
plt.hist(n1, bins=np.arange(0.3, 5., .1), histtype='step')
plt.title('Station %d: Number of particles in detector 1' % STATION)
plt.xlabel('number of particles (MIP)')
plt.ylabel('counts')

ph = event_tabel.col('pulseheights')
ph1 = ph[:, 0]
ph2 = ph[:, 1]

plt.figure(figsize=(10,4))
plt.hist(ph1, bins=np.arange(0, 2000., 20.), histtype='step', log=True)
plt.hist(ph2, bins=np.arange(0, 2000., 20.), histtype='step', log=True)
plt.title('Station %d: Pulseheights' % STATION)
plt.xlabel('Pulseheight (ADC)')
plt.ylabel('counts')
plt.legend(['detector 1', 'detector 2' ])
plt.ylim(10, 1e4)