Merge pull request #136 from cta-observatory/camorgan_v2

Camorgan v2

eventio/scripts/cut_eventio_file.py

@@ -0,0 +1,79 @@
+'''
+Cleanly cut an eventio file, so that the
+uncompressed size is at most `max_size`
+'''
+from argparse import ArgumentParser
+import re
+import gzip
+from ..base import EventIOFile
+
+
+parser = ArgumentParser(description=__doc__)
+parser.add_argument('inputfile', help='Input eventio file')
+parser.add_argument(
+    'outputfile',
+    help='Output file, if ending with .gz, it is written gzip compressed'
+)
+parser.add_argument(
+    'max_size',
+    help=(
+        'Maximimum, uncompressed output file size.'
+        ' You can use k, M, G suffixes.'
+        ' As many complete top-level objects are written, so that max_size is'
+        ' not exceeded'
+    )
+)
+
+
+def parse_size(size):
+    m = re.match(r'(\d+)([kMG])?', size)
+    if not m:
+        raise ValueError('Format not recognized, valid are #, #k, #M or #G')
+
+    size, unit = m.groups()
+    size = float(size)
+
+    if unit is None:
+        return int(size)
+
+    if unit == 'k':
+        return int(size * 1024)
+
+    if unit == 'M':
+        return int(size * 1024**2)
+
+    if unit == 'G':
+        return int(size * 1024**3)
+
+
+def main():
+    args = parser.parse_args()
+
+    max_size = parse_size(args.max_size)
+    bytes_to_read = 0
+    with EventIOFile(args.inputfile) as f:
+        last_o = None
+        for o in f:
+            last_byte = o.header.address + o.header.length
+            if last_byte < max_size:
+                bytes_to_read = last_byte
+                last_o = o
+            else:
+                break
+
+    print('Writing {:.2f}MB (uncompressed)'.format(bytes_to_read / 1024**2))
+    print('Last object is of type {}'.format(last_o.header.type))
+
+    if args.outputfile.endswith('.gz'):
+        open_file = gzip.open
+    else:
+        open_file = open
+
+    with open_file(args.outputfile, 'wb') as of:
+
+        with EventIOFile(args.inputfile) as f:
+            of.write(f.read(bytes_to_read))
+
+
+if __name__ == '__main__':
+    main()

eventio/simtel/objects.py

@@ -237,13 +237,20 @@ def parse(self):
 class CameraOrganization(TelescopeObject):
     eventio_type = 2003
 
-    from .parsing import read_sector_information
+    from .parsing import read_sector_information_v1
+    from ..var_int import read_sector_information_v2
 
     def parse(self):
-        assert_exact_version(self, supported_version=1)
+        assert_version_in(self, supported_versions={1, 2})
+
+        if self.header.version == 1:
+            return self.parse_v1()
+        elif self.header.version == 2:
+            return self.parse_v2()
+
+    def parse_v1(self):
         self.seek(0)
         byte_stream = BytesIO(self.read())
-
         num_pixels = read_int(byte_stream)
         num_drawers = read_int(byte_stream)
         num_gains = read_int(byte_stream)
@@ -262,7 +269,7 @@ def parse(self):
 
         data = read_remaining_with_check(byte_stream, self.header.length)
         pos = 0
-        sectors, bytes_read = CameraOrganization.read_sector_information(
+        sectors, bytes_read = CameraOrganization.read_sector_information_v1(
             data, num_pixels
         )
         pos += bytes_read
@@ -291,6 +298,67 @@ def parse(self):
             'sector_pixthresh': sector_data['pix_thresh'],
         }
 
+    def parse_v2(self):
+        self.seek(0)
+        byte_stream = BytesIO(self.read())
+        num_pixels = read_int(byte_stream)
+        num_drawers = read_int(byte_stream)
+        num_gains = read_int(byte_stream)
+        num_sectors = read_int(byte_stream)
+
+        data = read_remaining_with_check(byte_stream, self.header.length)
+        pos = 0
+
+        drawer, length = varint_array(data, num_pixels, offset=pos)
+        pos += length
+
+        card, length = varint_array(
+            data, num_pixels * num_gains, offset=pos
+        )
+        card.shape = (num_pixels, num_gains)
+        pos += length
+
+        chip, length = varint_array(
+            data, num_pixels * num_gains, offset=pos
+        )
+        chip.shape = (num_pixels, num_gains)
+        pos += length
+
+        channel, length = varint_array(
+            data, num_pixels * num_gains, offset=pos
+        )
+        channel.shape = (num_pixels, num_gains)
+        pos += length
+
+        sectors, bytes_read = CameraOrganization.read_sector_information_v2(
+            data, num_pixels, offset=pos,
+        )
+        pos += bytes_read
+
+        sector_data = np.frombuffer(
+            data,
+            dtype=[
+                ('type', 'uint8'),
+                ('thresh', 'float32'),
+                ('pix_thresh', 'float32')
+            ],
+            count=num_sectors,
+            offset=pos,
+        )
+
+        return {
+            'telescope_id': self.telescope_id,
+            'num_drawers': num_drawers,
+            'drawer': drawer,
+            'card': card,
+            'chip': chip,
+            'channel': channel,
+            'sectors': sectors,
+            'sector_type': sector_data['type'],
+            'sector_threshold': sector_data['thresh'],
+            'sector_pixthresh': sector_data['pix_thresh'],
+        }
+
 
 class PixelSettings(TelescopeObject):
     eventio_type = 2004

