Merge pull request #57 from cta-observatory/ctapipe_v0.17

Adapt to ctapipe v0.17

Author: aleberti

.github/workflows/ci.yml

@@ -12,33 +12,39 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: [3.7]
-        ctapipe-version: [v0.12.0]
+        python-version: ["3.8", "3.9", "3.10"]
+        ctapipe-version: ["v0.17.0",]
+
+    defaults:
+      run:
+        shell: bash -exl {0}
 
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v3
         with:
           fetch-depth: 0
 
-      - name: Set up Python ${{ matrix.python-version }}
-        uses: actions/setup-python@v2
-        with:
+      - name: Set python version
+        env:
           python-version: ${{ matrix.python-version }}
+        run: |
+          sed -i -e "s/- python=.*/- python=$PYTHON_VERSION/g" environment.yml
 
-      - name: Install dependencies
+      - name: Create and activate env
+        uses: conda-incubator/setup-miniconda@v2
+        with:
+          mamba-version: "*"
+          use-mamba: true
+          activate-environment: ci
+          environment-file: environment.yml
+
+      - name: install
         env:
-          PYTHON_VERSION: ${{ matrix.python-version }}
           CTAPIPE_VERSION: ${{ matrix.ctapipe-version }}
-
         run: |
-          . $CONDA/etc/profile.d/conda.sh
-          conda config --set always_yes yes --set changeps1 no
-          sed -i -e "s/- python=.*/- python=$PYTHON_VERSION/g" environment.yml
-          conda env create -n ci -f environment.yml
-          conda activate ci
+          pip install -e .
           # we install ctapipe using pip to be able to select any commit, e.g. the current master
           pip install pytest-cov "git+https://github.com/cta-observatory/ctapipe@$CTAPIPE_VERSION"
-          pip install -e .
           git describe --tags
 
       - name: Download test data
@@ -53,9 +59,7 @@ jobs:
         run: |
           # github actions starts a new shell for each "step", so we need to
           # activate our env again
-          source $CONDA/etc/profile.d/conda.sh
-          conda activate ci
-          coverage run -m pytest -v
-          coverage xml
+          python eventsource_subclasses.py | grep MAGICEventSource
+          pytest --cov=ctapipe_io_magic --cov-report=xml
 
       - uses: codecov/codecov-action@v1

README.md

@@ -1,32 +1,25 @@
 ## *ctapipe* MAGIC event source
 
