Merge pull request #39 from vmarandon/ctapipe018

Ctapipe019

.github/workflows/ci.yml

@@ -18,8 +18,8 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: [3.8]
-        ctapipe-version: [v0.12.0]
+        python-version: [3.9]
+        ctapipe-version: [v0.19.2]
 
     defaults:
       run:

environment.yml

@@ -2,22 +2,23 @@ channels:
   - conda-forge
   - default
 dependencies:
-  - astropy~=4.2
-  - python=3.8 # nail the python version, so conda does not try upgrading / dowgrading
-  - ctapipe=0.12
-  - eventio
-  - corsikaio
-  - protozfits=2.0
+  - ctapipe~=0.19
+  #- astropy~=4.2 # VIM : let the ctapipe version decide
+  #- python=3.8 # nail the python version, so conda does not try upgrading / dowgrading # VIM : let the ctapipe version decide
+  #- eventio
+  #- corsikaio
+  - protozfits
   - zeromq
   - ipython
   - numba
-  - numpy=1.22 # higher versions >=1.23 don't work due to astropy4 dependency in ctapipe0.12
+  #- numpy=1.22 # higher versions >=1.23 don't work due to astropy4 dependency in ctapipe0.12 # VIM : let the ctapipe version decide
   - numpydoc
   - pytest
-  - pytables>=3.7
+  #- pytables>=3.7 # VIM : let the ctapipe version decide
   - pyyaml
   - zlib
   - pip
   - h5py
-  - pip:
-      - pytest_runner
+  - pytest-runner
+#  - pip:
+#      - pytest_runner

setup.cfg

@@ -17,20 +17,20 @@ classifiers =
     Topic :: Scientific/Engineering :: Astronomy
     Topic :: Scientific/Engineering :: Physics
     Programming Language :: Python :: 3 :: Only
-    Programming Language :: Python :: 3.8
+    Programming Language :: Python :: 3.9
 
 [options]
 packages = find:
 package_dir =
-    = src
-python_requires = ~=3.8
+	     = src
+#python_requires = ~=3.8
 zip_safe = False
 install_requires =
