Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision

Target

Select target project
  • km3py/km3io
  • spenamartinez/km3io
2 results
Select Git revision
Show changes
Commits on Source (43)
Showing with 1109 additions and 197 deletions
......@@ -54,6 +54,14 @@ test-py3.7:
- make test
<<: *junit_definition
test-py3.8:
image: docker.km3net.de/base/python:3.8
stage: test
script:
- *virtualenv_definition
- make test
<<: *junit_definition
code-style:
image: docker.km3net.de/base/python:3.7
stage: test
......
Unreleased changes
------------------
Version 0
---------
0.8.1 / 2020-02-10
~~~~~~~~~~~~~~~~~~
* update of reco data from offline files
* Documentation on how to read DAQ data
0.8.0 / 2020-01-23
~~~~~~~~~~~~~~~~~~
* Offline file headers are now accessible
0.7.0 / 2020-01-23
~~~~~~~~~~~~~~~~~~
* Reading of summary slice status information is now supported
0.6.3 / 2020-01-09
~~~~~~~~~~~~~~~~~~
* Bugfixes
0.6.2 / 2019-12-22
~~~~~~~~~~~~~~~~~~
* Fixes slicing of ``OfflineTracks``
......@@ -32,7 +48,7 @@ Version 0
0.3.0 / 2019-11-19
~~~~~~~~~~~~~~~~~~~
* Preliminary Jpp timeslice reader prototype
* Updated ``AaetReader``
* Updated ``AanetReader``
* Updated docs
0.2.1 / 2019-11-15
......
The km3io Python package
========================
.. image:: https://git.km3net.de/km3py/km3io/badges/master/build.svg
.. image:: https://git.km3net.de/km3py/km3io/badges/master/pipeline.svg
:target: https://git.km3net.de/km3py/km3io/pipelines
.. image:: https://git.km3net.de/km3py/km3io/badges/master/coverage.svg
......@@ -56,6 +56,12 @@ Tutorial
* `DAQ files reader <#daq-files-reader>`__
* `Reading Events <#reading-events>`__
* `Reading SummarySlices <#reading-summaryslices>`__
* `Reading Timeslices <#reading-timeslices>`__
* `Offline files reader <#offline-file-reader>`__
* `reading events data <#reading-events-data>`__
......@@ -92,41 +98,140 @@ A jagged array, is a 2+ dimentional array with different arrays lengths. In othe
# <JaggedArray [[102 102 102 ... 11517 11518 11518] [] [101 101 102 ... 11518 11518 11518] ... [101 101 102 ... 11516 11516 11517] [] [101 101 101 ... 11517 11517 11518]] at 0x7f74b0ef8810>
Overview of daq files
Overview of DAQ files
"""""""""""""""""""""
# info needed here
DAQ files, or also called online files, are written by the DataWriter and
contain events, timeslics and summary slices.
Overview of offline files
"""""""""""""""""""""""""
# info needed here
Offline files contain data about events, hits and tracks.
DAQ files reader
----------------
# an update is needed here?
Currently only events (the ``KM3NET_EVENT`` tree) are supported but timeslices and summaryslices will be implemented very soon.
``km3io`` is able to read events, summary slices and timeslices (except the L0
slices, which is work in progress).
Let's have a look at some ORCA data (``KM3NeT_00000044_00005404.root``)
Reading Events
""""""""""""""
To get a lazy ragged array of the events:
.. code-block:: python3
import km3io as ki
events = ki.DAQReader("KM3NeT_00000044_00005404.root").events
import km3io
f = km3io.DAQReader("KM3NeT_00000044_00005404.root")
That's it! Now let's have a look at the hits data:
That's it, we created an object which gives access to all the events, but the
relevant data is still not loaded into the memory (lazy access)!
Now let's have a look at the hits data:
.. code-block:: python3
>>> events
>>> f.events
Number of events: 17023
>>> events[23].snapshot_hits.tot
>>> f.events[23].snapshot_hits.tot
array([28, 22, 17, 29, 5, 27, 24, 26, 21, 28, 26, 21, 26, 24, 17, 28, 23,29, 27, 24, 23, 26, 29, 25, 18, 28, 24, 28, 26, 20, 25, 31, 28, 23, 26, 21, 30, 33, 27, 16, 23, 24, 19, 24, 27, 22, 23, 21, 25, 16, 28, 22, 22, 29, 24, 29, 24, 24, 25, 25, 21, 31, 26, 28, 30, 42, 28], dtype=uint8)
The resulting arrays are numpy arrays.
Reading SummarySlices
"""""""""""""""""""""
The following example shows how to access summary slices, in particular the DOM
IDs of the slice with the index ``23``:
.. code-block:: python3
>>> f.summaryslices
<km3io.daq.SummmarySlices at 0x7effcc0e52b0>
>>> f.summaryslices.slices[23].dom_id
array([806451572, 806455814, 806465101, 806483369, 806487219, 806487226,
806487231, 808432835, 808435278, 808447180, 808447186, 808451904,
808451907, 808469129, 808472260, 808472265, 808488895, 808488990,
808489014, 808489117, 808493910, 808946818, 808949744, 808951460,
808956908, 808959411, 808961448, 808961480, 808961504, 808961655,
808964815, 808964852, 808964883, 808964908, 808969848, 808969857,
808972593, 808972598, 808972698, 808974758, 808974773, 808974811,
808974972, 808976377, 808979567, 808979721, 808979729, 808981510,
808981523, 808981672, 808981812, 808981864, 808982005, 808982018,
808982041, 808982066, 808982077, 808982547, 808984711, 808996773,
808997793, 809006037, 809007627, 809503416, 809521500, 809524432,
809526097, 809544058, 809544061], dtype=int32)
The ``.dtype`` attribute (or in general, <TAB> completion) is useful to find out
more about the field structure:
.. code-block:: python3
>>> f.summaryslices.headers.dtype
dtype([(' cnt', '<u4'), (' vers', '<u2'), (' cnt2', '<u4'), (' vers2',
'<u2'), (' cnt3', '<u4'), (' vers3', '<u2'), ('detector_id', '<i4'), ('run',
'<i4'), ('frame_index', '<i4'), (' cnt4', '<u4'), (' vers4', '<u2'),
('UTC_seconds', '<u4'), ('UTC_16nanosecondcycles', '<u4')])
>>> f.summaryslices.headers.frame_index
<ChunkedArray [162 163 173 ... 36001 36002 36003] at 0x7effccd4af10>
The resulting array is a ``ChunkedArray`` which is an extended version of a
numpy array and behaves like one.
Reading Timeslices
""""""""""""""""""
Timeslices are split into different streams since 2017 and ``km3io`` currently
supports everything except L0, i.e. L1, L2 and SN streams. The API is
work-in-progress and will be improved in future, however, all the data is
already accessible (although in ugly ways ;-)
To access the timeslice data:
.. code-block:: python3
>>> f.timeslices
Available timeslice streams: L1, SN
>>> f.timeslices.stream("L1", 24).frames
{806451572: <Table [<Row 1577843> <Row 1577844> ... <Row 1578147>],
806455814: <Table [<Row 1578148> <Row 1578149> ... <Row 1579446>],
806465101: <Table [<Row 1579447> <Row 1579448> ... <Row 1580885>],
...
}
The frames are represented by a dictionary where the key is the ``DOM ID`` and
the value a numpy array of hits, with the usual fields to access the PMT
channel, time and ToT:
.. code-block:: python3
>>> f.timeslices.stream("L1", 24).frames[806451572].dtype
dtype([('pmt', 'u1'), ('tdc', '<u4'), ('tot', 'u1')])
>>> f.timeslices.stream("L1", 24).frames[806451572].tot
array([29, 21, 8, 29, 22, 20, 1, 37, 11, 22, 11, 22, 12, 20, 29, 94, 26,
26, 18, 16, 13, 22, 6, 29, 24, 30, 14, 26, 12, 23, 4, 25, 6, 27,
5, 13, 21, 28, 30, 4, 25, 10, 5, 6, 5, 17, 4, 27, 24, 25, 27,
28, 32, 6, 3, 15, 3, 20, 33, 30, 30, 20, 28, 6, 7, 3, 14, 12,
25, 27, 26, 25, 22, 21, 23, 6, 20, 21, 4, 4, 10, 24, 29, 12, 30,
5, 3, 24, 15, 14, 25, 5, 27, 23, 26, 4, 28, 15, 34, 22, 4, 29,
24, 26, 29, 23, 25, 28, 14, 31, 27, 26, 27, 28, 23, 54, 4, 25, 11,
28, 25, 24, 7, 27, 28, 28, 18, 3, 13, 14, 38, 28, 4, 21, 16, 16,
4, 21, 26, 21, 28, 64, 21, 1, 24, 21, 26, 26, 25, 4, 28, 11, 31,
10, 24, 24, 28, 10, 6, 4, 20, 26, 18, 5, 18, 24, 5, 27, 23, 20,
29, 20, 6, 18, 5, 24, 17, 28, 24, 15, 26, 27, 25, 9, 3, 18, 3,
34, 29, 10, 25, 30, 28, 19, 26, 34, 27, 14, 17, 15, 26, 8, 19, 5,
27, 13, 5, 27, 46, 3, 25, 13, 30, 9, 21, 12, 1, 32, 25, 8, 30,
4, 24, 11, 3, 11, 27, 5, 13, 5, 16, 18, 3, 22, 10, 7, 32, 29,
15, 20, 18, 16, 27, 5, 22, 4, 33, 5, 29, 24, 30, 7, 7, 25, 33,
7, 20, 8, 30, 4, 4, 6, 26, 8, 24, 22, 12, 6, 3, 21, 28, 11,
24, 27, 27, 6, 29, 5, 18, 11, 26, 5, 19, 32, 25, 4, 20, 35, 30,
5, 3, 26, 30, 23, 28, 6, 25, 25, 5, 45, 23, 18, 29, 28, 23],
dtype=uint8)
Offline files reader
--------------------
......@@ -140,12 +245,59 @@ First, let's read our file:
.. code-block:: python3
>>> import km3io as ki
>>> file = 'datav6.0test.jchain.aanet.00005971.root'
>>> file = 'my_file.root'
>>> r = ki.OfflineReader(file)
<km3io.aanet.OfflineReader at 0x7f24cc2bd550>
<km3io.offline.OfflineReader at 0x7f24cc2bd550>
and that's it! Note that `file` can be either an str of your file path, or a path-like object.
To read the file header:
.. code-block:: python3
>>> r.header
DAQ 394
PDF 4 58
XSecFile
can 0 1027 888.4
can_user 0.00 1027.00 888.40
coord_origin 0 0 0
cut_in 0 0 0 0
cut_nu 100 1e+08 -1 1
cut_primary 0 0 0 0
cut_seamuon 0 0 0 0
decay doesnt happen
detector NOT
drawing Volume
end_event
genhencut 2000 0
genvol 0 1027 888.4 2.649e+09 100000
kcut 2
livetime 0 0
model 1 2 0 1 12
muon_desc_file
ngen 0.1000E+06
norma 0 0
nuflux 0 3 0 0.500E+00 0.000E+00 0.100E+01 0.300E+01
physics GENHEN 7.2-220514 181116 1138
seed GENHEN 3 305765867 0 0
simul JSirene 11012 11/17/18 07
sourcemode diffuse
spectrum -1.4
start_run 1
target isoscalar
usedetfile false
xlat_user 0.63297
xparam OFF
zed_user 0.00 3450.00
**Note:** not all file header types are supported, so don't be surprised when you get the following warning
.. code-block:: python3
/home/zineb/km3net/km3net/km3io/km3io/offline.py:341: UserWarning: Your file header has an unsupported format
warnings.warn("Your file header has an unsupported format")
To explore all the available branches in our offline file:
.. code-block:: python3
......@@ -458,7 +610,7 @@ to get all data of a specific hit (let's say hit 0) from event 0:
.. code-block:: python3
>>>r[0].hits[0]
>>> r[0].hits[0]
offline hit:
id : 0
dom_id : 806451572
......@@ -485,7 +637,7 @@ to get a specific value from hit 0 in event 0, let's say for example the dom id:
.. code-block:: python3
>>>r[0].hits[0].dom_id
>>> r[0].hits[0].dom_id
806451572
reading tracks data
......@@ -550,7 +702,7 @@ to get all data of a specific track (let's say track 0) from event 0:
.. code-block:: python3
>>>r[0].tracks[0]
>>> r[0].tracks[0]
offline track:
fUniqueID : 0
fBits : 33554432
......@@ -595,9 +747,93 @@ to get a specific value from track 0 in event 0, let's say for example the likli
.. code-block:: python3
>>>r[0].tracks[0].lik
>>> r[0].tracks[0].lik
294.6407542676734
to get the reconstruction parameters, first take a look at the available reconstruction keys:
.. code-block:: python3
>>> r.best_reco.dtype.names
['JGANDALF_BETA0_RAD',
'JGANDALF_BETA1_RAD',
'JGANDALF_CHI2',
'JGANDALF_NUMBER_OF_HITS',
'JENERGY_ENERGY',
'JENERGY_CHI2',
'JGANDALF_LAMBDA',
'JGANDALF_NUMBER_OF_ITERATIONS',
'JSTART_NPE_MIP',
'JSTART_NPE_MIP_TOTAL',
'JSTART_LENGTH_METRES',
'JVETO_NPE',
'JVETO_NUMBER_OF_HITS',
'JENERGY_MUON_RANGE_METRES',
'JENERGY_NOISE_LIKELIHOOD',
'JENERGY_NDF',
'JENERGY_NUMBER_OF_HITS']
the keys above can also be accessed with a tab completion:
.. image:: https://git.km3net.de/km3py/km3io/raw/master/examples/pictures/reco.png
to get a numpy `recarray <https://docs.scipy.org/doc/numpy/reference/generated/numpy.recarray.html>`__ of all fit data of the best reconstructed track:
.. code-block:: python3
>>> r.best_reco
to get an array of a parameter of interest, let's say `'JENERGY_ENERGY'`:
.. code-block:: python3
>>> r.best_reco['JENERGY_ENERGY']
array([1141.87137899, 4708.16378575, 499.7243005 , 103.54680875,
208.6103912 , 1336.52338666, 998.87632267, 1206.54345674,
16.28973662])
**Note**: In km3io, the best fit is defined as the track fit with the maximum reconstruction stages. When "nan" is returned, it means that the reconstruction parameter of interest is not found. for example, in the case of muon simulations: if `[1, 2]` are the reconstruction stages, then only the fit parameters corresponding to the stages `[1, 2]` are found in the Offline files, the remaining fit parameters corresponding to the stages `[3, 4, 5]` are all filled with nan.
to get a numpy recarray of the fit data of tracks with specific reconstruction stages, let's say `[1, 2, 3, 4, 5]` in the case of a muon track reconstruction:
.. code-block:: python3
>>> r.get_reco_fit([1, 2, 3, 4, 5])
again, to get the reconstruction parameters names:
.. code-block:: python3
>>> r.get_reco_fit([1, 2, 3, 4, 5]).dtype.names
('JGANDALF_BETA0_RAD',
'JGANDALF_BETA1_RAD',
'JGANDALF_CHI2',
'JGANDALF_NUMBER_OF_HITS',
'JENERGY_ENERGY',
'JENERGY_CHI2',
'JGANDALF_LAMBDA',
'JGANDALF_NUMBER_OF_ITERATIONS',
'JSTART_NPE_MIP',
'JSTART_NPE_MIP_TOTAL',
'JSTART_LENGTH_METRES',
'JVETO_NPE',
'JVETO_NUMBER_OF_HITS',
'JENERGY_MUON_RANGE_METRES',
'JENERGY_NOISE_LIKELIHOOD',
'JENERGY_NDF',
'JENERGY_NUMBER_OF_HITS')
to get the reconstruction data of interest, for example ['JENERGY_ENERGY']:
.. code-block:: python3
>>> r.get_reco_fit([1, 2, 3, 4, 5])['JENERGY_ENERGY']
array([1141.87137899, 4708.16378575, 499.7243005 , 103.54680875,
208.6103912 , 1336.52338666, 998.87632267, 1206.54345674,
16.28973662])
**Note**: When the reconstruction stages of interest are not found in all your data file, an error is raised.
reading mc hits data
""""""""""""""""""""
......@@ -606,7 +842,7 @@ to read mc hits data:
.. code-block:: python3
>>>r.mc_hits
>>> r.mc_hits
<OfflineHits: 10 parsed elements>
that's it! All branches in mc hits tree can be accessed in the exact same way described in the section `reading hits data <#reading-hits-data>`__ . All data is easily accesible and if you are stuck, hit tab key to see all the available branches:
......@@ -620,7 +856,7 @@ to read mc tracks data:
.. code-block:: python3
>>>r.mc_tracks
>>> r.mc_tracks
<OfflineTracks: 10 parsed elements>
that's it! All branches in mc tracks tree can be accessed in the exact same way described in the section `reading tracks data <#reading-tracks-data>`__ . All data is easily accesible and if you are stuck, hit tab key to see all the available branches:
......
......@@ -26,7 +26,7 @@ release = get_distribution('km3io').version
version = '.'.join(release.split('.')[:2])
project = 'km3io {}'.format(km3io.__version__)
copyright = '{0}, Zineb Aly and Tamas Gal'.format(date.today().year)
author = 'Zineb Aly, Tamas Gal'
author = 'Zineb Aly, Tamas Gal, Johannes Schumann'
# -- General configuration ---------------------------------------------------
......
examples/pictures/reco.png

15.1 KiB

......@@ -2,7 +2,7 @@ import uproot
import numpy as np
import numba as nb
TIMESLICE_FRAME_BASKET_CACHE_SIZE = 23 * 1024**2 # [byte]
TIMESLICE_FRAME_BASKET_CACHE_SIZE = 523 * 1024**2 # [byte]
SUMMARYSLICE_FRAME_BASKET_CACHE_SIZE = 523 * 1024**2 # [byte]
# Parameters for PMT rate conversions, since the rates in summary slices are
......@@ -13,6 +13,8 @@ MINIMAL_RATE_HZ = 2.0e3
MAXIMAL_RATE_HZ = 2.0e6
RATE_FACTOR = np.log(MAXIMAL_RATE_HZ / MINIMAL_RATE_HZ) / 255
CHANNEL_BITS_TEMPLATE = np.zeros(31, dtype=bool)
@nb.vectorize([
nb.int32(nb.int8),
......@@ -28,10 +30,83 @@ def get_rate(value):
return MINIMAL_RATE_HZ * np.exp(value * RATE_FACTOR)
@nb.guvectorize("void(i8, b1[:], b1[:])",
"(), (n) -> (n)",
target="parallel",
nopython=True)
def unpack_bits(value, bits_template, out):
"""Return a boolean array for a value's bit representation.
This function also accepts arrays as input, the output shape will be
NxM where N is the number of input values and M the length of the
``bits_template`` array, which is just a dummy array, due to the weird
signature system of numba.
Parameters
----------
value: int or np.array(int) with shape (N,)
The binary value of containing the bit information
bits_template: np.array() with shape (M,)
The template for the output array, the only important is its shape
Returns
-------
np.array(bool) either with shape (M,) or (N, M)
"""
for i in range(bits_template.shape[0]):
out[30 - i] = value & (1 << i) > 0
def get_channel_flags(value):
"""Returns the hrv/fifo flags for the PMT channels (hrv/fifo)
Parameters
----------
value : int32
The integer value to be parsed.
"""
channel_bits = np.bitwise_and(value, 0x3FFFFFFF)
flags = unpack_bits(channel_bits, CHANNEL_BITS_TEMPLATE)
return np.flip(flags, axis=-1)
def get_number_udp_packets(value):
"""Returns the number of received UDP packets (dq_status)
Parameters
----------
value : int32
The integer value to be parsed.
"""
return np.bitwise_and(value, 0x7FFF)
def get_udp_max_sequence_number(value):
"""Returns the maximum sequence number of the received UDP packets (dq_status)
Parameters
----------
value : int32
The integer value to be parsed.
"""
return np.right_shift(value, 16)
def has_udp_trailer(value):
"""Returns the UDP Trailer flag (fifo)
Parameters
----------
value : int32
The integer value to be parsed.
"""
return np.any(np.bitwise_and(value, np.left_shift(1, 31)))
class DAQReader:
"""Reader for DAQ ROOT files"""
def __init__(self, filename):
self.fobj = uproot.open(filename)
self._fobj = uproot.open(filename)
self._events = None
self._timeslices = None
self._summaryslices = None
......@@ -39,7 +114,7 @@ class DAQReader:
@property
def events(self):
if self._events is None:
tree = self.fobj["KM3NET_EVENT"]
tree = self._fobj["KM3NET_EVENT"]
headers = tree["KM3NETDAQ::JDAQEventHeader"].array(
uproot.interpret(tree["KM3NETDAQ::JDAQEventHeader"],
......@@ -62,20 +137,20 @@ class DAQReader:
@property
def timeslices(self):
if self._timeslices is None:
self._timeslices = DAQTimeslices(self.fobj)
self._timeslices = DAQTimeslices(self._fobj)
return self._timeslices
@property
def summaryslices(self):
if self._summaryslices is None:
self._summaryslices = SummmarySlices(self.fobj)
self._summaryslices = SummmarySlices(self._fobj)
return self._summaryslices
class SummmarySlices:
"""A wrapper for summary slices"""
def __init__(self, fobj):
self.fobj = fobj
self._fobj = fobj
self._slices = None
self._headers = None
self._rates = None
......@@ -101,7 +176,7 @@ class SummmarySlices:
def _read_summaryslices(self):
"""Reads a lazyarray of summary slices"""
tree = self.fobj[b'KM3NET_SUMMARYSLICE'][b'KM3NET_SUMMARYSLICE']
tree = self._fobj[b'KM3NET_SUMMARYSLICE'][b'KM3NET_SUMMARYSLICE']
return tree[b'vector<KM3NETDAQ::JDAQSummaryFrame>'].lazyarray(
uproot.asjagged(uproot.astable(
uproot.asdtype([("dom_id", "i4"), ("dq_status", "u4"),
......@@ -114,7 +189,7 @@ class SummmarySlices:
def _read_headers(self):
"""Reads a lazyarray of summary slice headers"""
tree = self.fobj[b'KM3NET_SUMMARYSLICE'][b'KM3NET_SUMMARYSLICE']
tree = self._fobj[b'KM3NET_SUMMARYSLICE'][b'KM3NET_SUMMARYSLICE']
return tree[b'KM3NETDAQ::JDAQSummarysliceHeader'].lazyarray(
uproot.interpret(tree[b'KM3NETDAQ::JDAQSummarysliceHeader'],
cntvers=True))
......@@ -123,42 +198,35 @@ class SummmarySlices:
class DAQTimeslices:
"""A simple wrapper for DAQ timeslices"""
def __init__(self, fobj):
self.fobj = fobj
self._fobj = fobj
self._timeslices = {}
self._read_default_stream()
self._read_streams()
def _read_default_stream(self):
"""Read the default KM3NET_TIMESLICE stream"""
tree = self.fobj[b'KM3NET_TIMESLICE'][b'KM3NET_TIMESLICE']
headers = tree[b'KM3NETDAQ::JDAQTimesliceHeader']
superframes = tree[b'vector<KM3NETDAQ::JDAQSuperFrame>']
self._timeslices['default'] = (headers, superframes)
def _read_streams(self):
"""Read the L0, L1, L2 and SN streams if available"""
streams = [
streams = set(
s.split(b"KM3NET_TIMESLICE_")[1].split(b';')[0]
for s in self.fobj.keys() if b"KM3NET_TIMESLICE_" in s
]
for s in self._fobj.keys() if b"KM3NET_TIMESLICE_" in s)
for stream in streams:
tree = self.fobj[b'KM3NET_TIMESLICE_' +
stream][b'KM3NETDAQ::JDAQTimeslice']
tree = self._fobj[b'KM3NET_TIMESLICE_' +
stream][b'KM3NETDAQ::JDAQTimeslice']
headers = tree[b'KM3NETDAQ::JDAQTimesliceHeader'][
b'KM3NETDAQ::JDAQHeader'][b'KM3NETDAQ::JDAQChronometer']
if len(headers) == 0:
continue
superframes = tree[b'vector<KM3NETDAQ::JDAQSuperFrame>']
hits_dtype = np.dtype([("pmt", "u1"), ("tdc", "<u4"),
("tot", "u1")])
hits_buffer = superframes[
b'vector<KM3NETDAQ::JDAQSuperFrame>.buffer'].lazyarray(
uproot.asjagged(uproot.astable(
uproot.asdtype([("pmt", "u1"), ("tdc", "u4"),
("tot", "u1")])),
uproot.asjagged(uproot.astable(uproot.asdtype(hits_dtype)),
skipbytes=6),
basketcache=uproot.cache.ThreadSafeArrayCache(
TIMESLICE_FRAME_BASKET_CACHE_SIZE))
self._timeslices[stream.decode("ascii")] = (headers, superframes,
hits_buffer)
setattr(self, stream.decode("ascii"),
DAQTimesliceStream(headers, superframes, hits_buffer))
def stream(self, stream, idx):
ts = self._timeslices[stream]
......@@ -172,6 +240,21 @@ class DAQTimeslices:
return str(self)
class DAQTimesliceStream:
def __init__(self, headers, superframes, hits_buffer):
# self.headers = headers.lazyarray(
# uproot.asjagged(uproot.astable(
# uproot.asdtype(
# np.dtype([('a', 'i4'), ('b', 'i4'), ('c', 'i4'),
# ('d', 'i4'), ('e', 'i4')]))),
# skipbytes=6),
# basketcache=uproot.cache.ThreadSafeArrayCache(
# TIMESLICE_FRAME_BASKET_CACHE_SIZE))
self.headers = headers
self.superframes = superframes
self._hits_buffer = hits_buffer
class DAQTimeslice:
"""A wrapper for a DAQ timeslice"""
def __init__(self, header, superframe, hits_buffer, idx, stream):
......
# -*- coding: utf-8 -*-
"""
KM3NeT Data Definitions v1.1.2
https://git.km3net.de/common/km3net-dataformat
"""
# fitparameters
data = {
"JGANDALF_BETA0_RAD": 0,
"JGANDALF_BETA1_RAD": 1,
"JGANDALF_CHI2": 2,
"JGANDALF_NUMBER_OF_HITS": 3,
"JENERGY_ENERGY": 4,
"JENERGY_CHI2": 5,
"JGANDALF_LAMBDA": 6,
"JGANDALF_NUMBER_OF_ITERATIONS": 7,
"JSTART_NPE_MIP": 8,
"JSTART_NPE_MIP_TOTAL": 9,
"JSTART_LENGTH_METRES": 10,
"JVETO_NPE": 11,
"JVETO_NUMBER_OF_HITS": 12,
"JENERGY_MUON_RANGE_METRES": 13,
"JENERGY_NOISE_LIKELIHOOD": 14,
"JENERGY_NDF": 15,
"JENERGY_NUMBER_OF_HITS": 16,
"JCOPY_Z_M": 17,
}
# -*- coding: utf-8 -*-
"""
KM3NeT Data Definitions v1.1.2
https://git.km3net.de/common/km3net-dataformat
"""
# reconstruction
data = {
"JPP_RECONSTRUCTION_TYPE": 4000,
"JMUONFIT": 0,
"JMUONBEGIN": 0,
"JMUONPREFIT": 1,
"JMUONSIMPLEX": 2,
"JMUONGANDALF": 3,
"JMUONENERGY": 4,
"JMUONSTART": 5,
"JLINEFIT": 6,
"JMUONEND": 99,
"JSHOWERFIT": 100,
"JSHOWERBEGIN": 100,
"JSHOWERPREFIT": 101,
"JSHOWERPOSITIONFIT": 102,
"JSHOWERCOMPLETEFIT": 103,
"JSHOWER_BJORKEN_Y": 104,
"JSHOWEREND": 199,
"DUSJSHOWERFIT": 200,
"DUSJBEGIN": 200,
"DUSJPREFIT": 201,
"DUSJPOSITIONFIT": 202,
"JDUSJCOMPLETEFIT": 203,
"DUSJEND": 299,
"AASHOWERFIT": 300,
"AASHOWERBEGIN": 300,
"AASHOWERCOMPLETEFIT": 301,
"AASHOWEREND": 399,
"JUSERBEGIN": 1000,
"JMUONVETO": 1001,
"JMUONPATH": 1003,
"JMCEVT": 1004,
"JUSEREND": 1099,
"RECTYPE_UNKNOWN": -1,
"RECSTAGE_UNKNOWN": -1,
}
# -*- coding: utf-8 -*-
"""
KM3NeT Data Definitions v1.1.2
https://git.km3net.de/common/km3net-dataformat
"""
# trigger
data = {
"JTRIGGER3DSHOWER": 1,
"JTRIGGERMXSHOWER": 2,
"JTRIGGER3DMUON": 4,
"JTRIGGERNB": 5,
}
# -*- coding: utf-8 -*-
"""
KM3NeT Data Definitions v1.0.1
https://git.km3net.de/common/data-format
"""
# fitparameters
JGANDALF_BETA0_RAD = 0
JGANDALF_BETA1_RAD = 1
JGANDALF_CHI2 = 2
JGANDALF_NUMBER_OF_HITS = 3
JENERGY_ENERGY = 4
JENERGY_CHI2 = 5
JGANDALF_LAMBDA = 6
JGANDALF_NUMBER_OF_ITERATIONS = 7
JSTART_NPE_MIP = 8
JSTART_NPE_MIP_TOTAL = 9
JSTART_LENGTH_METRES = 10
JVETO_NPE = 11
JVETO_NUMBER_OF_HITS = 12
JENERGY_MUON_RANGE_METRES = 13
JENERGY_NOISE_LIKELIHOOD = 14
JENERGY_NDF = 15
JENERGY_NUMBER_OF_HITS = 16
JCOPY_Z_M = 17
import uproot
import numpy as np
import warnings
import km3io.definitions.trigger
import km3io.definitions.fitparameters
import km3io.definitions.reconstruction
# 110 MB based on the size of the largest basket found so far in km3net
BASKET_CACHE_SIZE = 110 * 1024**2
......@@ -27,6 +32,9 @@ class OfflineKeys:
self._cut_hits_keys = None
self._cut_tracks_keys = None
self._cut_events_keys = None
self._trigger = None
self._fitparameters = None
self._reconstruction = None
def __str__(self):
return '\n'.join([
......@@ -169,17 +177,10 @@ class OfflineKeys:
list of all "trks.fitinf" keys.
"""
if self._fit_keys is None:
# these are hardcoded because they are not outsourced in offline
# files
self._fit_keys = [
'JGANDALF_BETA0_RAD', 'JGANDALF_BETA1_RAD', 'JGANDALF_CHI2',
'JGANDALF_NUMBER_OF_HITS', 'JENERGY_ENERGY', 'JENERGY_CHI2',
'JGANDALF_LAMBDA', 'JGANDALF_NUMBER_OF_ITERATIONS',
'JSTART_NPE_MIP', 'JSTART_NPE_MIP_TOTAL',
'JSTART_LENGTH_METRES', 'JVETO_NPE', 'JVETO_NUMBER_OF_HITS',
'JENERGY_MUON_RANGE_METRES', 'JENERGY_NOISE_LIKELIHOOD',
'JENERGY_NDF', 'JENERGY_NUMBER_OF_HITS', 'JCOPY_Z_M'
]
self._fit_keys = sorted(self.fitparameters,
key=self.fitparameters.get,
reverse=False)
# self._fit_keys = [*fit.keys()]
return self._fit_keys
@property
......@@ -227,6 +228,48 @@ class OfflineKeys:
]
return self._cut_events_keys
@property
def trigger(self):
"""trigger parameters and their index from km3net-Dataformat.
Returns
-------
dict
dictionary of trigger parameters and their index in an Offline
file.
"""
if self._trigger is None:
self._trigger = km3io.definitions.trigger.data
return self._trigger
@property
def reconstruction(self):
"""reconstruction parameters and their index from km3net-Dataformat.
Returns
-------
dict
dictionary of reconstruction parameters and their index in an
Offline file.
"""
if self._reconstruction is None:
self._reconstruction = km3io.definitions.reconstruction.data
return self._reconstruction
@property
def fitparameters(self):
"""fit parameters parameters and their index from km3net-Dataformat.
Returns
-------
dict
dictionary of fit parameters and their index in an Offline
file.
"""
if self._fitparameters is None:
self._fitparameters = km3io.definitions.fitparameters.data
return self._fitparameters
class Reader:
"""Reader for one offline ROOT file"""
......@@ -319,6 +362,8 @@ class OfflineReader:
self._mc_hits = None
self._mc_tracks = None
self._keys = None
self._best_reco = None
self._header = None
def __getitem__(self, item):
return OfflineReader(file_path=self._file_path, data=self._data[item])
......@@ -326,6 +371,20 @@ class OfflineReader:
def __len__(self):
return len(self._data)
@property
def header(self):
if self._header is None:
fobj = uproot.open(self._file_path)
if b'Head;1' in fobj.keys():
self._header = {}
for n, x in fobj['Head']._map_3c_string_2c_string_3e_.items():
print("{:15s} {}".format(n.decode("utf-8"),
x.decode("utf-8")))
self._header[n.decode("utf-8")] = x.decode("utf-8")
if b'Header;1' in fobj.keys():
warnings.warn("Your file header has an unsupported format")
return self._header
@property
def keys(self):
"""wrapper for all keys in an offline file.
......@@ -382,7 +441,7 @@ class OfflineReader:
self._tracks = OfflineTracks(
self.keys.cut_tracks_keys,
[self._data[key] for key in self.keys.tracks_keys],
fit_keys=self.keys.fit_keys)
fitparameters=self.keys.fitparameters)
return self._tracks
@property
......@@ -412,9 +471,171 @@ class OfflineReader:
if self._mc_tracks is None:
self._mc_tracks = OfflineTracks(
self.keys.cut_tracks_keys,
[self._data[key] for key in self.keys.mc_tracks_keys])
[self._data[key] for key in self.keys.mc_tracks_keys],
fitparameters=self.keys.fitparameters)
return self._mc_tracks
@property
def best_reco(self):
"""returns the best reconstructed track fit data. The best fit is defined
as the track fit with the maximum reconstruction stages. When "nan" is
returned, it means that the reconstruction parameter of interest is not
found. for example, in the case of muon simulations: if [1, 2] are the
reconstruction stages, then only the fit parameters corresponding to the
stages [1, 2] are found in the Offline files, the remaining fit parameters
corresponding to the stages 3, 4, 5 are all filled with nan.
Returns
-------
numpy recarray
a recarray of the best track fit data (reconstruction data).
"""
if self._best_reco is None:
keys = ", ".join(self.keys.fit_keys[:-1])
empty_fit_info = np.array(
[match for match in self._find_empty(self.tracks.fitinf)])
fit_info = [
i for i, j in zip(self.tracks.fitinf, empty_fit_info[:, 1])
if j is not None
]
stages = self._get_max_reco_stages(self.tracks.rec_stages)
fit_data = np.array([i[j] for i, j in zip(fit_info, stages[:, 2])])
rows_size = len(max(fit_data, key=len))
equal_size_data = np.vstack([
np.hstack([i, np.zeros(rows_size - len(i)) + np.nan])
for i in fit_data
])
self._best_reco = np.core.records.fromarrays(
equal_size_data.transpose(), names=keys)
return self._best_reco
def _get_max_reco_stages(self, reco_stages):
"""find the longest reconstructed track based on the maximum size of
reconstructed stages.
Parameters
----------
reco_stages : chunked array
chunked array of all the reconstruction stages of all tracks.
In km3io, it is accessed with
km3io.OfflineReader(my_file).tracks.rec_stages .
Returns
-------
numpy array
array with 3 columns: *list of the maximum reco_stages
*lentgh of the maximum reco_stages
*position of the maximum reco_stages
"""
empty_reco_stages = np.array(
[match for match in self._find_empty(reco_stages)])
max_reco_stages = np.array(
[[max(i, key=len),
len(max(i, key=len)),
i.index(max(i, key=len))]
for i, j in zip(reco_stages, empty_reco_stages[:, 1])
if j is not None])
return max_reco_stages
def get_reco_fit(self, stages):
"""construct a numpy recarray of the fit information (reconstruction
data) of the tracks reconstructed following the reconstruction stages
of interest.
Parameters
----------
stages : list
list of reconstruction stages of interest. for example
[1, 2, 3, 4, 5].
Returns
-------
numpy recarray
a recarray of the fit information (reconstruction data) of
the tracks of interest.
Raises
------
ValueError
ValueError raised when the reconstruction stages of interest
are not found in the file.
"""
keys = ", ".join(self.keys.fit_keys[:-1])
fit_info = self.tracks.fitinf
rec_stages = np.array(
[match for match in self._find_rec_stages(stages)])
if np.all(rec_stages[:, 1] == None):
raise ValueError(
"The stages {} are not found in your file.".format(
str(stages)))
else:
fit_data = np.array([
i[k]
for i, j, k in zip(fit_info, rec_stages[:, 0], rec_stages[:,
1])
if k is not None
])
rec_array = np.core.records.fromarrays(fit_data.transpose(),
names=keys)
return rec_array
def _find_rec_stages(self, stages):
"""find the index of reconstruction stages of interest in a
list of multiple reconstruction stages.
Parameters
----------
stages : list
list of reconstruction stages of interest. for example
[1, 2, 3, 4, 5].
Yieldsma
------
generator
the track id and the index of the reconstruction stages of
interest if found. If the reconstruction stages of interest
are not found, None is returned as the stages index.
"""
for trk_index, rec_stages in enumerate(self.tracks.rec_stages):
try:
stages_index = rec_stages.index(stages)
except ValueError:
stages_index = None
yield trk_index, stages_index
continue
yield trk_index, stages_index
def _find_empty(self, array):
"""finds empty lists/arrays in an awkward array
Parameters
----------
array : awkward array
Awkward array of data of interest. For example:
km3io.OfflineReader(my_file).tracks.fitinf .
Yields
------
generator
the empty list id and the index of the empty list. When
data structure (list) is simply empty, None is written in the
corresponding index. However, when data structure (list) is not
empty and does not contain an empty list, then False is written in the
corresponding index.
"""
for i, rs in enumerate(array):
try:
if len(rs) == 0:
j = None
if len(rs) != 0:
j = rs.index([])
except ValueError:
j = False # rs not empty but [] not found
yield i, j
continue
yield i, j
class OfflineEvents:
"""wrapper for offline events"""
......@@ -549,7 +770,7 @@ class OfflineHit:
class OfflineTracks:
"""wrapper for offline tracks"""
def __init__(self, keys, values, fit_keys=None):
def __init__(self, keys, values, fitparameters=None):
"""wrapper for offline tracks
Parameters
......@@ -558,27 +779,24 @@ class OfflineTracks:
list of cropped tracks keys.
values : list of arrays
list of arrays containting tracks data.
fit_keys : None, optional
list of tracks fit information (not yet outsourced in offline
fitparameters : None, optional
dictionary of tracks fit information (not yet outsourced in offline
files).
"""
self._keys = keys
self._values = values
if fit_keys is not None:
self._fit_keys = fit_keys
if fitparameters is not None:
self._fitparameters = fitparameters
for k, v in zip(self._keys, self._values):
setattr(self, k, v)
def __getitem__(self, item):
if isinstance(item, int):
return OfflineTrack(self._keys, [v[item] for v in self._values],
fit_keys=self._fit_keys)
fitparameters=self._fitparameters)
else:
return OfflineTracks(
self._keys,
[v[item] for v in self._values],
fit_keys=self._fit_keys
)
return OfflineTracks(self._keys, [v[item] for v in self._values],
fitparameters=self._fitparameters)
def __len__(self):
try:
......@@ -596,7 +814,7 @@ class OfflineTracks:
class OfflineTrack:
"""wrapper for an offline track"""
def __init__(self, keys, values, fit_keys=None):
def __init__(self, keys, values, fitparameters=None):
"""wrapper for one offline track.
Parameters
......@@ -605,14 +823,14 @@ class OfflineTrack:
list of cropped tracks keys.
values : list of arrays
list of arrays containting track data.
fit_keys : None, optional
list of tracks fit information (not yet outsourced in offline
fitparameters : None, optional
dictionary of tracks fit information (not yet outsourced in offline
files).
"""
self._keys = keys
self._values = values
if fit_keys is not None:
self._fit_keys = fit_keys
if fitparameters is not None:
self._fitparameters = fitparameters
for k, v in zip(self._keys, self._values):
setattr(self, k, v)
......@@ -621,10 +839,11 @@ class OfflineTrack:
"{:30} {:^2} {:>26}".format(k, ':', str(v))
for k, v in zip(self._keys, self._values) if k not in ['fitinf']
]) + "\n\t" + "\n\t".join([
"{:30} {:^2} {:>26}".format(k, ':', str(v))
for k, v in zip(self._fit_keys, self._values[18]
) # I don't like 18 being explicit here
])
"{:30} {:^2} {:>26}".format(k, ':', str(
getattr(self, 'fitinf')[v]))
for k, v in self._fitparameters.items()
if len(getattr(self, 'fitinf')) > v
]) # I don't like 18 being explicit here
def __getitem__(self, item):
return self._values[item]
......
# -*- coding: utf-8 -*-
"""
KM3NeT Data Definitions v1.0.1
https://git.km3net.de/common/data-format
"""
# reconstruction
JMUONFIT = 0
JMUONBEGIN = 0
JMUONPREFIT = 1
JMUONSIMPLEX = 2
JMUONGANDALF = 3
JMUONENERGY = 4
JMUONSTART = 5
JLINEFIT = 6
JMUONEND = 99
JSHOWERFIT = 100
JSHOWERBEGIN = 100
JSHOWERPREFIT = 101
JSHOWERPOSITIONFIT = 102
JSHOWERCOMPLETEFIT = 103
JSHOWER_BJORKEN_Y = 104
JSHOWEREND = 199
DUSJSHOWERFIT = 200
DUSJBEGIN = 200
DUSJPREFIT = 201
DUSJPOSITIONFIT = 202
JDUSJCOMPLETEFIT = 203
DUSJEND = 299
AASHOWERFIT = 300
AASHOWERBEGIN = 300
AASHOWERCOMPLETEFIT = 301
AASHOWEREND = 399
JUSERBEGIN = 1000
JMUONVETO = 1001
JMUONPATH = 1003
JMCEVT = 1004
JUSEREND = 1099
RECTYPE_UNKNOWN = -1
RECSTAGE_UNKNOWN = -1
# -*- coding: utf-8 -*-
"""
KM3NeT Data Definitions v1.0.1
https://git.km3net.de/common/data-format
"""
# trigger
JTRIGGER3DSHOWER = 1
JTRIGGERMXSHOWER = 2
JTRIGGER3DMUON = 4
JTRIGGERNB = 5
%% Cell type:code id: tags:
``` python
# Add file to current python path
from pathlib import Path
import sys
sys.path.append(str(Path.cwd().parent))
Path.cwd()
# test samples directory - aanet test file
files_path = Path.cwd().parent / 'tests/samples'
aanet_file = files_path / 'aanet_v2.0.0.root'
```
%% Cell type:markdown id: tags:
# Read offline files (aanet)
%% Cell type:code id: tags:
``` python
import km3io as ki
r = ki.OfflineReader(aanet_file)
r
```
%% Output
<km3io.aanet.OfflineReader at 0x7f87bd4ce310>
<km3io.offline.OfflineReader at 0x7fb60c2dc590>
%% Cell type:markdown id: tags:
To explore all the available branches in our offline file:
%% Cell type:code id: tags:
``` python
r.keys
```
%% Output
Events keys are:
id
det_id
mc_id
run_id
mc_run_id
frame_index
trigger_mask
trigger_counter
overlays
hits
trks
w
w2list
w3list
mc_t
mc_hits
mc_trks
comment
index
flags
t.fSec
t.fNanoSec
Hits keys are:
hits.id
hits.dom_id
hits.channel_id
hits.tdc
hits.tot
hits.trig
hits.pmt_id
hits.t
hits.a
hits.pos.x
hits.pos.y
hits.pos.z
hits.dir.x
hits.dir.y
hits.dir.z
hits.pure_t
hits.pure_a
hits.type
hits.origin
hits.pattern_flags
Tracks keys are:
trks.fUniqueID
trks.fBits
trks.id
trks.pos.x
trks.pos.y
trks.pos.z
trks.dir.x
trks.dir.y
trks.dir.z
trks.t
trks.E
trks.len
trks.lik
trks.type
trks.rec_type
trks.rec_stages
trks.status
trks.mother_id
trks.fitinf
trks.hit_ids
trks.error_matrix
trks.comment
Mc hits keys are:
mc_hits.id
mc_hits.dom_id
mc_hits.channel_id
mc_hits.tdc
mc_hits.tot
mc_hits.trig
mc_hits.pmt_id
mc_hits.t
mc_hits.a
mc_hits.pos.x
mc_hits.pos.y
mc_hits.pos.z
mc_hits.dir.x
mc_hits.dir.y
mc_hits.dir.z
mc_hits.pure_t
mc_hits.pure_a
mc_hits.type
mc_hits.origin
mc_hits.pattern_flags
Mc tracks keys are:
mc_trks.fUniqueID
mc_trks.fBits
mc_trks.id
mc_trks.pos.x
mc_trks.pos.y
mc_trks.pos.z
mc_trks.dir.x
mc_trks.dir.y
mc_trks.dir.z
mc_trks.t
mc_trks.E
mc_trks.len
mc_trks.lik
mc_trks.type
mc_trks.rec_type
mc_trks.rec_stages
mc_trks.status
mc_trks.mother_id
mc_trks.fitinf
mc_trks.hit_ids
mc_trks.error_matrix
mc_trks.comment
%% Cell type:markdown id: tags:
# read events data
%% Cell type:code id: tags:
``` python
r.events
```
%% Output
<OfflineEvents: 10 parsed events>
%% Cell type:markdown id: tags:
number of events:
%% Cell type:code id: tags:
``` python
len(r.events)
```
%% Output
10
%% Cell type:code id: tags:
``` python
r.events.id
```
%% Output
<ChunkedArray [1 2 3 ... 8 9 10] at 0x7f87906e86d0>
<ChunkedArray [1 2 3 ... 8 9 10] at 0x7fb5da350b50>
%% Cell type:code id: tags:
``` python
r.events.det_id
```
%% Output
<ChunkedArray [44 44 44 ... 44 44 44] at 0x7f87906e87d0>
<ChunkedArray [44 44 44 ... 44 44 44] at 0x7fb5da3506d0>
%% Cell type:code id: tags:
``` python
r.events.frame_index
```
%% Output
<ChunkedArray [182 183 202 ... 185 185 204] at 0x7f87906e8c10>
<ChunkedArray [182 183 202 ... 185 185 204] at 0x7fb5da350950>
%% Cell type:code id: tags:
``` python
r.events.hits
```
%% Output
<ChunkedArray [176 125 318 ... 84 255 105] at 0x7f87906f9210>
<ChunkedArray [176 125 318 ... 84 255 105] at 0x7fb5da350350>
%% Cell type:markdown id: tags:
lazyarrays can be used with any Numpy universal function. For example:
%% Cell type:code id: tags:
``` python
import numpy as np
np.log(r.events.hits)
```
%% Output
<ChunkedArray [5.170483995038151 4.8283137373023015 5.762051382780177 ... 4.430816798843313 5.541263545158426 4.653960350157523] at 0x7f879052d4d0>
<ChunkedArray [5.170483995038151 4.8283137373023015 5.762051382780177 ... 4.430816798843313 5.541263545158426 4.653960350157523] at 0x7fb60fcdd2d0>
%% Cell type:markdown id: tags:
let's look at event 0:
%% Cell type:code id: tags:
``` python
r.events[0]
```
%% Output
offline event:
id : 1
det_id : 44
mc_id : 0
run_id : 5971
mc_run_id : 0
frame_index : 182
trigger_mask : 22
trigger_counter : 0
overlays : 60
hits : 176
trks : 56
w : []
w2list : []
w3list : []
mc_t : 0.0
mc_hits : 0
mc_trks : 0
comment : b''
index : 0
flags : 0
t_fSec : 1567036818
t_fNanoSec : 200000000
%% Cell type:code id: tags:
``` python
r.events[0].overlays
```
%% Output
60
%% Cell type:code id: tags:
``` python
r.events[0].hits
```
%% Output
176
%% Cell type:markdown id: tags:
# read hits data:
%% Cell type:code id: tags:
``` python
r.hits.dom_id
```
%% Output
<ChunkedArray [[806451572 806451572 806451572 ... 809544061 809544061 809544061] [806451572 806451572 806451572 ... 809524432 809526097 809544061] [806451572 806451572 806451572 ... 809544061 809544061 809544061] ... [806451572 806455814 806465101 ... 809526097 809544058 809544061] [806455814 806455814 806455814 ... 809544061 809544061 809544061] [806455814 806455814 806455814 ... 809544058 809544058 809544061]] at 0x7f87906f9450>
<ChunkedArray [[806451572 806451572 806451572 ... 809544061 809544061 809544061] [806451572 806451572 806451572 ... 809524432 809526097 809544061] [806451572 806451572 806451572 ... 809544061 809544061 809544061] ... [806451572 806455814 806465101 ... 809526097 809544058 809544061] [806455814 806455814 806455814 ... 809544061 809544061 809544061] [806455814 806455814 806455814 ... 809544058 809544058 809544061]] at 0x7fb5da350fd0>
%% Cell type:code id: tags:
``` python
r.hits.tot
```
%% Output
<ChunkedArray [[24 30 22 ... 38 26 23] [29 26 22 ... 26 28 24] [27 19 13 ... 27 24 16] ... [22 22 9 ... 27 32 27] [30 32 17 ... 30 24 29] [27 41 36 ... 29 24 28]] at 0x7f87906f9810>
<ChunkedArray [[24 30 22 ... 38 26 23] [29 26 22 ... 26 28 24] [27 19 13 ... 27 24 16] ... [22 22 9 ... 27 32 27] [30 32 17 ... 30 24 29] [27 41 36 ... 29 24 28]] at 0x7fb5da3a3310>
%% Cell type:code id: tags:
``` python
r[0].hits
```
%% Output
<OfflineHits: 176 parsed elements>
%% Cell type:code id: tags:
``` python
r[0].hits.dom_id
```
%% Output
array([806451572, 806451572, 806451572, 806451572, 806455814, 806455814,
806455814, 806483369, 806483369, 806483369, 806483369, 806483369,
806483369, 806483369, 806483369, 806483369, 806483369, 806487219,
806487226, 806487231, 806487231, 808432835, 808435278, 808435278,
808435278, 808435278, 808435278, 808447180, 808447180, 808447180,
808447180, 808447180, 808447180, 808447180, 808447180, 808447186,
808451904, 808451904, 808472265, 808472265, 808472265, 808472265,
808472265, 808472265, 808472265, 808472265, 808488895, 808488990,
808488990, 808488990, 808488990, 808488990, 808489014, 808489014,
808489117, 808489117, 808489117, 808489117, 808493910, 808946818,
808949744, 808951460, 808951460, 808951460, 808951460, 808951460,
808956908, 808956908, 808959411, 808959411, 808959411, 808961448,
808961448, 808961504, 808961504, 808961655, 808961655, 808961655,
808964815, 808964815, 808964852, 808964908, 808969857, 808969857,
808969857, 808969857, 808969857, 808972593, 808972698, 808972698,
808972698, 808974758, 808974758, 808974758, 808974758, 808974758,
808974758, 808974758, 808974758, 808974758, 808974758, 808974758,
808974773, 808974773, 808974773, 808974773, 808974773, 808974972,
808974972, 808976377, 808976377, 808976377, 808979567, 808979567,
808979567, 808979721, 808979721, 808979721, 808979721, 808979721,
808979721, 808979721, 808979729, 808979729, 808979729, 808981510,
808981510, 808981510, 808981510, 808981672, 808981672, 808981672,
808981672, 808981672, 808981672, 808981672, 808981672, 808981672,
808981672, 808981672, 808981672, 808981672, 808981672, 808981672,
808981672, 808981672, 808981812, 808981812, 808981812, 808981864,
808981864, 808982005, 808982005, 808982005, 808982018, 808982018,
808982018, 808982041, 808982041, 808982077, 808982077, 808982547,
808982547, 808982547, 808997793, 809006037, 809524432, 809526097,
809526097, 809544061, 809544061, 809544061, 809544061, 809544061,
809544061, 809544061], dtype=int32)
%% Cell type:code id: tags:
``` python
r[0].hits[0]
```
%% Output
offline hit:
id : 0
dom_id : 806451572
channel_id : 8
tdc : 0
tot : 24
trig : 1
pmt_id : 0
t : 70104010.0
a : 0.0
pos_x : 0.0
pos_y : 0.0
pos_z : 0.0
dir_x : 0.0
dir_y : 0.0
dir_z : 0.0
pure_t : 0.0
pure_a : 0.0
type : 0
origin : 0
pattern_flags : 0
%% Cell type:code id: tags:
``` python
r[0].hits[0].dom_id
```
%% Output
806451572
%% Cell type:markdown id: tags:
# read tracks data:
%% Cell type:code id: tags:
``` python
r.tracks
```
%% Output
<OfflineTracks: 10 parsed elements>
%% Cell type:code id: tags:
``` python
r.tracks.lik
```
%% Output
<ChunkedArray [[294.6407542676734 294.6407542676734 294.6407542676734 ... 67.81221253265059 67.7756405143316 67.77250505700384] [96.75133289411137 96.75133289411137 96.75133289411137 ... 39.21916536442286 39.184645826013806 38.870325146341884] [560.2775306614813 560.2775306614813 560.2775306614813 ... 118.88577278801066 118.72271313687405 117.80785995187605] ... [71.03251451148226 71.03251451148226 71.03251451148226 ... 16.714140573909347 16.444395245214945 16.34639241716669] [326.440133294878 326.440133294878 326.440133294878 ... 87.79818671079849 87.75488082571873 87.74839444768625] [159.77779654216795 159.77779654216795 159.77779654216795 ... 33.8669134999348 33.821631538334984 33.77240735670646]] at 0x7f8790683d50>
<ChunkedArray [[294.6407542676734 294.6407542676734 294.6407542676734 ... 67.81221253265059 67.7756405143316 67.77250505700384] [96.75133289411137 96.75133289411137 96.75133289411137 ... 39.21916536442286 39.184645826013806 38.870325146341884] [560.2775306614813 560.2775306614813 560.2775306614813 ... 118.88577278801066 118.72271313687405 117.80785995187605] ... [71.03251451148226 71.03251451148226 71.03251451148226 ... 16.714140573909347 16.444395245214945 16.34639241716669] [326.440133294878 326.440133294878 326.440133294878 ... 87.79818671079849 87.75488082571873 87.74839444768625] [159.77779654216795 159.77779654216795 159.77779654216795 ... 33.8669134999348 33.821631538334984 33.77240735670646]] at 0x7fb5da302b50>
%% Cell type:code id: tags:
``` python
r[0].tracks
```
%% Output
<OfflineTracks: 56 parsed elements>
%% Cell type:code id: tags:
``` python
r[0].tracks.lik
```
%% Output
array([294.64075427, 294.64075427, 294.64075427, 291.64653113,
291.27392663, 290.69031512, 289.19290546, 289.08449217,
289.03373947, 288.19030836, 282.92343367, 282.71527118,
282.10762402, 280.20553861, 275.93183966, 273.01809111,
257.46433694, 220.94357656, 194.99426403, 190.47809685,
79.95235686, 78.94389763, 78.90791169, 77.96122466,
77.9579604 , 76.90769883, 75.97546175, 74.91530508,
74.9059469 , 72.94007716, 72.90467038, 72.8629316 ,
72.81280833, 72.80229533, 72.78899435, 71.82404165,
71.80085542, 71.71028058, 70.91130096, 70.89150223,
70.85845637, 70.79081796, 70.76929743, 69.80667603,
69.64058976, 68.93085058, 68.84304037, 68.83154232,
68.79944298, 68.79019375, 68.78581291, 68.72340328,
67.86628937, 67.81221253, 67.77564051, 67.77250506])
%% Cell type:code id: tags:
``` python
r[0].tracks[0]
```
%% Output
offline track:
fUniqueID : 0
fBits : 33554432
id : 1
pos_x : 445.835395997812
pos_y : 615.1089636184813
pos_z : 125.1448339836911
dir_x : 0.0368711082700674
dir_y : -0.48653048395923415
dir_z : -0.872885221293917
t : 70311446.46401498
E : 99.10458562488608
len : 0.0
lik : 294.6407542676734
type : 0
rec_type : 4000
rec_stages : [1, 3, 5, 4]
status : 0
mother_id : -1
hit_ids : []
error_matrix : []
comment : 0
JGANDALF_BETA0_RAD : 0.004957442219414389
JGANDALF_BETA1_RAD : 0.003417848024252858
JGANDALF_CHI2 : -294.6407542676734
JGANDALF_NUMBER_OF_HITS : 142.0
JENERGY_ENERGY : 99.10458562488608
JENERGY_CHI2 : 1.7976931348623157e+308
JGANDALF_LAMBDA : 4.2409761837248484e-12
JGANDALF_NUMBER_OF_ITERATIONS : 10.0
JSTART_NPE_MIP : 24.88469697331908
JSTART_NPE_MIP_TOTAL : 55.88169412579765
JSTART_LENGTH_METRES : 98.89582506402911
JVETO_NPE : 0.0
JVETO_NUMBER_OF_HITS : 0.0
JENERGY_MUON_RANGE_METRES : 344.9767431592819
JENERGY_NOISE_LIKELIHOOD : -333.87773581129136
JENERGY_NDF : 1471.0
JENERGY_NUMBER_OF_HITS : 101.0
%% Cell type:code id: tags:
``` python
r[0].tracks[0].lik
```
%% Output
294.6407542676734
%% Cell type:code id: tags:
``` python
r.tracks.reco.dtype.names
```
%% Output
('JGANDALF_BETA0_RAD',
'JGANDALF_BETA1_RAD',
'JGANDALF_CHI2',
'JGANDALF_NUMBER_OF_HITS',
'JENERGY_ENERGY',
'JENERGY_CHI2',
'JGANDALF_LAMBDA',
'JGANDALF_NUMBER_OF_ITERATIONS',
'JSTART_NPE_MIP',
'JSTART_NPE_MIP_TOTAL',
'JSTART_LENGTH_METRES',
'JVETO_NPE',
'JVETO_NUMBER_OF_HITS',
'JENERGY_MUON_RANGE_METRES',
'JENERGY_NOISE_LIKELIHOOD',
'JENERGY_NDF',
'JENERGY_NUMBER_OF_HITS')
%% Cell type:code id: tags:
``` python
r.tracks.reco["JENERGY_ENERGY"]
```
%% Output
array([99.10458562, 99.10458562, 99.10458562, 37.85515249, 99.10458562,
7.16916787, 99.10458562, 99.10458562, 49.13672986, 20.35137468])
%% Cell type:markdown id: tags:
# read mc hits:
%% Cell type:code id: tags:
``` python
r.mc_hits
```
%% Output
<OfflineHits: 10 parsed elements>
%% Cell type:markdown id: tags:
# read mc tracks:
%% Cell type:code id: tags:
``` python
r.mc_tracks
```
%% Output
<OfflineTracks: 10 parsed elements>
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
......
......@@ -10,4 +10,4 @@ sphinx-autoapi
sphinx-gallery
sphinx-rtd-theme
sphinxcontrib-versioning
yapf
yapf>=0.29.0
......@@ -20,8 +20,8 @@ setup(
url='http://git.km3net.de/km3py/km3io',
description='KM3NeT I/O without ROOT',
long_description=long_description,
author='Zineb Aly, Tamas Gal',
author_email='zaly@km3net.de, tgal@km3net.de',
author='Zineb Aly, Tamas Gal, Johannes Schumann',
author_email='zaly@km3net.de, tgal@km3net.de, johannes.schumann@fau.de',
packages=['km3io'],
include_package_data=True,
platforms='any',
......@@ -33,9 +33,7 @@ setup(
install_requires=requirements,
python_requires='>=3.5',
entry_points={
'console_scripts': [
'KPrintTree=km3io.utils.kprinttree:main'
]
'console_scripts': ['KPrintTree=km3io.utils.kprinttree:main']
},
classifiers=[
'Intended Audience :: Developers',
......@@ -44,4 +42,4 @@ setup(
],
)
__author__ = 'Zineb Aly and Tamas Gal'
__author__ = 'Zineb Aly, Tamas Gal and Johannes Schumann'
......@@ -2,7 +2,7 @@ import os
import re
import unittest
from km3io.daq import DAQReader, get_rate
from km3io.daq import DAQReader, get_rate, has_udp_trailer, get_udp_max_sequence_number, get_channel_flags, get_number_udp_packets
SAMPLES_DIR = os.path.join(os.path.dirname(__file__), "samples")
......@@ -122,7 +122,6 @@ class TestDAQTimeslices(unittest.TestCase):
"daq_v1.0.0.root")).timeslices
def test_data_lengths(self):
assert 3 == len(self.ts._timeslices["default"][0])
assert 3 == len(self.ts._timeslices["L1"][0])
assert 3 == len(self.ts._timeslices["SN"][0])
with self.assertRaises(KeyError):
......@@ -130,12 +129,15 @@ class TestDAQTimeslices(unittest.TestCase):
with self.assertRaises(KeyError):
assert 0 == len(self.ts._timeslices["L0"][0])
def test_streams(self):
self.ts.stream("L1", 0)
self.ts.stream("SN", 0)
def test_reading_frames(self):
assert 8 == len(self.ts.stream("SN", 1).frames[808447186])
def test_str(self):
s = str(self.ts)
assert "default" in s
assert "L1" in s
assert "SN" in s
......@@ -173,8 +175,227 @@ class TestSummaryslices(unittest.TestCase):
def test_rates(self):
assert 3 == len(self.ss.rates)
class TestGetReate(unittest.TestCase):
def test_fifo(self):
s = self.ss.slices[0]
dct_fifo_stat = {
808981510: True,
808981523: False,
808981672: False,
808974773: False
}
for dom_id, fifo_status in dct_fifo_stat.items():
frame = s[s.dom_id == dom_id]
assert any(get_channel_flags(frame.fifo[0])) == fifo_status
def test_has_udp_trailer(self):
s = self.ss.slices[0]
dct_udp_trailer = {
806451572: True,
806455814: True,
806465101: True,
806483369: True,
806487219: True,
806487226: True,
806487231: True,
808432835: True,
808435278: True,
808447180: True,
808447186: True
}
for dom_id, udp_trailer in dct_udp_trailer.items():
frame = s[s.dom_id == dom_id]
assert has_udp_trailer(frame.fifo[0]) == udp_trailer
def test_high_rate_veto(self):
s = self.ss.slices[0]
dct_high_rate_veto = {
808489014: True,
808489117: False,
808493910: True,
808946818: True,
808951460: True,
808956908: True,
808959411: True,
808961448: True,
808961480: True,
808961504: True,
808961655: False,
808964815: False,
808964852: True,
808969848: False,
808969857: True,
808972593: True,
808972598: True,
808972698: False,
808974758: False,
808974773: True,
808974811: True,
808974972: True,
808976377: True,
808979567: False,
808979721: False,
808979729: False,
808981510: True,
808981523: True,
808981672: False,
808981812: True,
808981864: False,
808982018: False
}
for dom_id, high_rate_veto in dct_high_rate_veto.items():
frame = s[s.dom_id == dom_id]
assert any(get_channel_flags(frame.hrv[0])) == high_rate_veto
def test_max_sequence_number(self):
s = self.ss.slices[0]
dct_seq_numbers = {
808974758: 18,
808974773: 26,
808974811: 25,
808974972: 41,
808976377: 35,
808979567: 20,
808979721: 17,
808979729: 25,
808981510: 35,
808981523: 27,
808981672: 17,
808981812: 34,
808981864: 18,
808982018: 21,
808982041: 27,
808982077: 32,
808982547: 20,
808984711: 26,
808996773: 31,
808997793: 21,
809006037: 26,
809007627: 18,
809503416: 28,
809521500: 31,
809524432: 21,
809526097: 23,
809544058: 21,
809544061: 23
}
for dom_id, max_sequence_number in dct_seq_numbers.items():
frame = s[s.dom_id == dom_id]
assert get_udp_max_sequence_number(
frame.dq_status[0]) == max_sequence_number
def test_number_udp_packets(self):
s = self.ss.slices[0]
dct_n_packets = {
808451904: 27,
808451907: 22,
808469129: 20,
808472260: 21,
808472265: 22,
808488895: 20,
808488990: 20,
808489014: 28,
808489117: 22,
808493910: 26,
808946818: 23,
808951460: 37,
808956908: 33,
808959411: 36,
808961448: 28,
808961480: 24,
808961504: 28,
808961655: 20,
808964815: 20,
808964852: 28,
808969848: 21
}
for dom_id, n_udp_packets in dct_n_packets.items():
frame = s[s.dom_id == dom_id]
assert get_number_udp_packets(frame.dq_status[0]) == n_udp_packets
def test_hrv_flags(self):
s = self.ss.slices[0]
dct_hrv_flags = {
809524432: [
False, False, False, False, False, False, False, False, False,
False, False, False, False, False, False, False, False, False,
False, False, False, False, False, False, False, False, False,
False, False, False, False
],
809526097: [
False, False, False, False, False, False, False, False, False,
False, False, False, False, False, False, False, False, False,
True, False, False, False, False, False, False, False, True,
False, False, False, False
],
809544058: [
False, False, False, False, False, False, False, False, False,
False, False, False, False, False, False, False, False, False,
False, False, False, False, False, False, False, False, False,
False, False, False, False
],
809544061: [
False, True, False, False, False, True, False, False, False,
False, False, False, False, False, False, True, False, False,
False, False, False, True, False, False, False, False, False,
False, False, False, False
]
}
for dom_id, hrv_flags in dct_hrv_flags.items():
frame = s[s.dom_id == dom_id]
assert any([
a == b
for a, b in zip(get_channel_flags(frame.hrv[0]), hrv_flags)
])
def test_fifo_flags(self):
s = self.ss.slices[0]
dct_fifo_flags = {
808982547: [
False, False, False, False, False, False, False, False, False,
False, False, False, False, False, False, False, False, False,
False, False, False, False, False, False, False, False, False,
False, False, False, False
],
808984711: [
False, False, False, False, False, False, False, False, False,
False, False, False, False, False, False, False, False, False,
False, False, False, False, False, False, False, False, False,
False, False, False, False
],
808996773: [
False, False, False, False, False, False, False, False, False,
False, False, False, False, False, False, False, False, False,
False, False, False, False, False, False, False, False, False,
False, False, False, False
],
808997793: [
False, False, False, False, False, False, False, False, False,
False, False, False, False, False, False, False, False, False,
False, False, False, False, False, False, False, False, False,
False, False, False, False
],
809006037: [
False, False, False, False, False, False, False, False, False,
False, False, False, False, False, False, False, False, False,
False, False, False, False, False, False, False, False, False,
False, False, False, False
],
808981510: [
False, False, False, False, False, False, False, False, False,
False, False, False, False, False, True, True, False, False,
False, True, False, True, True, True, True, True, True, False,
False, True, False
]
}
for dom_id, fifo_flags in dct_fifo_flags.items():
frame = s[s.dom_id == dom_id]
assert any([
a == b for a, b in zip(
get_channel_flags(frame.fifo[0]), fifo_flags)
])
class TestGetRate(unittest.TestCase):
def test_zero(self):
assert 0 == get_rate(0)
......@@ -185,4 +406,4 @@ class TestGetReate(unittest.TestCase):
def test_vectorized_input(self):
self.assertListEqual([2054], list(get_rate([1])))
self.assertListEqual([2054, 2111, 2169], list(get_rate([1,2,3])))
self.assertListEqual([2054, 2111, 2169], list(get_rate([1, 2, 3])))
import unittest
import numpy as np
from pathlib import Path
from km3io.offline import Reader, OfflineEvents, OfflineHits, OfflineTracks
......@@ -45,6 +46,40 @@ class TestOfflineKeys(unittest.TestCase):
# there are 18 fit keys
self.assertEqual(len(self.keys.fit_keys), 18)
def test_trigger(self):
# there are 4 trigger keys in v1.1.2 of km3net-Dataformat
trigger = self.keys.trigger
keys = [
'JTRIGGER3DSHOWER', 'JTRIGGERMXSHOWER', 'JTRIGGER3DMUON',
'JTRIGGERNB'
]
values = [1, 2, 4, 5]
for k, v in zip(keys, values):
self.assertEqual(v, trigger[k])
def test_reconstruction(self):
# there are 34 parameters in v1.1.2 of km3net-Dataformat
reco = self.keys.reconstruction
keys = [
'JPP_RECONSTRUCTION_TYPE', 'JMUONFIT', 'JMUONBEGIN', 'JMUONPREFIT',
'JMUONSIMPLEX', 'JMUONGANDALF', 'JMUONENERGY', 'JMUONSTART'
]
values = [4000, 0, 0, 1, 2, 3, 4, 5]
self.assertEqual(34, len([*reco.keys()]))
for k, v in zip(keys, values):
self.assertEqual(v, reco[k])
def test_fitparameters(self):
# there are 18 parameters in v1.1.2 of km3net-Dataformat
fit = self.keys.fitparameters
values = [i for i in range(18)]
self.assertEqual(18, len([*fit.keys()]))
for k, v in fit.items():
self.assertEqual(values[v], fit[k])
class TestReader(unittest.TestCase):
def setUp(self):
......@@ -115,6 +150,7 @@ class TestReader(unittest.TestCase):
class TestOfflineReader(unittest.TestCase):
def setUp(self):
self.r = OfflineReader(OFFLINE_FILE)
self.nu = OfflineReader(OFFLINE_NUMUCC)
self.Nevents = 10
self.selected_data = OfflineReader(OFFLINE_FILE,
data=self.r._data[0])._data
......@@ -132,6 +168,66 @@ class TestOfflineReader(unittest.TestCase):
# check that there are 10 events
self.assertEqual(Nevents, self.Nevents)
def test_find_empty(self):
fitinf = self.nu.tracks.fitinf
rec_stages = self.nu.tracks.rec_stages
empty_fitinf = np.array(
[match for match in self.nu._find_empty(fitinf)])
empty_stages = np.array(
[match for match in self.nu._find_empty(rec_stages)])
self.assertListEqual(empty_fitinf[:5, 1].tolist(),
[23, 14, 14, 4, None])
self.assertListEqual(empty_stages[:5, 1].tolist(),
[False, False, False, False, None])
def test_find_rec_stages(self):
stages = np.array(
[match for match in self.nu._find_rec_stages([1, 2, 3, 4, 5])])
self.assertListEqual(stages[:5, 1].tolist(), [0, 0, 0, 0, None])
def test_get_reco_fit(self):
JGANDALF_BETA0_RAD = [
0.0020367251782607574, 0.003306725805622178, 0.0057877124222254885,
0.015581698352185896
]
reco_fit = self.nu.get_reco_fit([1, 2, 3, 4, 5])['JGANDALF_BETA0_RAD']
self.assertListEqual(JGANDALF_BETA0_RAD, reco_fit[:4].tolist())
with self.assertRaises(ValueError):
self.nu.get_reco_fit([1000, 4512, 5625])
def test_get_max_reco_stages(self):
rec_stages = self.nu.tracks.rec_stages
max_reco = self.nu._get_max_reco_stages(rec_stages)
self.assertEqual(len(max_reco.tolist()), 9)
self.assertListEqual(max_reco[0].tolist(), [[1, 2, 3, 4, 5], 5, 0])
def test_best_reco(self):
JGANDALF_BETA1_RAD = [
0.0014177681261476852, 0.002094094517471032, 0.003923368624980349,
0.009491461076780453
]
best = self.nu.best_reco
self.assertEqual(best.size, 9)
self.assertEqual(best['JGANDALF_BETA1_RAD'][:4].tolist(),
JGANDALF_BETA1_RAD)
def test_reading_header(self):
# head is the supported format
head = OfflineReader(OFFLINE_NUMUCC).header
self.assertEqual(float(head['DAQ']), 394)
self.assertEqual(float(head['kcut']), 2)
# test the warning for unsupported fheader format
with self.assertWarns(UserWarning):
self.r.header
class TestOfflineEvents(unittest.TestCase):
def setUp(self):
......@@ -314,16 +410,14 @@ class TestOfflineTracks(unittest.TestCase):
track_selection_2 = tracks[1:3]
assert 2 == len(track_selection_2)
for _slice in [
slice(0, 0),
slice(0, 1),
slice(0, 2),
slice(1, 5),
slice(3, -2)
slice(0, 0),
slice(0, 1),
slice(0, 2),
slice(1, 5),
slice(3, -2)
]:
self.assertListEqual(
list(tracks.E[:,0][_slice]),
list(tracks[_slice].E[:,0])
)
self.assertListEqual(list(tracks.E[:, 0][_slice]),
list(tracks[_slice].E[:, 0]))
class TestOfflineTrack(unittest.TestCase):
......@@ -336,6 +430,4 @@ class TestOfflineTrack(unittest.TestCase):
def test_str(self):
self.assertEqual(repr(self.track).split('\n\t')[0], 'offline track:')
self.assertEqual(
repr(self.track).split('\n\t')[28],
'JGANDALF_LAMBDA : 4.2409761837248484e-12')
self.assertTrue("JGANDALF_LAMBDA" in repr(self.track))