-EventSource plugin for *ctapipe*, needed to read the calibrated data of the MAGIC telescope system. It requires the [*ctapipe*](https://github.com/cta-observatory/ctapipe) (v0.12.0) and [*uproot*](https://github.com/scikit-hep/uproot4) (>=4.1) packages to run.
+EventSource plugin for *ctapipe*, needed to read the calibrated data of the MAGIC telescope system. It requires the [*ctapipe*](https://github.com/cta-observatory/ctapipe) (v0.17.0) and [*uproot*](https://github.com/scikit-hep/uproot4) (>=5) packages to run.
 
 #### Installation
 
-Provided that *ctapipe* is already installed, the installation can be done via *pip* (the module is available in PyPI):
+If *ctapipe* is already installed, the installation can be done via *pip* (the module is available in PyPI):
 
 ```bash
 pip install ctapipe_io_magic
 ```
 
 Alternatively, you can always clone the repository and install like in the following:
 
-```bash
-git clone https://github.com/cta-observatory/ctapipe_io_magic.git
-pip install ./ctapipe_io_magic/
-```
-
-This installation via *pip* (provided, *pip* is installed) has the advantage to be nicely controlled for belonging to a given conda environment (and to be uninstalled). Alternatively, do
-
 ```bash
 git clone https://github.com/cta-observatory/ctapipe_io_magic.git
 cd ctapipe_io_magic
-python setup.py install --user
+conda env create -n ctapipe-io_magic -f environment.yml
+conda activate ctapipe-io_magic
+pip install .
 ```
 
-In all cases, using *pip* will check if the version of *ctapipe* and *uproot* is compatible with the requested version of *ctapipe_io_magic*.
-
 #### Usage
 
 ```python
@@ -118,3 +111,4 @@ Some general information about the simulated data, useful for IRF calculation, a
 -   v0.4.5: fixed automatic tests, add possibility to choose between effective and nominal focal length
 -   v0.4.6: add support to read in data taken in mono mode (full for real data, partial for MCs). Fixed bug in recognition of mono/stereo or standard trigger/SumT data (added also for MC)
 -   v0.4.7: add full support to read in real and MC data taken in mono mode, and with SumT. Added treatment of unsuitable pixels for MC data. Added readout of true XMax value from MC data (usually not available, filled with 0 otherwise)
+-   v0.5.0: release compatible with ctapipe 0.17. Also, the equivalent focal length is set to the correct value used in MAGIC simulations (i.e. 16.97 meters)

ctapipe_io_magic/__init__.py

@@ -18,14 +18,17 @@
 
 from ctapipe.io import EventSource, DataLevel
 from ctapipe.core import Provenance
-from ctapipe.core.traits import Bool, CaselessStrEnum
+from ctapipe.core.traits import Bool, UseEnum
 from ctapipe.coordinates import CameraFrame
 
 from ctapipe.containers import (
     EventType,
     ArrayEventContainer,
     SimulationConfigContainer,
     SimulatedEventContainer,
+    SchedulingBlockContainer,
+    ObservationBlockContainer,
+    PointingMode,
 )
 
 from ctapipe.instrument import (
@@ -35,6 +38,9 @@
     CameraDescription,
     CameraGeometry,
     CameraReadout,
+    SizeType,
+    ReflectorShape,
+    FocalLengthKind,
 )
 
 from .mars_datalevels import MARSDataLevel
@@ -112,7 +118,7 @@ class MAGICEventSource(EventSource):
         Data level according to MARS convention
     metadata : dict
         Dictionary containing metadata
-    run_numbers : int
+    run_id : int
         Run number of the file
     simulation_config : SimulationConfigContainer
         Container filled with the information about the simulation
@@ -137,9 +143,9 @@ class MAGICEventSource(EventSource):
         help='Use mono events in MC stereo data (needed for mono analysis).'
     ).tag(config=True)
 
-    focal_length_choice = CaselessStrEnum(
-        ['nominal', 'effective'],
-        default_value='effective',
+    focal_length_choice = UseEnum(
+        FocalLengthKind,
+        default_value=FocalLengthKind.EFFECTIVE,
         help='Which focal length to use when constructing the SubarrayDescription.',
     ).tag(config=True)
 
@@ -199,7 +205,7 @@ def __init__(self, input_url=None, config=None, parent=None, **kwargs):
         self.files_ = [uproot.open(rootf) for rootf in self.file_list]
         run_info = self.parse_run_info()
 
-        self.run_numbers = run_info[0]
+        self.run_id = run_info[0][0]
         self.is_mc = run_info[1][0]
         self.telescope = run_info[2][0]
         self.mars_datalevel = run_info[3][0]
@@ -210,7 +216,7 @@ def __init__(self, input_url=None, config=None, parent=None, **kwargs):
         self.datalevel = DataLevel.DL0
 
         if self.is_simulation:
-            self.simulation_config = self.parse_simulation_header()
+            self._simulation_config = self.parse_simulation_header()
 
         self.is_stereo, self.is_sumt = self.parse_data_info()
 
@@ -234,6 +240,24 @@ def __init__(self, input_url=None, config=None, parent=None, **kwargs):
             # Get the arrival time differences
             self.event_time_diffs = self.get_event_time_difference()
 
+        pointing_mode = PointingMode.TRACK
+
+        self._scheduling_blocks = {
+            self.run_id: SchedulingBlockContainer(
+                sb_id=np.uint64(self.run_id),
+                producer_id=f"MAGIC-{self.telescope}",
+                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"MAGIC-{self.telescope}",
+            )
+        }
+
     def __exit__(self, exc_type, exc_val, exc_tb):
         """
         Releases resources (e.g. open files).
@@ -620,23 +644,51 @@ def prepare_subarray_info(self):
             2: [-34.99, 24.02, 0.00] * u.m
         }
 
-        if self.focal_length_choice == 'effective':
-            # Use the effective focal length that the coma aberration is corrected
-            focal_length = u.Quantity(17*1.0713, u.m)
-        else:
-            focal_length = u.Quantity(17, u.m)
+        equivalent_focal_length = u.Quantity(16.97, u.m)
+        effective_focal_length = u.Quantity(17*1.0713, u.m)
 
         OPTICS = OpticsDescription(
-            'MAGIC',
-            num_mirrors=1,
-            equivalent_focal_length=focal_length,
+            name='MAGIC',
+            size_type=SizeType.LST,
+            n_mirrors=1,
+            n_mirror_tiles=964,
+            reflector_shape=ReflectorShape.PARABOLIC,
+            equivalent_focal_length=equivalent_focal_length,
+            effective_focal_length=effective_focal_length,
             mirror_area=u.Quantity(239.0, u.m**2),
-            num_mirror_tiles=964,
         )
 
+        if self.focal_length_choice is FocalLengthKind.EFFECTIVE:
+            focal_length = effective_focal_length
+        elif self.focal_length_choice is FocalLengthKind.EQUIVALENT:
+            focal_length = equivalent_focal_length
+        else:
+            raise ValueError(
+                f"Invalid focal length choice: {self.focal_length_choice}"
+            )
+
         # camera info from MAGICCam.camgeom.fits.gz file
         camera_geom = load_camera_geometry()
 
+        n_pixels = camera_geom.n_pixels
+
+        n_samples_array_list = ["MRawRunHeader.fNumSamplesHiGain"]
+        n_samples_list = []
+
+        for rootf in self.files_:
+            nsample_info = rootf['RunHeaders'].arrays(n_samples_array_list, library="np")
+            n_samples_file = int(nsample_info['MRawRunHeader.fNumSamplesHiGain'][0])
+            n_samples_list.append(n_samples_file)
+
+        n_samples_list = np.unique(n_samples_list).tolist()
+
+        if len(n_samples_list) > 1:
+            raise ValueError(
+                "Loaded files contain different number of readout samples. \
+                 Please load files with the same readout configuration.")
+
+        n_samples = n_samples_list[0]
+
         pulse_shape_lo_gain = np.array([0., 1., 2., 1., 0.])
         pulse_shape_hi_gain = np.array([1., 2., 3., 2., 1.])
         pulse_shape = np.vstack((pulse_shape_lo_gain, pulse_shape_hi_gain))
@@ -645,18 +697,21 @@ def prepare_subarray_info(self):
             u.GHz
         )
         camera_readout = CameraReadout(
-            camera_name='MAGICCam',
+            name='MAGICCam',
             sampling_rate=sampling_speed,
             reference_pulse_shape=pulse_shape,
-            reference_pulse_sample_width=u.Quantity(0.5, u.ns)
+            reference_pulse_sample_width=u.Quantity(0.5, u.ns),
+            n_channels=1,
+            n_pixels=n_pixels,
+            n_samples=n_samples,
         )
 
         camera = CameraDescription('MAGICCam', camera_geom, camera_readout)
 
-        camera.geometry.frame = CameraFrame(focal_length=OPTICS.equivalent_focal_length)
+        camera.geometry.frame = CameraFrame(focal_length=focal_length)
 
         MAGIC_TEL_DESCRIPTION = TelescopeDescription(
-            name='MAGIC', tel_type='MAGIC', optics=OPTICS, camera=camera
+            name='MAGIC', optics=OPTICS, camera=camera
         )
 
         MAGIC_TEL_DESCRIPTIONS = {1: MAGIC_TEL_DESCRIPTION, 2: MAGIC_TEL_DESCRIPTION}
@@ -731,7 +786,7 @@ def parse_metadata_info(self):
 
         metadata = dict()
         metadata["file_list"] = self.file_list
-        metadata['run_numbers'] = self.run_numbers
+        metadata['run_numbers'] = self.run_id
         metadata['is_simulation'] = self.is_simulation
         metadata['telescope'] = self.telescope
         metadata['subrun_number'] = []
@@ -818,7 +873,7 @@ def parse_simulation_header(self):
 
         simulation_config = dict()
 
-        for run_number, rootf in zip(self.run_numbers, self.files_):
+        for rootf in self.files_:
 
             run_header_tree = rootf['RunHeaders']
             spectral_index = run_header_tree['MMcCorsikaRunHeader.fSlopeSpec'].array(library="np")[0]
@@ -848,13 +903,13 @@ def parse_simulation_header(self):
                 max_wavelength = run_header_tree['MMcRunHeader_1.fCWaveUpper'].array(library="np")[0]
                 min_wavelength = run_header_tree['MMcRunHeader_1.fCWaveLower'].array(library="np")[0]
 
-            simulation_config[run_number] = SimulationConfigContainer(
+            simulation_config[self.run_id] = SimulationConfigContainer(
                 corsika_version=corsika_version,
                 energy_range_min=u.Quantity(e_low, u.GeV).to(u.TeV),
                 energy_range_max=u.Quantity(e_high, u.GeV).to(u.TeV),
                 prod_site_alt=u.Quantity(site_height, u.cm).to(u.m),
                 spectral_index=spectral_index,
-                num_showers=n_showers,
+                n_showers=n_showers,
                 shower_reuse=1,
                 # shower_reuse not written in the magic root file, but since the
                 # sim_events already include shower reuse we artificially set it
@@ -968,14 +1023,26 @@ def subarray(self):
     def is_simulation(self):
         return self.is_mc
 
+    @property
+    def observation_blocks(self):
+        return self._observation_blocks
+
+    @property
+    def scheduling_blocks(self):
+        return self._scheduling_blocks
+
+    @property
+    def simulation_config(self):
+        return self._simulation_config
+
     @property
     def datalevels(self):
         return (self.datalevel, )
 
     @property
     def obs_ids(self):
         # ToCheck: will this be compatible in the future, e.g. with merged MC files
-        return self.run_numbers
+        return list(self.observation_blocks)
 
     def _get_badrmspixel_mask(self, event):
         """

ctapipe_io_magic/resources/MAGICCam.camgeom.fits.gz

@@ -171,7 +171,7 @@ def test_run_info(dataset):
         is_mc = [i[1] for i in run_info][0]
         telescope = [i[2] for i in run_info][0]
         datalevel = [i[3] for i in run_info][0]
-        assert run_numbers == source.run_numbers
+        assert run_numbers == [source.run_id]
         assert is_mc == source.is_simulation
         assert telescope == source.telescope
         assert datalevel == source.mars_datalevel
@@ -239,16 +239,17 @@ def test_geom(dataset):
 def test_focal_length_choice(dataset):
     from astropy import units as u
     from ctapipe_io_magic import MAGICEventSource
+    from ctapipe.instrument import FocalLengthKind
 
-    with MAGICEventSource(input_url=dataset, process_run=False, focal_length_choice='nominal') as source:
-        assert source.subarray.tel[1].optics.equivalent_focal_length == u.Quantity(17, u.m)
-        assert source.subarray.tel[2].optics.equivalent_focal_length == u.Quantity(17, u.m)
-        assert source.subarray.tel[1].camera.geometry.frame.focal_length == u.Quantity(17, u.m)
-        assert source.subarray.tel[2].camera.geometry.frame.focal_length == u.Quantity(17, u.m)
+    with MAGICEventSource(input_url=dataset, process_run=False, focal_length_choice=FocalLengthKind.EQUIVALENT) as source:
+        assert source.subarray.tel[1].optics.equivalent_focal_length == u.Quantity(16.97, u.m)
+        assert source.subarray.tel[2].optics.equivalent_focal_length == u.Quantity(16.97, u.m)
+        assert source.subarray.tel[1].camera.geometry.frame.focal_length == u.Quantity(16.97, u.m)
+        assert source.subarray.tel[2].camera.geometry.frame.focal_length == u.Quantity(16.97, u.m)
 
-    with MAGICEventSource(input_url=dataset, process_run=False, focal_length_choice='effective') as source:
-        assert source.subarray.tel[1].optics.equivalent_focal_length == u.Quantity(17*1.0713, u.m)
-        assert source.subarray.tel[2].optics.equivalent_focal_length == u.Quantity(17*1.0713, u.m)
+    with MAGICEventSource(input_url=dataset, process_run=False, focal_length_choice=FocalLengthKind.EFFECTIVE) as source:
+        assert source.subarray.tel[1].optics.effective_focal_length == u.Quantity(17*1.0713, u.m)
+        assert source.subarray.tel[2].optics.effective_focal_length == u.Quantity(17*1.0713, u.m)
         assert source.subarray.tel[1].camera.geometry.frame.focal_length == u.Quantity(17*1.0713, u.m)
         assert source.subarray.tel[2].camera.geometry.frame.focal_length == u.Quantity(17*1.0713, u.m)