-    astropy~=4.2
-    ctapipe~=0.12
-    numpy~=1.22
-    protozfits~=2.0
-    tables>=3.7
+#    astropy~=4.2 # VIM : ctapipe also use this so let's use the same version
+     ctapipe~=0.19
+#    numpy~=1.22 # numpy is required by ctapipe, so let use the same # VIM : ctapipe also use this so let's use the same version
+     protozfits # ~=2.0  
+#    tables>=3.7 # VIM : ctapipe also use this so let's use the same version
 
 [options.package_data]
 * = resources/*
@@ -49,6 +49,11 @@ all =
     %(tests)s
     %(dev)s
 
+# copy from LST Chain... No idea what it is doing
+[options.entry_points]
+ctapipe_io =
+    NectarCAMEventSource = ctapipe_io_nectarcam:NectarCAMEventSource
+
 [tool:pytest]
 minversion = 3.0
 addopts = -v
@@ -57,4 +62,4 @@ addopts = -v
 test = pytest
 
 [flake8]
-max-line-length = 88
+max-line-length = 88

src/ctapipe_io_nectarcam/__init__.py

@@ -13,6 +13,7 @@
 from astropy.time import Time
 from ctapipe.containers import (
     PixelStatusContainer, EventType, R0CameraContainer, R1CameraContainer,
+    CoordinateFrameType, PointingMode, SchedulingBlockContainer, ObservationBlockContainer, # VIM : line of things in lst equivalent
 )
 from ctapipe.coordinates import CameraFrame
 from ctapipe.core import Provenance
@@ -24,8 +25,11 @@
     CameraReadout,
     CameraGeometry,
     OpticsDescription,
+    SizeType,
+    ReflectorShape
 )
-from ctapipe.io import EventSource
+
+from ctapipe.io import EventSource, DataLevel
 from enum import IntFlag, auto
 from pkg_resources import resource_filename
 
@@ -46,7 +50,6 @@
 
 S_TO_NS = np.uint64(1e9)
 
-
 class TriggerBits(IntFlag):
     '''
     See TIB User manual
@@ -86,10 +89,13 @@ class PixelStatus(IntFlag):
     # https://gitlab.cta-observatory.org/cta-science/simulations/simulation-model/verification/verification-process/mst-structure/-/blob/master/Appendix-MST-Structure.pdf
     # version from 20 Jan 2022
     'MST',
+    size_type=SizeType.MST,
     equivalent_focal_length=u.Quantity(16., u.m),
-    num_mirrors=1,
+    effective_focal_length=u.Quantity(16.44505,u.m), # from https://www.mpi-hd.mpg.de/hfm/CTA/MC/Prod5/Config/PSF/fov_prod4.pdf mentioned in the link above
+    n_mirrors=1,
     mirror_area=u.Quantity(106., u.m ** 2),  # no shadowing, uncertainty is 0.5 m2
-    num_mirror_tiles=86,  # Garczarczyk 2017
+    n_mirror_tiles=86,  # Garczarczyk 2017
+    reflector_shape=ReflectorShape.HYBRID #according to Dan Parsons
 )
 
 
@@ -203,22 +209,74 @@ def __init__(self, **kwargs):
             self.file_list.sort()
             kwargs['input_url'] = self.file_list[0]
             super().__init__(**kwargs)
+        elif 'input_filelist' in kwargs.keys():
+            input_filelist = kwargs['input_filelist']
+            if isinstance(input_filelist,str):
+                self.file_list = [input_filelist]
+            else:
+                self.file_list = list(input_filelist)
+            self.file_list.sort()
+            kwargs['input_url'] = self.file_list[0]
+            del kwargs['input_filelist']
+            super().__init__(**kwargs)
         else:
             super().__init__(**kwargs)
             self.file_list = [self.input_url]
 
         self.multi_file = MultiFiles(self.file_list)
         self.camera_config = self.multi_file.camera_config
         self.log.info("Read {} input files".format(self.multi_file.num_inputs()))
-        self._tel_id = self.camera_config.telescope_id
+        self.run_id = self.camera_config.configuration_id
+        self.tel_id = self.camera_config.telescope_id #VIM : Change the _tel_id to tel_id to be consistent with lst
+        self.run_start = Time(self.camera_config.date, format='unix')
         self.geometry_version = 3
-        self._subarray = self.create_subarray(self.geometry_version, self._tel_id)
+        self._subarray = self.create_subarray(self.geometry_version, self.tel_id)
         self.r0_r1_calibrator = NectarCAMR0Corrections(
-            subarray=self._subarray, parent=self
-        )
-        self.nectarcam_service = self.fill_nectarcam_service_container_from_zfile(self._tel_id,
+             subarray=self._subarray, parent=self
+         )
+        self.nectarcam_service = self.fill_nectarcam_service_container_from_zfile(self.tel_id,
                                                                                   self.camera_config)
 
+        target_info = {}
+
+        # VIM : Put some pointing info to be filled as a reminder that it should be done at some point
+        #self.pointing_source = PointingSource(subarray=self.subarray, parent=self)
+        pointing_mode = PointingMode.UNKNOWN
+        #if self.pointing_information:
+        #    target = self.pointing_source.get_target(tel_id=self.tel_id, time=self.run_start)
+        #    if target is not None:
+        #        target_info["subarray_pointing_lon"] = target["ra"]
+        #        target_info["subarray_pointing_lat"] = target["dec"]
+        #        target_info["subarray_pointing_frame"] = CoordinateFrameType.ICRS
+        #        pointing_mode = PointingMode.TRACK
+
+        self._scheduling_blocks = {
+            self.run_id: SchedulingBlockContainer(
+                sb_id=np.uint64(self.run_id),
+                producer_id=f"MST-{self.tel_id}",
+                pointing_mode=pointing_mode,
+            )
+        }
+
+        self._observation_blocks = {
+            self.run_id: ObservationBlockContainer(
+                obs_id=np.uint64(self.run_id),
+                sb_id=np.uint64(self.run_id),
+                producer_id=f"MST-{self.tel_id}",
+                actual_start_time=self.run_start,
+                **target_info
+            )
+        }
+
+    def get_entries(self):
+        tot_events = 0
+        for v in self.multi_file._events_table.values():
+            tot_events += len(v)
+        return tot_events
+
+    def __len__(self):
+        return self.get_entries()
+
     @property
     def is_simulation(self):
         return False
@@ -249,16 +307,19 @@ def create_subarray(geometry_version, tel_id=0):
         daq_time_per_sample, pulse_shape_time_step, pulse_shapes = read_pulse_shapes()
 
         camera_readout = CameraReadout(
-            'NectarCam',
-            1 / daq_time_per_sample,
-            pulse_shapes,
-            pulse_shape_time_step,
+            name = 'NectarCam',
+            n_pixels=N_PIXELS,
+            n_channels=N_GAINS,
+            n_samples=N_SAMPLES,
+            sampling_rate = (1 / daq_time_per_sample).to(u.GHz),
+            reference_pulse_shape = pulse_shapes,
+            reference_pulse_sample_width = pulse_shape_time_step,
         )
 
         camera = CameraDescription('NectarCam', camera_geom, camera_readout)
 
         mst_tel_descr = TelescopeDescription(
-            name='NectarCam', tel_type='MST', optics=OPTICS, camera=camera
+            name='NectarCam', optics=OPTICS, camera=camera
         )
 
         tel_descriptions = {tel_id: mst_tel_descr}
@@ -274,6 +335,15 @@ def create_subarray(geometry_version, tel_id=0):
 
         return subarray
 
+
+    @property
+    def observation_blocks(self):
+        return self._observation_blocks
+
+    @property
+    def scheduling_blocks(self):
+        return self._scheduling_blocks
+
     @property
     def obs_ids(self):
         # currently no obs id is available from the input files
@@ -288,7 +358,7 @@ def _generator(self):
         array_event.meta['origin'] = 'NectarCAM'
 
         # also add service container to the event section
-        array_event.nectarcam.tel[self._tel_id].svc = self.nectarcam_service
+        array_event.nectarcam.tel[self.tel_id].svc = self.nectarcam_service
 
         # initialize general monitoring container
         self.initialize_mon_container(array_event)
@@ -378,7 +448,7 @@ def fill_nectarcam_service_container_from_zfile(tel_id, camera_config):
 
     def fill_nectarcam_event_container_from_zfile(self, array_event, event):
 
-        tel_id = self._tel_id
+        tel_id = self.tel_id
         event_container = NectarCAMEventContainer()
         array_event.nectarcam.tel[tel_id].evt = event_container
 
@@ -429,7 +499,7 @@ def fill_nectarcam_event_container_from_zfile(self, array_event, event):
         self.unpack_feb_data(event_container, event)
 
     def fill_trigger_info(self, array_event):
-        tel_id = self._tel_id
+        tel_id = self.tel_id
 
         nectarcam = array_event.nectarcam.tel[tel_id]
         tib_available = nectarcam.evt.extdevices_presence & 1
@@ -619,8 +689,8 @@ def fill_r0r1_container(self, array_event, zfits_event):
         Fill with R0Container
         """
         r0, r1 = self.fill_r0r1_camera_container(zfits_event)
-        array_event.r0.tel[self._tel_id] = r0
-        array_event.r1.tel[self._tel_id] = r1
+        array_event.r0.tel[self.tel_id] = r0
+        array_event.r1.tel[self.tel_id] = r1
 
     def initialize_mon_container(self, array_event):
         """
@@ -629,7 +699,7 @@ def initialize_mon_container(self, array_event):
 
         """
         container = array_event.mon
-        mon_camera_container = container.tel[self._tel_id]
+        mon_camera_container = container.tel[self.tel_id]
 
         # initialize the container
         status_container = PixelStatusContainer()
@@ -648,7 +718,7 @@ def fill_mon_container_from_zfile(self, array_event, event):
 
         """
 
-        status_container = array_event.mon.tel[self._tel_id].pixel_status
+        status_container = array_event.mon.tel[self.tel_id].pixel_status
 
         # reorder the array
         pixel_status = np.zeros(N_PIXELS)

src/ctapipe_io_nectarcam/calibration.py

@@ -1,7 +1,7 @@
 import numpy as np
 import tables
 from ctapipe.calib.camera.gainselection import ThresholdGainSelector
-from ctapipe.containers import MonitoringContainer
+from ctapipe.containers import MonitoringContainer, MonitoringCameraContainer, FlatFieldContainer, WaveformCalibrationContainer
 from ctapipe.core import TelescopeComponent
 from ctapipe.core.traits import (
     Path, FloatTelescopeParameter, Bool, Float
@@ -31,7 +31,7 @@ class NectarCAMR0Corrections(TelescopeComponent):
     calibration_path = Path(
         default_value=resource_filename(
             'ctapipe_io_nectarcam',
-            'resources/calibrationfile_run3255_pedrun3255_gainrun3155.hdf5'
+            'resources/calibrationfile_run3255_pedrun3255_gainrun3155_ctapipe018.hdf5'
         ),
         exists=True, directory_ok=False, allow_none=True,
         help='Path to calibration file',
@@ -151,6 +151,7 @@ def _read_calibration_file(path):
         Read the correction from hdf5 calibration file
         Only calibration and flatfield containers are filled
         """
+
         mon = MonitoringContainer()
 
         with tables.open_file(path) as f:
@@ -161,14 +162,13 @@ def _read_calibration_file(path):
 
         for tel_id in tel_ids:
             with HDF5TableReader(path) as h5_table:
-                base = f'/tel_{tel_id}'
-                # read the calibration data
-                table = base + '/calibration'
-                next(h5_table.read(table, mon.tel[tel_id].calibration))
-
-                # read flat-field data
-                table = base + '/flatfield'
-                next(h5_table.read(table, mon.tel[tel_id].flatfield))
+
+                mon.tel[tel_id] = MonitoringCameraContainer(
+                    calibration=next(h5_table.read(f'/tel_{tel_id}/calibration', WaveformCalibrationContainer)),
+                    flatfield=next(h5_table.read(f'/tel_{tel_id}/flatfield', FlatFieldContainer)),
+                    #pedestal=next(h5_table.read(f'/{base}/pedestal', PedestalContainer)),
+                    #pixel_status=next(h5_table.read(f"/{base}/pixel_status", PixelStatusContainer)),
+                )
 
         return mon