eventio/simtel/parsing.pyx

@@ -9,7 +9,7 @@ ctypedef np.int16_t INT16_t
 
 
 @cython.wraparound(False)  # disable negative indexing
-cpdef read_sector_information(
+cpdef read_sector_information_v1(
     const unsigned char[:] data,
     unsigned long n_pixels,
     unsigned long offset = 0,
@@ -21,16 +21,16 @@ cpdef read_sector_information(
     cdef list sectors = []
     cdef array.array sector
 
-    cdef unsigned long pos = 0
+    cdef unsigned long pos = offset
     for i in range(n_pixels):
 
-        short_ptr = <INT16_t*> &data[pos + offset]
+        short_ptr = <INT16_t*> &data[pos]
         n = short_ptr[0]
         pos += 2
 
         sector = array.array('h')
         for j in range(n):
-            short_ptr = <INT16_t*> &data[pos + offset]
+            short_ptr = <INT16_t*> &data[pos]
             sector.append(short_ptr[0])
             pos += 2
 
@@ -46,7 +46,7 @@ cpdef read_sector_information(
         # I will check for it in the tests.
         sectors.append(sector)
 
-    return sectors, pos
+    return sectors, pos - offset
 
 
 

eventio/var_int.pyx

@@ -171,12 +171,12 @@ cpdef varint_array(
     cdef unsigned long bytes_read
     cdef np.ndarray[INT64_t, ndim=1] output = np.empty(n_elements, dtype=INT64)
 
-    cdef UINT64_t val
+    cdef INT64_t val
     cdef unsigned int length
 
-    cdef unsigned long one = 1;
+    cdef int one = 1;
     cdef unsigned long i
-    cdef unsigned long pos = 0
+    cdef unsigned long pos = offset
 
     for i in range(n_elements):
         length = get_length_of_varint(data[pos])
@@ -188,7 +188,7 @@ cpdef varint_array(
         else:
             output[i] = val >> one
 
-    return output, pos
+    return output, pos - offset
 
 
 @cython.wraparound(False)  # disable negative indexing
@@ -497,3 +497,30 @@ cpdef simtel_pixel_timing_parse_list_type_1(
         'pulse_sum_glob': pulse_sum_glob,
         'pulse_sum_loc': pulse_sum_loc,
     }, pos
+
+
+
+@cython.wraparound(False)  # disable negative indexing
+cpdef read_sector_information_v2(
+    const unsigned char[:] data,
+    unsigned long n_pixels,
+    unsigned long offset = 0,
+):
+    cdef unsigned long i
+    cdef unsigned long length
+    cdef int n
+
+    cdef list sectors = []
+    cdef np.ndarray[INT64_t, ndim=1] sector
+
+    cdef unsigned long pos = offset
+    for i in range(n_pixels):
+
+        n, length = varint(data, offset=pos)
+        pos += length
+
+        sector, length = varint_array(data, n, offset=pos)
+        pos += length
+        sectors.append(sector)
+
+    return sectors, pos - offset

setup.py

@@ -50,7 +50,7 @@ def finalize_options(self):
 
 setup(
     name='eventio',
-    version='0.14.0',
+    version='0.15.0',
     description='Python read-only implementation of the EventIO file format',
     long_description=long_description,
     url='https://github.com/fact-project/pyeventio',
@@ -76,6 +76,7 @@ def finalize_options(self):
             'eventio_print_structure = eventio.scripts.print_structure:main',
             'eventio_plot_histograms = eventio.scripts.plot_hists:main',
             'eventio_print_object_information = eventio.scripts.print_object_information:main',
+            'eventio_cut_file = eventio.scripts.cut_eventio_file:main',
         ]
     },
     setup_requires=['pytest-runner', 'numpy'],

tests/simtel/test_simtel_objects.py

@@ -8,9 +8,11 @@
 )
 
 prod2_file = 'tests/resources/gamma_test.simtel.gz'
+camorgan_v2_file = 'tests/resources/test_camorganv2.simtel.gz'
 prod4b_sst1m_file = 'tests/resources/gamma_20deg_0deg_run102___cta-prod4-sst-1m_desert-2150m-Paranal-sst-1m.simtel.gz'
 prod4b_astri_file = 'tests/resources/gamma_20deg_0deg_run103___cta-prod4-sst-astri_desert-2150m-Paranal-sst-astri.simtel.gz'
 
+
 test_files = [
     EventIOFile(path) for path in
     [prod2_file, prod4b_astri_file, prod4b_sst1m_file]
@@ -140,6 +142,24 @@ def test_2003_3_objects():
                 # print(pixel_id, sector)
 
 
+def test_2003_v2():
+    from eventio.simtel.objects import CameraOrganization
+
+    with EventIOFile(camorgan_v2_file) as f:
+        for i, o in enumerate(yield_toplevel_of_type(f, CameraOrganization)):
+            if o.header.version != 2:
+                continue
+
+            cam_organ = parse_and_assert_consumption(o, limit=1)
+            assert cam_organ['telescope_id'] == i + 1
+
+            for sector in cam_organ['sectors']:
+                # sector must never contain a zero, unless it is in the
+                # very first element
+                assert all(s != 0 for s in sector[1:])
+                # print(pixel_id, sector)
+
+
 def test_2004_3_objects():
     from eventio.simtel.objects import PixelSettings