__init__.py

"""
# Event source for MAGIC calibrated data files.
# Requires uproot package (https://github.com/scikit-hep/uproot).
"""

import logging

import glob
import re
import os.path
from pathlib import Path

import numpy as np

import scipy
import scipy.interpolate

from astropy import units as u
from astropy.time import Time

from ctapipe.io.eventsource import EventSource
from ctapipe.io.datalevels import DataLevel

from ctapipe.core import Container
from ctapipe.core import Field

from ctapipe.containers import (
    ArrayEventContainer,
    SimulatedShowerContainer,
    PointingContainer,
    TelescopePointingContainer,
    TelescopeTriggerContainer,
    MonitoringCameraContainer,
    PedestalContainer,
)

from ctapipe.instrument import (
    TelescopeDescription,
    SubarrayDescription,
    OpticsDescription,
    CameraDescription,
)

__all__ = ['MAGICEventSource']

LOGGER = logging.getLogger(__name__)


# MAGIC telescope positions in m wrt. to the center of CTA simulations
#MAGIC_TEL_POSITIONS = {
#    1: [-27.24, -146.66, 50.00] * u.m,
#    2: [-96.44, -96.77, 51.00] * u.m
#}

# MAGIC telescope positions in m wrt. to the center of MAGIC simulations, from reflector input card
MAGIC_TEL_POSITIONS = {
    1: [31.80, -28.10, 0.00] * u.m,
    2: [-31.80, 28.10, 0.00] * u.m
}

# MAGIC telescope description
OPTICS = OpticsDescription.from_name('MAGIC')
MAGICCAM = CameraDescription.from_name("MAGICCam")
GEOM = MAGICCAM.geometry
MAGIC_TEL_DESCRIPTION = TelescopeDescription(
    name='MAGIC', tel_type='MAGIC', optics=OPTICS, camera=MAGICCAM)
MAGIC_TEL_DESCRIPTIONS = {1: MAGIC_TEL_DESCRIPTION, 2: MAGIC_TEL_DESCRIPTION}

# trigger patterns:
MC_TRIGGER_PATTERN = 1
PEDESTAL_TRIGGER_PATTERN = 8
DATA_TRIGGER_PATTERN = 128

class L3JumpError(Exception):
    """
    Exception raised when L3 trigger number jumps backward.
    """

    def __init__(self, message):
        self.message = message

class MissingDriveReportError(Exception):
    """
    Exception raised when a subrun does not have drive reports.
    """

    def __init__(self, message):
        self.message = message

