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)
Na enige tijd verschijnt hierboven een regel zoals:
"100%|############################################################|Time: 0:00:06
"
Soms is de download zo snel dat deze regel niet wordt afgedrukt.
De variabele "data
" bevat nu een set meetgegevens. Deze set is af te drukken.
# 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)
Het "data
" bestand heeft een hierarchise opbouw. In "data
" zit een
RootGroup, deze is te benaderen met "data.root
". Hierin zit weer een groep
"s102
", deze is te benaderen met "data.root.s102
". Hierin zit een tabel
"events
".
Voor het gemak maken we een variable
events
die naar de eventstabel van het station wijst:
De tabel heeft een
bepaalde plaats het HDF5 data bestand: /s????/events
waarbij ????
staat voor
het station nummer. Deze plaats heet een node
:
node_naam = '/s%d/events' % STATION
node_naam
event_tabel = data.get_node(node_naam)
event_tabel
Dit is een tabel tienduizenden regels. Elke regel is een event.
De informatie van het eerste event is op te halen met:
event_tabel[0]
Het tweede event: (Let op, python telt vanaf 0 en niet vanaf 1)
event_tabel[1]
De informatie in een event bestaat uit een lijst getallen. Deze getallen hebben de volgende betekenis:
-1
" betekent dat er geen detector was.-1
" betekent ook hier dat er geen detector was.In het
werkblad https://docs.hisparc.nl/infopakket/pdf/traces.pdf
wordt de natuurkundige betekenis van deze
getallen beschreven. De afbeeldingen in dit werkblad zijn afkomstig uit het
interactieve werkblad https://data.hisparc.nl/media/jsparc/jsparc.html.
Let op, computers tellen vanaf "0
" en niet vanaf "1
".
Een kolom zoals 'event_id', 'timestamp' of 't1' kan opgevraagd worden door de index van de kolom (0, 1, 2, ...) of door de kolomnaam. Door gebruik te maken van de kolomnaam wordt de code veel beter leesbaar:
first_event = event_tabel[0]
first_event['timestamp']
Het aantal gereconstrueerde deeltjes in detector 1 (het zevende getal) bij het eerste event is dus te vinden met:
event_tabel[0][6] # 7de kolom van 1ste rij
en:
first_event = event_tabel[0]
first_event['n1']
De tweede code is weliswaar langer, maar veel beter leesbaar.
print(first_event['n1'])
print(event_tabel[0][6])
Een array met pulshoogten in ADC-waarden is in dit geval te vinden met:
first_event['pulseheights']
Merk op dat de pulshoogtes van detector 3 en 4 de waarde '-1' hebben. De waarde '-1' betekent dat de pulsehoogte niet bepaald kon worden; Station 102 heeft slechts twee detectoren.
De eerste pulshoogte is te vinden met:
print("pulshoogte detector 1: %d ADC (eerste event)" % first_event['pulseheights'][0])
Vaak is het eenvoudiger om een hele kolom bijvoorbeeld timestamp
in een keer te bekijken.
Eerst lezen we de hele tabel in het geheugen. Het object events
is de gehele tabel:
events = event_tabel.read()
De variabele ts
wijst naar de kolom timestamp
en we bekijken de eerste 30 regels (events):
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')
Maak een histogram van de pulshoogtes van detector 1 en 2 van het station.
Een voorbeeld is hier te zien: https://data.hisparc.nl/show/stations/15
ph = event_tabel.col('pulseheights')
ph1 = ph[:, 0]
ph2 = ph[:, 1]
'pulseheights' is een matrix:
[:, 0]
is de gehele eerste rij, dwz de pulshoogtes per event van detector 0[:, 1]
is de gehele tweede rij, dwz de pulshoogtes per event van detector 1plt.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)