class MAGICEventSource(EventSource):
    """
    EventSource for MAGIC calibrated data.

    This class operates with the MAGIC data run-wise. This means that the files
    corresponding to the same data run are loaded and processed together.
    """
    _count = 0

    def __init__(self, input_url=None, config=None, parent=None, **kwargs):
        """
        Constructor

        Parameters
        ----------
        config: traitlets.loader.Config
            Configuration specified by config file or cmdline arguments.
            Used to set traitlet values.
            Set to None if no configuration to pass.
        parent : ctapipe.core.Tool
            Tool executable that is calling this component.
            Passes the correct logger to the component.
            Set to None if no Tool to pass.
        kwargs: dict
            Additional parameters to be passed.
            NOTE: The file mask of the data to read can be passed with
            the 'input_url' parameter.
        """

        super().__init__(input_url=input_url, config=config, parent=parent, **kwargs)

        # Retrieving the list of run numbers corresponding to the data files
        run_info = self.get_run_info_from_name(str(input_url))
        run_numbers = run_info[0]
        is_mc_runs = run_info[1]

        self.run_numbers = run_numbers
        self.is_mc = is_mc_runs

        # Retrieving the data level (so far HARDCODED Sorcerer)
        self.datalevel = DataLevel.DL1_IMAGES

        # # Setting up the current run with the first run present in the data
        # self.current_run = self._set_active_run(run_number=0)
        self.current_run = None

        self._subarray_info = SubarrayDescription(
            'MAGIC', MAGIC_TEL_POSITIONS, MAGIC_TEL_DESCRIPTIONS)

    @staticmethod
    def is_compatible(file_mask):
        """
        This method checks if the specified file mask corresponds
        to MAGIC data files. The result will be True only if all
        the files are of ROOT format and contain an 'Events' tree.

        Parameters
        ----------
        file_mask: str
            A file mask to check

        Returns
        -------
        bool:
            True if the masked files are MAGIC data runs, False otherwise.

        """

        is_magic_root_file = True

        file_list = glob.glob(str(file_mask))

        for file_path in file_list:
            try:
                import uproot

                try:
                    with uproot.open(file_path) as input_data:
                        if 'Events' not in input_data:
                            is_magic_root_file = False
                except ValueError:
                    # uproot raises ValueError if the file is not a ROOT file
                    is_magic_root_file = False

            except ImportError:
                if re.match(r'.+_m\d_.+root', file_path.lower()) is None:
                    is_magic_root_file = False

        return is_magic_root_file

    @staticmethod
    def get_run_info_from_name(file_name):
        """
        This internal method extracts the run number and
        type (data/MC) from the specified file name.

        Parameters
        ----------
        file_name: str
            A file name to process.

        Returns
        -------
        int:
            A run number of the file.
        """

        mask_data = r".*\d+_M\d+_(\d+)\.\d+_Y_.*"
        mask_mc = r".*_M\d_za\d+to\d+_\d_(\d+)_Y_.*"
        mask_mc_alt = r".*_M\d_\d_(\d+)_.*"
        if re.findall(mask_data, file_name):
            parsed_info = re.findall(mask_data, file_name)
            is_mc = False
        elif re.findall(mask_mc, file_name):
            parsed_info = re.findall(mask_mc, file_name)
            is_mc = True
        else:
            parsed_info = re.findall(mask_mc_alt, file_name)
            is_mc = True

        try:
            run_number = int(parsed_info[0])
        except IndexError:
            raise IndexError(
                'Can not identify the run number and type (data/MC) of the file'
                '{:s}'.format(file_name))

        return run_number, is_mc

    def _set_active_run(self, run_number):
        """
        This internal method sets the run that will be used for data loading.

        Parameters
        ----------
        run_number: int
            The run number to use.

        Returns
        -------
        MarsRun:
            The run to use
        """

        run = dict()
        run['number'] = run_number
        run['read_events'] = 0
        run['data'] = MarsRun(run_file_mask=str(self.input_url))

        return run

    @property
    def subarray(self):
        return self._subarray_info

    @property
    def is_simulation(self):
        return self.is_mc

    @property
    def datalevels(self):
        return (self.datalevel, )

    @property
    def obs_ids(self):
        return self.run_numbers

    def _generator(self):
        """
        The default event generator. Return the stereo event
        generator instance.

        Returns
        -------

        """

        return self._stereo_event_generator()

    def _stereo_event_generator(self):
        """
        Stereo event generator. Yields DataContainer instances, filled
        with the read event data.

        Returns
        -------

        """

        counter = 0

        # Data container - is initialized once, and data is replaced within it after each yield
        data = ArrayEventContainer()

        # Telescopes with data:
        tels_in_file = ["m1", "m2"]
        tels_with_data = {1, 2}

        # Loop over the available data runs
        for run_number in self.run_numbers:

            # Removing the previously read data run from memory
            if self.current_run is not None:
                if 'data' in self.current_run:
                    del self.current_run['data']

            # Setting the new active run (class MarsRun object)
            self.current_run = self._set_active_run(run_number)
            
            # Set monitoring data:
            if not self.is_mc:

                monitoring_data = self.current_run['data'].monitoring_data

                for tel_i, tel_id in enumerate(tels_in_file):
                    monitoring_camera = MonitoringCameraContainer()
                    pedestal_info = PedestalContainer()
                    badpixel_info = PixelStatusContainer()

                    time_tmp = Time(monitoring_data['M{:d}'.format(
                        tel_i + 1)]['PedestalMJD'], scale='utc', format='mjd')
                    pedestal_info.sample_time = Time(
                        time_tmp, format='unix', scale='utc', precision=9)
                    # hardcoded number of pedestal events averaged over:
                    pedestal_info.n_events = 500
                    pedestal_info.charge_mean = []
                    pedestal_info.charge_mean.append(
                        monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFundamental']['Mean'])
                    pedestal_info.charge_mean.append(
                        monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFromExtractor']['Mean'])
                    pedestal_info.charge_mean.append(monitoring_data['M{:d}'.format(
                        tel_i + 1)]['PedestalFromExtractorRndm']['Mean'])
                    pedestal_info.charge_std = []
                    pedestal_info.charge_std.append(
                        monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFundamental']['Rms'])
                    pedestal_info.charge_std.append(
                        monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFromExtractor']['Rms'])
                    pedestal_info.charge_std.append(
                        monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFromExtractorRndm']['Rms'])

                    t_range = Time(monitoring_data['M{:d}'.format(
                        tel_i + 1)]['badpixelinfoMJDrange'], scale='utc', format='mjd')

                    badpixel_info.hardware_failing_pixels = monitoring_data['M{:d}'.format(
                        tel_i + 1)]['badpixelinfo']
                    badpixel_info.sample_time_range = t_range

                    monitoring_camera.pedestal = pedestal_info
                    monitoring_camera.pixel_status = badpixel_info

                    data.mon.tels_with_data = {1, 2}
                    data.mon.tel[tel_i + 1] = monitoring_camera
            else:
                assert self.current_run['data'].mcheader_data['M1'] == self.current_run['data'].mcheader_data['M2'], "Simulation configurations are different for M1 and M2 !!!"
                data.mcheader.num_showers = self.current_run['data'].mcheader_data['M1']['sim_nevents']
                data.mcheader.shower_reuse = self.current_run['data'].mcheader_data['M1']['sim_reuse']
                data.mcheader.energy_range_min = (self.current_run['data'].mcheader_data['M1']['sim_emin']).to(u.TeV) # GeV->TeV
                data.mcheader.energy_range_max = (self.current_run['data'].mcheader_data['M1']['sim_emax']).to(u.TeV) # GeV->TeV
                data.mcheader.spectral_index = self.current_run['data'].mcheader_data['M1']['sim_eslope']
                data.mcheader.max_scatter_range = (self.current_run['data'].mcheader_data['M1']['sim_max_impact']).to(u.m) # cm->m
                data.mcheader.max_viewcone_radius = (self.current_run['data'].mcheader_data['M1']['sim_conesemiangle']).to(u.deg)# deg->deg
                if data.mcheader.max_viewcone_radius != 0.:
                    data.mcheader.diffuse = True
                else:
                    data.mcheader.diffuse = False
                    

            # Loop over the events
            for event_i in range(self.current_run['data'].n_stereo_events):
                # Event and run ids
                event_order_number = self.current_run['data'].stereo_ids[event_i][0]
                event_id = self.current_run['data'].event_data['M1']['stereo_event_number'][event_order_number]
                obs_id = self.current_run['number']

                # Reading event data
                event_data = self.current_run['data'].get_stereo_event_data(
                    event_i)

                data.meta = event_data['mars_meta']

                # Event counter
                data.count = counter
                data.index.obs_id = obs_id
                data.index.event_id = event_id

                # Setting up the R0 container
                data.r0.tel.clear()

                # Setting up the R1 container
                data.r1.tel.clear()

                # Setting up the DL0 container
                data.dl0.tel.clear()

                # Setting up the DL1 container
                data.dl1.tel.clear()

                pointing = PointingContainer()
                # Filling the DL1 container with the event data
                for tel_i, tel_id in enumerate(tels_in_file):
                    # Creating the telescope pointing container
                    pointing_tel = TelescopePointingContainer()
                    
                    pointing_tel.azimuth = np.deg2rad(
                        event_data['{:s}_pointing_az'.format(tel_id)]) * u.rad
                    
                    pointing_tel.altitude = np.deg2rad(
                        90 - event_data['{:s}_pointing_zd'.format(tel_id)]) * u.rad
                    
                    # pointing.ra = np.deg2rad(
                    #    event_data['{:s}_pointing_ra'.format(tel_id)]) * u.rad
                    # pointing.dec = np.deg2rad(
                    #    event_data['{:s}_pointing_dec'.format(tel_id)]) * u.rad

                    pointing.tel[tel_i + 1] = pointing_tel

                    # Adding trigger id (MAGIC nomenclature)
                    data.r0.tel[tel_i + 1].trigger_type = self.current_run['data'].event_data['M1']['trigger_pattern'][event_order_number]
                    data.r1.tel[tel_i + 1].trigger_type = self.current_run['data'].event_data['M1']['trigger_pattern'][event_order_number]
                    data.dl0.tel[tel_i + 1].trigger_type = self.current_run['data'].event_data['M1']['trigger_pattern'][event_order_number]

                    # Adding event charge and peak positions per pixel
                    data.dl1.tel[tel_i +
                                 1].image = event_data['{:s}_image'.format(tel_id)]
                    data.dl1.tel[tel_i +
                                 1].peak_time = event_data['{:s}_pulse_time'.format(tel_id)]
                
                pointing.array_azimuth = np.deg2rad(event_data['m1_pointing_az']) * u.rad
                pointing.array_altitude = np.deg2rad(90 - event_data['m1_pointing_zd']) * u.rad
                pointing.array_ra = np.deg2rad(event_data['m1_pointing_ra']) * u.rad
                pointing.array_dec = np.deg2rad(90 - event_data['m1_pointing_dec']) * u.rad
                data.pointing = pointing

                if not self.is_mc:
                    # Adding the event arrival time
                    time_tmp = Time(
                        event_data['mjd'], scale='utc', format='mjd')
                    data.trigger.time = Time(
                        time_tmp, format='unix', scale='utc', precision=9)
                else:
                    data.mc.energy = event_data['true_energy'] * u.GeV
                    data.mc.alt = (np.pi/2 - event_data['true_zd']) * u.rad
                    # check meaning of 7deg transformation (I.Vovk)
                    data.mc.az = -1 * \
                        (event_data['true_az'] - np.deg2rad(180 - 7)) * u.rad
                    data.mc.shower_primary_id = 1 - \
                        event_data['true_shower_primary_id']
                    data.mc.h_first_int = event_data['true_h_first_int'] * u.cm
                    
                    # adding a 7deg rotation between the orientation of corsika (x axis = magnetic north) and MARS (x axis = geographical north) frames
                    # magnetic north is 7 deg westward w.r.t. geographical north
                    rot_corsika = 7 *u.deg
                    data.mc.core_x = (event_data['true_core_x']*np.cos(rot_corsika) - event_data['true_core_y']*np.sin(rot_corsika))* u.cm
                    data.mc.core_y = (event_data['true_core_x']*np.sin(rot_corsika) + event_data['true_core_y']*np.cos(rot_corsika))* u.cm

                # Setting the telescopes with data
                data.r0.tels_with_data = tels_with_data
                data.r1.tels_with_data = tels_with_data
                data.dl0.tels_with_data = tels_with_data
                data.trigger.tels_with_trigger = tels_with_data

                yield data
                counter += 1

        return

    def _mono_event_generator(self, telescope):
        """
        Mono event generator. Yields DataContainer instances, filled
        with the read event data.

        Parameters
        ----------
        telescope: str
            The telescope for which to return events. Can be either "M1" or "M2".

        Returns
        -------

        """

        counter = 0
        telescope = telescope.upper()

        # Data container - is initialized once, and data is replaced within it after each yield
        data = ArrayEventContainer()

        # Telescopes with data:
        tels_in_file = ["M1", "M2"]

        if telescope not in tels_in_file:
            raise ValueError("Specified telescope {:s} is not in the allowed list {}".format(
                telescope, tels_in_file))

        tel_i = tels_in_file.index(telescope)
        tels_with_data = {tel_i + 1, }

        # Removing the previously read data run from memory
        if self.current_run is not None:
            if 'data' in self.current_run:
                del self.current_run['data']

        # Setting the new active run
        self.current_run = self._set_active_run(self.run_numbers)

        # Set monitoring data:
        if not self.is_mc:

            monitoring_data = self.current_run['data'].monitoring_data

            monitoring_camera = MonitoringCameraContainer()
            pedestal_info = PedestalContainer()
            badpixel_info = PixelStatusContainer()

            time_tmp = Time(monitoring_data['M{:d}'.format(
                tel_i + 1)]['PedestalMJD'], scale='utc', format='mjd')
            pedestal_info.sample_time = Time(
                time_tmp, format='unix', scale='utc', precision=9)
            pedestal_info.n_events = 500 # hardcoded number of pedestal events averaged over
            pedestal_info.charge_mean = []
            pedestal_info.charge_mean.append(
                monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFundamental']['Mean'])
            pedestal_info.charge_mean.append(
                monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFromExtractor']['Mean'])
            pedestal_info.charge_mean.append(monitoring_data['M{:d}'.format(
                tel_i + 1)]['PedestalFromExtractorRndm']['Mean'])
            pedestal_info.charge_std = []
            pedestal_info.charge_std.append(
                monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFundamental']['Rms'])
            pedestal_info.charge_std.append(
                monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFromExtractor']['Rms'])
            pedestal_info.charge_std.append(
                monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFromExtractorRndm']['Rms'])

            t_range = Time(monitoring_data['M{:d}'.format(
                tel_i + 1)]['badpixelinfoMJDrange'], scale='utc', format='mjd')

            badpixel_info.hardware_failing_pixels = monitoring_data['M{:d}'.format(
                tel_i + 1)]['badpixelinfo']
            badpixel_info.sample_time_range = t_range

            monitoring_camera.pedestal = pedestal_info
            monitoring_camera.pixel_status = badpixel_info

            data.mon.tel[tel_i + 1] = monitoring_camera

        #fdp (mono version not fully tested)
        else:
            data.mcheader.num_showers = self.current_run['data'].mcheader_data[telescope]['sim_nevents'] # total, including reuse
            data.mcheader.shower_reuse = self.current_run['data'].mcheader_data[telescope]['sim_reuse']
            data.mcheader.energy_range_min = (self.current_run['data'].mcheader_data[telescope]['sim_emin']).to(u.TeV) # GeV->TeV
            data.mcheader.energy_range_max = (self.current_run['data'].mcheader_data[telescope]['sim_emax']).to(u.TeV) # GeV->TeV
            data.mcheader.spectral_index = self.current_run['data'].mcheader_data[telescope]['sim_eslope']
            data.mcheader.max_scatter_range = (self.current_run['data'].mcheader_data[telescope]['sim_max_impact']).to(u.m) # cm->m
            data.mcheader.max_viewcone_radius = (self.current_run['data'].mcheader_data[telescope]['sim_conesemiangle']).to(u.deg) # deg->deg

        if telescope == 'M1':
            n_events = self.current_run['data'].n_mono_events_m1
        else:
            n_events = self.current_run['data'].n_mono_events_m2

        # Loop over the events
        for event_i in range(n_events):
            # Event and run ids
            event_order_number = self.current_run['data'].mono_ids[telescope][event_i]
            event_id = self.current_run['data'].event_data[telescope]['stereo_event_number'][event_order_number]
            obs_id = self.current_run['number']

            # Reading event data
            event_data = self.current_run['data'].get_mono_event_data(
                event_i, telescope=telescope)

            data.meta = event_data['mars_meta']
            data.meta["max_events"] = self.max_events

            data.trigger.event_type = self.current_run['data'].event_data[telescope]['trigger_pattern'][event_order_number]
            data.trigger.tels_with_trigger = tels_with_data
            if self.allowed_tels:
                data.trigger.tels_with_trigger = np.intersect1d(
                    data.trigger.tels_with_trigger,
                    self.subarray.tel_ids,
                    assume_unique=True,)
            if not self.is_mc:
                # Adding the event arrival time
                time_tmp = Time(
                    event_data['mjd'], scale='utc', format='mjd')
                data.trigger.tel[tel_i + 1] = TelescopeTriggerContainer(time=Time(
                    time_tmp, format='unix', scale='utc', precision=9))

            # Event counter
            data.count = counter
            data.index.obs_id = obs_id
            data.index.event_id = event_id

            # Setting up the R0 container
            data.r0.tel.clear()
            data.r1.tel.clear()
            data.dl0.tel.clear()
            data.dl1.tel.clear()
            data.pointing.tel.clear()

            # Creating the telescope pointing container
            pointing = PointingContainer()
            pointing_tel = TelescopePointingContainer(
                azimuth = np.deg2rad(
                event_data['pointing_az']) * u.rad,
                altitude = np.deg2rad(
                90 - event_data['pointing_zd']) * u.rad,)

            pointing.tel[tel_i + 1] = pointing_tel

            pointing.array_azimuth  = np.deg2rad(event_data['pointing_az']) * u.rad
            pointing.array_altitude = np.deg2rad(90 - event_data['pointing_zd']) * u.rad
            pointing.array_ra       = np.deg2rad(event_data['pointing_ra']) * u.rad
            pointing.array_dec      = np.deg2rad(90 - event_data['pointing_dec']) * u.rad

            data.pointing = pointing

            # Adding event charge and peak positions per pixel
            data.dl1.tel[tel_i + 1].image = event_data['image']
            data.dl1.tel[tel_i + 1].peak_time = event_data['pulse_time']

            if self.is_mc:
                rot_corsika = 7 *u.deg
                # check meaning of 7deg transformation (I.Vovk)
                # adding a 7deg rotation between the orientation of corsika (x axis = magnetic north) and MARS (x axis = geographical north) frames
                # magnetic north is 7 deg westward w.r.t. geographical north
                data.simulation.shower = SimulatedShowerContainer(
                    energy = u.Quantity(event_data['true_energy'], u.GeV),
                    alt = Angle((np.pi/2 - event_data['true_zd']), u.rad),
                    az = Angle(-1 * (event_data['true_az'] - np.deg2rad(180 - 7)), u.rad),
                    shower_primary_id = 1 - event_data['true_shower_primary_id'],
                    h_first_int = u.Quantity(event_data['true_h_first_int'], u.cm),
                    core_x = u.Quantity((event_data['true_core_x']*np.cos(rot_corsika) - event_data['true_core_y']*np.sin(rot_corsika)), u.cm),
                    core_y = u.Quantity((event_data['true_core_x']*np.sin(rot_corsika) + event_data['true_core_y']*np.cos(rot_corsika)), u.cm),
                )

            yield data
            counter += 1

        return

    def _pedestal_event_generator(self, telescope):
        """
        Pedestal event generator. Yields DataContainer instances, filled
        with the read event data.

        Parameters
        ----------
        telescope: str
            The telescope for which to return events. Can be either "M1" or "M2".

        Returns
        -------

        """

        counter = 0
        telescope = telescope.upper()

        # Data container - is initialized once, and data is replaced within it after each yield
        data = EventAndMonDataContainer()

        # Telescopes with data:
        tels_in_file = ["M1", "M2"]

        if telescope not in tels_in_file:
            raise ValueError("Specified telescope {:s} is not in the allowed list {}".format(
                telescope, tels_in_file))

        tel_i = tels_in_file.index(telescope)
        tels_with_data = {tel_i + 1, }

        # Loop over the available data runs
        for run_number in self.run_numbers:

            # Removing the previously read data run from memory
            if self.current_run is not None:
                if 'data' in self.current_run:
                    del self.current_run['data']

            # Setting the new active run
            self.current_run = self._set_active_run(run_number)

            monitoring_data = self.current_run['data'].monitoring_data

            monitoring_camera = MonitoringCameraContainer()
            pedestal_info = PedestalContainer()
            badpixel_info = PixelStatusContainer()

            time_tmp = Time(monitoring_data['M{:d}'.format(
                tel_i + 1)]['PedestalMJD'], scale='utc', format='mjd')
            pedestal_info.sample_time = Time(
                time_tmp, format='unix', scale='utc', precision=9)
            pedestal_info.n_events = 500 # hardcoded number of pedestal events averaged over
            pedestal_info.charge_mean = []
            pedestal_info.charge_mean.append(
                monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFundamental']['Mean'])
            pedestal_info.charge_mean.append(
                monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFromExtractor']['Mean'])
            pedestal_info.charge_mean.append(monitoring_data['M{:d}'.format(
                tel_i + 1)]['PedestalFromExtractorRndm']['Mean'])
            pedestal_info.charge_std = []
            pedestal_info.charge_std.append(
                monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFundamental']['Rms'])
            pedestal_info.charge_std.append(
                monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFromExtractor']['Rms'])
            pedestal_info.charge_std.append(
                monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFromExtractorRndm']['Rms'])

            t_range = Time(monitoring_data['M{:d}'.format(
                tel_i + 1)]['badpixelinfoMJDrange'], scale='utc', format='mjd')

            badpixel_info.hardware_failing_pixels = monitoring_data['M{:d}'.format(
                tel_i + 1)]['badpixelinfo']
            badpixel_info.sample_time_range = t_range

            monitoring_camera.pedestal = pedestal_info
            monitoring_camera.pixel_status = badpixel_info

            data.mon.tels_with_data = tels_with_data
            data.mon.tel[tel_i + 1] = monitoring_camera

            if telescope == 'M1':
                n_events = self.current_run['data'].n_pedestal_events_m1
            else:
                n_events = self.current_run['data'].n_pedestal_events_m2

            # Loop over the events
            for event_i in range(n_events):
                # Event and run ids
                event_order_number = self.current_run['data'].pedestal_ids[telescope][event_i]
                event_id = self.current_run['data'].event_data[telescope]['stereo_event_number'][event_order_number]
                obs_id = self.current_run['number']

                # Reading event data
                event_data = self.current_run['data'].get_pedestal_event_data(
                    event_i, telescope=telescope)

                data.meta = event_data['mars_meta']

                # Event counter
                data.count = counter
                data.index.obs_id = obs_id
                data.index.event_id = event_id

                # Setting up the R0 container
                data.r0.tel.clear()
                data.r0.tel[tel_i + 1].trigger_type = self.current_run['data'].event_data[telescope]['trigger_pattern'][event_order_number]

                # Setting up the R1 container
                data.r1.tel.clear()
                data.r1.tel[tel_i + 1].trigger_type = self.current_run['data'].event_data[telescope]['trigger_pattern'][event_order_number]

                # Setting up the DL0 container
                data.dl0.tel.clear()
                data.dl0.tel[tel_i + 1].trigger_type = self.current_run['data'].event_data[telescope]['trigger_pattern'][event_order_number]

                # Creating the telescope pointing container
                pointing = TelescopePointingContainer()
                pointing.azimuth = np.deg2rad(
                    event_data['pointing_az']) * u.rad
                pointing.altitude = np.deg2rad(
                    90 - event_data['pointing_zd']) * u.rad
                #pointing.ra = np.deg2rad(event_data['pointing_ra']) * u.rad
                #pointing.dec = np.deg2rad(event_data['pointing_dec']) * u.rad

                # Adding the pointing container to the event data
                data.pointing.tel[tel_i + 1] = pointing

                # Adding event charge and peak positions per pixel
                data.dl1.tel[tel_i + 1].image = event_data['image']
                data.dl1.tel[tel_i + 1].peak_time = event_data['pulse_time']

                # Adding the event arrival time
                time_tmp = Time(event_data['mjd'], scale='utc', format='mjd')
                data.trigger.time = Time(
                    time_tmp, format='unix', scale='utc', precision=9)

                # Setting the telescopes with data
                data.r0.tels_with_data = tels_with_data
                data.r1.tels_with_data = tels_with_data
                data.dl0.tels_with_data = tels_with_data
                data.trigger.tels_with_trigger = tels_with_data

                yield data
                counter += 1

        return


class MarsRun:
    """
    This class implements reading of the event data from a single MAGIC data run.
    """

    def __init__(self, run_file_mask, filter_list=None):
        """
        Constructor of the class. Defines the run to use and the camera pixel arrangement.

        Parameters
        ----------
        run_file_mask: str
            A path mask for files belonging to the run. Must correspond to a single run
            or an exception will be raised. Must correspond to calibrated ("sorcerer"-level)
            data.
        filter_list: list, optional
            A list of files, to which the run_file_mask should be applied. If None, all the
            files satisfying run_file_mask will be used. Defaults to None.
        """

        self.n_camera_pixels = GEOM.n_pixels

        self.run_file_mask = run_file_mask

        # Preparing the lists of M1/2 data files
        file_list = glob.glob(run_file_mask)

        # Filtering out extra files if necessary
        if filter_list is not None:
            file_list = list(set(file_list) & set(filter_list))

        self.m1_file_list = list(
            filter(lambda name: '_M1_' in name, file_list))
        self.m2_file_list = list(
            filter(lambda name: '_M2_' in name, file_list))
        self.m1_file_list.sort()
        self.m2_file_list.sort()

        # Retrieving the list of run numbers corresponding to the data files
        run_info = list(
            map(MAGICEventSource.get_run_info_from_name, file_list))
        run_numbers = [i[0] for i in run_info]
        is_mc_runs = [i[1] for i in run_info]

        run_numbers = np.unique(run_numbers)
        is_mc_runs = np.unique(is_mc_runs)
        # Checking if run type (data/MC) is consistent:
        if len(is_mc_runs) > 1:
            raise ValueError(
                "Run type is not consistently data or MC: {}".format(is_mc_runs))

        self.is_mc = is_mc_runs[0]

        # Checking if a single run is going to be read
        if len(run_numbers) > 1:
            raise ValueError(
                "Run mask corresponds to more than one run: {}".format(run_numbers))

        # Reading the data
        m1_data = self.load_events(
            self.m1_file_list, self.is_mc, self.n_camera_pixels)
        m2_data = self.load_events(
            self.m2_file_list, self.is_mc, self.n_camera_pixels)

        # Getting the event data
        self.event_data = dict()
        self.event_data['M1'] = m1_data[0]
        self.event_data['M2'] = m2_data[0]

        # Getting the monitoring data
        self.monitoring_data = dict()
        self.monitoring_data['M1'] = m1_data[1]
        self.monitoring_data['M2'] = m2_data[1]

        # Getting the run-wise MC header data
        if self.is_mc:
            self.mcheader_data = dict()
            self.mcheader_data['M1'] = m1_data[2]
            self.mcheader_data['M2'] = m2_data[2]

        # Detecting pedestal events
        self.pedestal_ids = self._find_pedestal_events()
        # Detecting stereo events
        self.stereo_ids = self._find_stereo_events()
        # Detecting mono events
        self.mono_ids = self._find_mono_events()

    @property
    def n_events_m1(self):
        return len(self.event_data['M1']['MJD'])

    @property
    def n_events_m2(self):
        return len(self.event_data['M2']['MJD'])

    @property
    def n_stereo_events(self):
        return len(self.stereo_ids)

    @property
    def n_mono_events_m1(self):
        return len(self.mono_ids['M1'])

    @property
    def n_mono_events_m2(self):
        return len(self.mono_ids['M2'])

    @property
    def n_pedestal_events_m1(self):
        return len(self.pedestal_ids['M1'])

    @property
    def n_pedestal_events_m2(self):
        return len(self.pedestal_ids['M2'])

    @staticmethod
    def load_events(file_list, is_mc, n_camera_pixels):
        """
        This method loads events and monitoring data from the pre-defiled file and returns them as a dictionary.

        Parameters
        ----------
        file_name: str
            Name of the MAGIC calibrated file to use.
        is_mc: boolean
            Specify whether Monte Carlo (True) or data (False) events are read
        n_camera_pixels: int
            Number of MAGIC camera pixels (not hardcoded, but specified solely via ctapipe.instrument.CameraGeometry)

        Returns
        -------
        dict:
            A dictionary with the even properties: charge / arrival time data, trigger, direction etc.
        """

        try:
            import uproot
        except ImportError:
            msg = "The `uproot` python module is required to access the MAGIC data"
            raise ImportError(msg)

        event_data = dict()

        event_data['charge'] = []
        event_data['arrival_time'] = []
        event_data['trigger_pattern'] = np.array([], dtype=np.int32)
        event_data['stereo_event_number'] = np.array([], dtype=np.int32)
        event_data['pointing_zd'] = np.array([])
        event_data['pointing_az'] = np.array([])
        event_data['pointing_ra'] = np.array([])
        event_data['pointing_dec'] = np.array([])
        event_data['MJD'] = np.array([])
        event_data['mars_meta'] = []

        # monitoring information (updated from time to time)
        monitoring_data = dict()

        monitoring_data['badpixelinfo'] = []
        monitoring_data['badpixelinfoMJDrange'] = []
        monitoring_data['PedestalMJD'] = np.array([])
        monitoring_data['PedestalFundamental'] = dict()
        monitoring_data['PedestalFundamental']['Mean'] = []
        monitoring_data['PedestalFundamental']['Rms'] = []
        monitoring_data['PedestalFromExtractor'] = dict()
        monitoring_data['PedestalFromExtractor']['Mean'] = []
        monitoring_data['PedestalFromExtractor']['Rms'] = []
        monitoring_data['PedestalFromExtractorRndm'] = dict()
        monitoring_data['PedestalFromExtractorRndm']['Mean'] = []
        monitoring_data['PedestalFromExtractorRndm']['Rms'] = []

        #MC Header information, dictionary always created, but filled only in case of MC run
        mcheader_data = dict()

        event_data['file_edges'] = [0]

        # if no file in the list (e.g. when reading mono information), then simply
        # return empty dicts/array
        if len(file_list) == 0:
            return event_data, monitoring_data, mcheader_data

        drive_data = dict()
        drive_data['mjd'] = np.array([])
        drive_data['zd']  = np.array([])
        drive_data['az']  = np.array([])
        drive_data['ra']  = np.array([])
        drive_data['dec'] = np.array([])

        degrees_per_hour = 15.0
        seconds_per_day = 86400.0
        seconds_per_hour = 3600.

        evt_common_list = [
            'MCerPhotEvt.fPixels.fPhot',
            'MArrivalTime.fData',
            'MTriggerPattern.fPrescaled',
            'MRawEvtHeader.fStereoEvtNumber',
            'MRawEvtHeader.fDAQEvtNumber',
        ]

        # Separately, because only used with pre-processed MARS data
        # to create MPointingPos container
        pointing_array_list = [
            'MPointingPos.fZd',
            'MPointingPos.fAz',
            'MPointingPos.fRa',
            'MPointingPos.fDec',
            'MPointingPos.fDevZd',
            'MPointingPos.fDevAz',
            'MPointingPos.fDevHa',
            'MPointingPos.fDevDec',
        ]

        # Info only applicable for data:
        time_array_list = [
            'MTime.fMjd',
            'MTime.fTime.fMilliSec',
            'MTime.fNanoSec',
        ]

        drive_array_list = [
            'MReportDrive.fMjd',
            'MReportDrive.fCurrentZd',
            'MReportDrive.fCurrentAz',
            'MReportDrive.fRa',
            'MReportDrive.fDec'
        ]

        pedestal_array_list = [
            'MTimePedestals.fMjd',
            'MTimePedestals.fTime.fMilliSec',
            'MTimePedestals.fNanoSec',
            'MPedPhotFundamental.fArray.fMean',
            'MPedPhotFundamental.fArray.fRms',
            'MPedPhotFromExtractor.fArray.fMean',
            'MPedPhotFromExtractor.fArray.fRms',
            'MPedPhotFromExtractorRndm.fArray.fMean',
            'MPedPhotFromExtractorRndm.fArray.fRms'
        ]

        # Info only applicable for MC:
        mc_list = [
            'MMcEvt.fEnergy',
            'MMcEvt.fTheta',
            'MMcEvt.fPhi',
            'MMcEvt.fPartId',
            'MMcEvt.fZFirstInteraction',
            'MMcEvt.fCoreX',
            'MMcEvt.fCoreY'
        ]
        
        mcheader_list = [
            #'MMcRunHeader.fNumSimulatedShowers',
            'MMcRunHeader.fNumEvents',
            'MMcCorsikaRunHeader.fELowLim', #GeV
            'MMcCorsikaRunHeader.fEUppLim', #GeV
            'MMcCorsikaRunHeader.fSlopeSpec',
            'MMcRunHeader.fImpactMax', #cm
            #'MMcCorsikaRunHeader.fViewconeAngles',
            'MMcRunHeader.fRandomPointingConeSemiAngle' # deg
        ]

        # Metadata, currently not strictly required
        metainfo_array_list = [
            'MRawRunHeader.fRunNumber',
            'MRawRunHeader.fRunType',
            'MRawRunHeader.fSubRunIndex',
            'MRawRunHeader.fSourceRA',
            'MRawRunHeader.fSourceDEC',
            'MRawRunHeader.fTelescopeNumber']

        for file_name in file_list:

            input_file = uproot.open(file_name)

            events = input_file['Events'].arrays(evt_common_list, library="np")

            # Reading the info common to MC and real data
            charge = events['MCerPhotEvt.fPixels.fPhot']
            arrival_time = events['MArrivalTime.fData']
            trigger_pattern = events['MTriggerPattern.fPrescaled']
            stereo_event_number = events['MRawEvtHeader.fStereoEvtNumber']

            # Reading run-wise meta information (same for all events in subrun):
            mars_meta = dict()

            mars_meta['is_simulation'] = is_mc

            if is_mc:
                mc_header_info = input_file['RunHeaders'].arrays(mcheader_list, library="np")
                mcheader_data['sim_nevents']=int(mc_header_info['MMcRunHeader.fNumEvents'][0]) #std: 5000
                mcheader_data['sim_reuse']= 1 #this value is not written in the magic root file, but since the sim_events already include shower reuse we artificially set it to 1 (actually every shower reused 5 times for std MAGIC MC)
                mcheader_data['sim_emin']=mc_header_info['MMcCorsikaRunHeader.fELowLim'][0]*u.GeV
                mcheader_data['sim_emax']=mc_header_info['MMcCorsikaRunHeader.fEUppLim'][0]*u.GeV
                mcheader_data['sim_eslope']=mc_header_info['MMcCorsikaRunHeader.fSlopeSpec'][0] #std: -1.6
                mcheader_data['sim_max_impact']=mc_header_info['MMcRunHeader.fImpactMax'][0]*u.cm
                mcheader_data['sim_conesemiangle']=mc_header_info['MMcRunHeader.fRandomPointingConeSemiAngle'][0]*u.deg #std: 2.5 deg, also corsika viewcone is defined by "half of the cone angle".

            # Reading event timing information:
            if not is_mc:
                event_times = input_file['Events'].arrays(time_array_list, library="np")
                # Computing the event arrival time

                event_mjd = event_times['MTime.fMjd']
                event_millisec = event_times['MTime.fTime.fMilliSec']
                event_nanosec = event_times['MTime.fNanoSec']

                event_mjd = event_mjd + \
                    (event_millisec / 1e3 + event_nanosec / 1e9) / seconds_per_day
                event_data['MJD'] = np.concatenate(
                    (event_data['MJD'], event_mjd))

            # try to read RunHeaders tree (soft fail if not present, to pass current tests)
            try:
                meta_info = input_file['RunHeaders'].arrays(
                    metainfo_array_list, library="np")

                mars_meta['origin'] = "MAGIC"
                mars_meta['input_url'] = file_name

                mars_meta['number'] = int(
                    meta_info['MRawRunHeader.fRunNumber'][0])
                mars_meta['number_subrun'] = int(
                    meta_info['MRawRunHeader.fSubRunIndex'][0])
                mars_meta['source_ra'] = meta_info['MRawRunHeader.fSourceRA'][0] / \
                    seconds_per_hour * degrees_per_hour * u.deg
                mars_meta['source_dec'] = meta_info['MRawRunHeader.fSourceDEC'][0] / \
                    seconds_per_hour * u.deg

                is_mc_check = int(meta_info['MRawRunHeader.fRunType'][0])
                if is_mc_check == 0:
                    is_mc_check = False
                elif is_mc_check == 256:
                    is_mc_check = True
                else:
                    msg = "Run type (Data or MC) of MAGIC data file not recognised."
                    LOGGER.error(msg)
                    raise ValueError(msg)
                if is_mc_check != is_mc:
                    msg = "Inconsistent run type (data or MC) between file name and runheader content."
                    LOGGER.error(msg)
                    raise ValueError(msg)

                # Reading the info only contained in real data
                if not is_mc:
                    badpixelinfo = input_file['RunHeaders']['MBadPixelsCam.fArray.fInfo'].array(
                        uproot.interpretation.jagged.AsJagged(uproot.interpretation.numerical.AsDtype(np.dtype('>i4'))), library="np")[0].reshape((4, 1183), order='F')
                    # now we have 4 axes:
                    # 0st axis: empty (?)
                    # 1st axis: Unsuitable pixels
                    # 2nd axis: Uncalibrated pixels (says why pixel is unsuitable)
                    # 3rd axis: Bad hardware pixels (says why pixel is unsuitable)
                    # Each axis cointains a 32bit integer encoding more information about the specific problem, see MARS software, MBADPixelsPix.h
                    # take first axis
                    unsuitable_pix_bitinfo = badpixelinfo[1][:n_camera_pixels]
                    # extract unsuitable bit:
                    unsuitable_pix = np.zeros(n_camera_pixels, dtype=np.bool)
                    for i in range(n_camera_pixels):
                        unsuitable_pix[i] = int('\t{0:08b}'.format(
                            unsuitable_pix_bitinfo[i] & 0xff)[-2])
                    monitoring_data['badpixelinfo'].append(unsuitable_pix)
                    # save time interval of badpixel info:
                    monitoring_data['badpixelinfoMJDrange'].append(
                        [event_mjd[0], event_mjd[-1]])

            except KeyError:
                LOGGER.warning(
                    "RunHeaders tree not present in file. Cannot read meta information - will assume it is a real data run.")
                is_mc = False

            # try to read Pedestals tree (soft fail if not present)
            if not is_mc:
                try:
                    pedestal_info = input_file['Pedestals'].arrays(
                        pedestal_array_list, library="np")