From fb1339998549fbb1837598cb1c0315cf0f70d960 Mon Sep 17 00:00:00 2001
From: Ievgen Vovk <vovk@icrr.u-tokyo.ac.jp>
Date: Wed, 18 Oct 2023 18:40:08 +0900
Subject: [PATCH 01/33] Magic events & MC info reading prototype

---
 src/srcsim/magic/mc/__init__.py |   2 +
 src/srcsim/magic/mc/events.py   | 253 ++++++++++++++++++++++++++++++++
 src/srcsim/magic/mc/info.py     |  85 +++++++++++
 3 files changed, 340 insertions(+)
 create mode 100644 src/srcsim/magic/mc/__init__.py
 create mode 100644 src/srcsim/magic/mc/events.py
 create mode 100644 src/srcsim/magic/mc/info.py

diff --git a/src/srcsim/magic/mc/__init__.py b/src/srcsim/magic/mc/__init__.py
new file mode 100644
index 0000000..f1cd4c4
--- /dev/null
+++ b/src/srcsim/magic/mc/__init__.py
@@ -0,0 +1,2 @@
+from .info import MagicMcInfo
+from .events import MagicStereoEvents
\ No newline at end of file
diff --git a/src/srcsim/magic/mc/events.py b/src/srcsim/magic/mc/events.py
new file mode 100644
index 0000000..453b720
--- /dev/null
+++ b/src/srcsim/magic/mc/events.py
@@ -0,0 +1,253 @@
+import numpy
+import uproot
+import pandas
+import astropy.units as u
+from dataclasses import dataclass
+
+
+@dataclass(frozen=True)
+class EventSample:
+    reco_alt: u.Quantity
+    reco_az: u.Quantity
+    reco_ra: u.Quantity
+    reco_dec: u.Quantity
+    reco_energy: u.Quantity
+    src_x: u.Quantity
+    src_y: u.Quantity
+    ra_tel: u.Quantity
+    dec_tel: u.Quantity
+    az_tel: u.Quantity
+    alt_tel: u.Quantity
+    delta_t: u.Quantity
+    gammaness: u.Quantity
+    mjd: u.Quantity = u.Quantity([], unit='d')
+    mc_alt: u.Quantity = u.Quantity([], unit='deg')
+    mc_az: u.Quantity = u.Quantity([], unit='deg')
+    mc_energy: u.Quantity = u.Quantity([], unit='TeV')
+
+
+@dataclass(frozen=True)
+class MagicStereoEvents(EventSample):
+    file_name: str = ''
+
+    def __str__(self) -> str:
+        summary = \
+f"""{type(self).__name__} instance
+    {'File name':.<20s}: {self.file_name}
+    {'MC':.<20s}: {self.is_mc}
+    {'Events':.<20s}: {self.n_events}
+"""
+        return summary
+
+    def __repr__(self):
+        print(str(self))
+
+        return super().__repr__()
+    
+    @property
+    def n_events(self):
+        return len(self.reco_alt)
+    
+    @property
+    def is_mc(self):
+        return len(self.mc_alt) > 0
+
+    @classmethod
+    def from_file(cls, file_name):
+        north_offset = 7
+
+        event_data = dict()
+
+        data_array_list = [
+            'MTriggerPattern_1.fPrescaled',
+            'MRawEvtHeader_1.fStereoEvtNumber',
+            'MRawEvtHeader_1.fDAQEvtNumber',
+            'MRawEvtHeader_1.fTimeDiff',
+            'MPointingPos_1.fZd',
+            'MPointingPos_1.fAz',
+            'MPointingPos_1.fRa',
+            'MPointingPos_1.fDec',
+            'MEnergyEst.fEnergy',
+            'MHadronness.fHadronness',
+            'MRawEvtHeader_1.fDAQEvtNumber',
+            'MSrcPosCam_1.fX',
+            'MSrcPosCam_1.fY',
+            'MStereoParDisp.fDirectionRA',
+            'MStereoParDisp.fDirectionDec',
+            'MStereoParDisp.fDirectionAz',
+            'MStereoParDisp.fDirectionZd',
+        ]
+
+        time_array_list = ['MTime_1.fMjd', 'MTime_1.fTime.fMilliSec', 'MTime_1.fNanoSec']
+
+        mc_array_list = [
+            'MMcEvt_1.fEnergy',
+            'MMcEvt_1.fTheta',
+            'MMcEvt_1.fPhi',
+        ]
+
+        names_mapping = {
+            'MTriggerPattern_1.fPrescaled': 'trigger_pattern',
+            'MRawEvtHeader_1.fStereoEvtNumber': 'stereo_event_number',
+            'MRawEvtHeader_1.fDAQEvtNumber': 'daq_event_number',
+            'MRawEvtHeader_1.fTimeDiff': 'delta_t',
+            'MPointingPos_1.fZd': 'zd_tel',
+            'MPointingPos_1.fAz': 'az_tel',
+            'MPointingPos_1.fRa': 'ra_tel',
+            'MPointingPos_1.fDec': 'dec_tel',
+            'MSrcPosCam_1.fX': 'src_x',
+            'MSrcPosCam_1.fY': 'src_y',
+            'MStereoParDisp.fDirectionRA': 'reco_ra',
+            'MStereoParDisp.fDirectionDec': 'reco_dec',
+            'MStereoParDisp.fDirectionAz': 'reco_az',
+            'MStereoParDisp.fDirectionZd': 'reco_zd',
+            'MHadronness.fHadronness': 'hadronness',
+            'MEnergyEst.fEnergy': 'reco_energy',
+            'MMcEvt_1.fEnergy': 'mc_energy',
+            'MMcEvt_1.fTheta': 'mc_zd',
+            'MMcEvt_1.fPhi': 'mc_az'
+        }
+        
+        data_units = {
+            'delta_t': u.s,
+            'src_x': u.mm,
+            'src_y': u.mm,
+            'reco_alt': u.deg,
+            'reco_az': u.deg,
+            'reco_zd': u.deg,
+            'reco_ra': u.deg,
+            'reco_dec': u.deg,
+            'reco_energy': u.GeV,
+            'gammaness': u.one,
+            'hadronness': u.one,
+            'mjd': u.d,
+            'alt_tel': u.deg,
+            'az_tel': u.deg,
+            'zd_tel': u.deg,
+            'ra_tel': u.deg,
+            'dec_tel':u.deg,
+            'mc_alt':u.deg,
+            'mc_az':u.deg,
+            'mc_zd':u.deg,
+        }
+
+        with uproot.open(file_name) as input_file:
+            if 'Events' in input_file:
+                data = input_file['Events'].arrays(data_array_list)
+
+                # Mapping the read structure to the alternative names
+                for key in data.fields:
+                    name = names_mapping[key]
+                    event_data[name] = data[key]
+
+                is_mc = 'MMcEvt_1.' in input_file['Events']
+                if is_mc:
+                    data = input_file['Events'].arrays(mc_array_list)
+
+                    # Mapping the read structure to the alternative names
+                    for key in data.fields:
+                        name = names_mapping[key]
+                        event_data[name] = data[key]
+
+                    # Post processing
+                    event_data['mc_zd'] = numpy.degrees(event_data['mc_zd'])
+                    event_data['mc_az'] = numpy.degrees(event_data['mc_az'])
+                    # Transformation from Monte Carlo to usual azimuth
+                    event_data['mc_az'] = -1 * (event_data['mc_az'] - 180 + north_offset)
+                else:
+                    # Reading the event arrival time information
+                    data = input_file['Events'].arrays(time_array_list)
+
+                    # Computing the event arrival time
+                    mjd = data['MTime_1.fMjd']
+                    millisec = data['MTime_1.fTime.fMilliSec']
+                    nanosec = data['MTime_1.fNanoSec']
+
+                    event_data['mjd'] = mjd + (millisec / 1e3 + nanosec / 1e9) / 86400.0
+
+            else:
+                # The file is likely corrupted, so return empty arrays
+                for key in names_mapping:
+                    name = names_mapping[key]
+                    event_data[name] = numpy.zeros(0)
+                    
+        event_data['mc_alt'] = 90 - event_data['mc_zd']
+        event_data['alt_tel'] = 90 - event_data['zd_tel']
+        event_data['reco_alt'] = 90 - event_data['reco_zd']
+        event_data['gammaness'] = 1 - event_data['hadronness']
+        
+        for key in event_data:
+            event_data[key] = event_data[key].to_numpy()
+            
+        for key in event_data:
+            if key in data_units:
+                event_data[key] = event_data[key] * data_units[key]
+
+        is_mc = 'mc_energy' in event_data
+
+        if is_mc:
+            events = MagicStereoEvents(
+                reco_alt = event_data['reco_alt'],
+                reco_az = event_data['reco_az'],
+                reco_ra = event_data['reco_ra'],
+                reco_dec = event_data['reco_dec'],
+                reco_energy = event_data['reco_energy'],
+                src_x = event_data['src_x'],
+                src_y = event_data['src_y'],
+                ra_tel = event_data['ra_tel'],
+                dec_tel = event_data['dec_tel'],
+                az_tel = event_data['az_tel'],
+                alt_tel = event_data['alt_tel'],
+                delta_t = event_data['delta_t'],
+                gammaness = event_data['gammaness'],
+                mc_alt = event_data['mc_alt'],
+                mc_az = event_data['mc_az'],
+                mc_energy = event_data['mc_energy'],
+                file_name = file_name
+            )
+        else:
+            events = MagicStereoEvents(
+                reco_alt = event_data['reco_alt'],
+                reco_az = event_data['reco_az'],
+                reco_ra = event_data['reco_ra'],
+                reco_dec = event_data['reco_dec'],
+                reco_energy = event_data['reco_energy'],
+                src_x = event_data['src_x'],
+                src_y = event_data['src_y'],
+                ra_tel = event_data['ra_tel'],
+                dec_tel = event_data['dec_tel'],
+                az_tel = event_data['az_tel'],
+                alt_tel = event_data['alt_tel'],
+                delta_t = event_data['delta_t'],
+                mjd = event_data['mjd'],
+                gammaness = event_data['gammaness'],
+                file_name = file_name
+            )
+
+        return events
+    
+    def to_df(self):
+        units = dict(
+            # delta_t = u.s,
+            src_x = u.mm,
+            src_y = u.mm,
+            reco_alt = u.deg,
+            reco_az = u.deg,
+            reco_ra = u.deg,
+            reco_dec = u.deg,
+            reco_energy = u.GeV,
+            gammaness = u.one,
+            mjd = u.d,
+            alt_tel = u.deg,
+            az_tel = u.deg,
+            ra_tel = u.deg,
+            dec_tel =u.deg,
+            mc_alt =u.deg,
+            mc_az =u.deg,
+        )
+        data = {
+            key: self.__getattribute__(key).to(units[key]).value
+            for key in units
+            if len(self.__getattribute__(key))
+        }
+        return pandas.DataFrame(data=data)
diff --git a/src/srcsim/magic/mc/info.py b/src/srcsim/magic/mc/info.py
new file mode 100644
index 0000000..82f24b7
--- /dev/null
+++ b/src/srcsim/magic/mc/info.py
@@ -0,0 +1,85 @@
+import uproot
+import pandas
+import astropy.units as u
+from dataclasses import dataclass
+
+
+@dataclass(frozen=True)
+class McInfo:
+    num_showers: u.Quantity
+    energy_range_min: u.Quantity
+    energy_range_max: u.Quantity
+    spectral_index: u.Quantity
+    max_scatter_range: u.Quantity
+    min_scatter_range: u.Quantity
+    max_viewcone_radius: u.Quantity
+    min_viewcone_radius: u.Quantity
+
+
+@dataclass(frozen=True)
+class MagicMcInfo(McInfo):
+    file_name: str = ''
+
+    def __repr__(self):
+        print(
+f"""{type(self).__name__} instance
+    {'File name':.<20s}: {self.file_name}
+    {'Events':.<20s}: {self.num_showers}
+    {'Energy min':.<20s}: {self.energy_range_min}
+    {'Energy max':.<20s}: {self.energy_range_max}
+    {'Spectral index':.<20s}: {self.spectral_index}
+    {'Scatter range min':.<20s}: {self.min_scatter_range}
+    {'Scatter range max':.<20s}: {self.max_scatter_range}
+    {'Viewcone min':.<20s}: {self.min_viewcone_radius}
+    {'Viewcone max':.<20s}: {self.max_viewcone_radius}
+"""
+        )
+
+        return super().__repr__()
+
+    @classmethod
+    def from_file(cls, file_name):
+        meta = dict()
+        with uproot.open(file_name) as input_file:
+            meta['num_showers'] = u.one * int(
+                input_file['RunHeaders']['MMcRunHeader_1./MMcRunHeader_1.fNumSimulatedShowers'].array()[0]
+            )
+            meta['energy_range_min'] = u.GeV * input_file['RunHeaders']['MMcCorsikaRunHeader./MMcCorsikaRunHeader.fELowLim'].array()[0]
+            meta['energy_range_max'] = u.GeV * input_file['RunHeaders']['MMcCorsikaRunHeader./MMcCorsikaRunHeader.fEUppLim'].array()[0]
+            meta['spectral_index'] = u.one * input_file['RunHeaders']['MMcRunHeader_1.fSlopeSpec'].array()[0]
+            meta['max_scatter_range'] = u.cm * input_file['RunHeaders']['MMcRunHeader_1.fImpactMax'].array()[0]
+            meta['min_scatter_range'] = u.cm * 0
+            meta['max_viewcone_radius'] = u.deg * input_file['RunHeaders']['MMcRunHeader_1./MMcRunHeader_1.fRandomPointingConeSemiAngle'].array()[0]
+            meta['min_viewcone_radius'] = u.deg * 0
+
+        info = MagicMcInfo(
+            num_showers = meta['num_showers'],
+            energy_range_min = meta['energy_range_min'],
+            energy_range_max = meta['energy_range_max'],
+            spectral_index = meta['spectral_index'],
+            max_scatter_range = meta['max_scatter_range'],
+            min_scatter_range = meta['min_scatter_range'],
+            max_viewcone_radius = meta['max_viewcone_radius'],
+            min_viewcone_radius = meta['min_viewcone_radius']
+        )
+
+        return info
+    
+    def to_df(self):
+        units = dict(
+            num_showers = u.one,
+            energy_range_min = u.TeV,
+            energy_range_max = u.TeV,
+            spectral_index = u.one,
+            max_scatter_range = u.m,
+            min_scatter_range = u.m,
+            max_viewcone_radius = u.deg,
+            min_viewcone_radius = u.deg,
+        )
+        data = {
+            key: [
+                self.__getattribute__(key).to(units[key]).value
+            ]
+            for key in units
+        }
+        return pandas.DataFrame(data=data)

From 1278f6cc4dc4bd1f2a4a1acd213577cf8fe9f825 Mon Sep 17 00:00:00 2001
From: Ievgen Vovk <vovk@icrr.u-tokyo.ac.jp>
Date: Thu, 19 Oct 2023 11:03:02 +0900
Subject: [PATCH 02/33] Magic "original" MC tree reader.

---
 src/srcsim/magic/mc/__init__.py |  2 +-
 src/srcsim/magic/mc/events.py   | 85 +++++++++++++++++++++++++++++++++
 2 files changed, 86 insertions(+), 1 deletion(-)

diff --git a/src/srcsim/magic/mc/__init__.py b/src/srcsim/magic/mc/__init__.py
index f1cd4c4..88d8473 100644
--- a/src/srcsim/magic/mc/__init__.py
+++ b/src/srcsim/magic/mc/__init__.py
@@ -1,2 +1,2 @@
 from .info import MagicMcInfo
-from .events import MagicStereoEvents
\ No newline at end of file
+from .events import MagicStereoEvents, MagicMcOrigEvents
\ No newline at end of file
diff --git a/src/srcsim/magic/mc/events.py b/src/srcsim/magic/mc/events.py
index 453b720..7797302 100644
--- a/src/srcsim/magic/mc/events.py
+++ b/src/srcsim/magic/mc/events.py
@@ -26,6 +26,17 @@ class EventSample:
     mc_energy: u.Quantity = u.Quantity([], unit='TeV')
 
 
+@dataclass(frozen=True)
+class McOrigEventSample:
+    az_tel: u.Quantity
+    alt_tel: u.Quantity
+    energy: u.Quantity
+    src_x: u.Quantity = u.Quantity([], unit='m')
+    src_y: u.Quantity = u.Quantity([], unit='m')
+    alt: u.Quantity = u.Quantity([], unit='deg')
+    az: u.Quantity = u.Quantity([], unit='deg')
+
+
 @dataclass(frozen=True)
 class MagicStereoEvents(EventSample):
     file_name: str = ''
@@ -251,3 +262,77 @@ def to_df(self):
             if len(self.__getattribute__(key))
         }
         return pandas.DataFrame(data=data)
+
+
+@dataclass(frozen=True)
+class MagicMcOrigEvents(McOrigEventSample):
+    file_name: str = ''
+
+    def __str__(self) -> str:
+        summary = \
+f"""{type(self).__name__} instance
+    {'File name':.<20s}: {self.file_name}
+    {'Events':.<20s}: {self.n_events}
+"""
+        return summary
+
+    def __repr__(self):
+        print(str(self))
+
+        return super().__repr__()
+
+    @property
+    def n_events(self):
+        return len(self.az_tel)
+
+    @classmethod
+    def from_file(cls, file_name):
+        north_offset = 7 * u.deg
+
+        data_array_list = [
+            'MMcEvtBasic_1.fEnergy',
+            'MMcEvtBasic_1.fTelescopePhi',
+            'MMcEvtBasic_1.fTelescopeTheta',
+            'MSrcPosCam_1.fX',
+            'MSrcPosCam_1.fY',
+        ]
+
+        data_units = {
+            'MMcEvtBasic_1.fEnergy': u.GeV,
+            'MSrcPosCam_1.fX': u.mm,
+            'MSrcPosCam_1.fY': u.mm,
+            'MMcEvtBasic_1.fTelescopePhi': u.rad,
+            'MMcEvtBasic_1.fTelescopeTheta': u.rad,
+        }
+
+        with uproot.open(file_name) as input_file:
+            if 'OriginalMC' in input_file:
+                data = input_file['OriginalMC'].arrays(data_array_list)
+            else:
+                # The file is likely corrupted, so return empty arrays
+                data = {
+                    key: numpy.zeros(0)
+                    for key in data_array_list
+                }
+
+        event_data = {
+            key: data[key].to_numpy() * data_units[key]
+            for key in data.fields
+        }
+
+        event_data['MMcEvtBasic_1.fTelescopeAlt'] = numpy.pi/2 * u.rad - event_data['MMcEvtBasic_1.fTelescopeTheta']
+        # Transformation from Monte Carlo to usual azimuth
+        data['MMcEvtBasic_1.fTelescopePhi'] = -1 * (
+            data['MMcEvtBasic_1.fTelescopePhi'] - (180 * u.deg + north_offset).to(data_units['MMcEvtBasic_1.fTelescopePhi'])
+        )
+
+        events = MagicMcOrigEvents(
+            az_tel = event_data['MMcEvtBasic_1.fTelescopePhi'],
+            alt_tel = event_data['MMcEvtBasic_1.fTelescopeAlt'],
+            energy = event_data['MMcEvtBasic_1.fEnergy'],
+            src_x = event_data['MSrcPosCam_1.fX'],
+            src_y = event_data['MSrcPosCam_1.fY'],
+            file_name = file_name
+        )
+
+        return events
\ No newline at end of file

From 23cfdb6472f50103cdef1f23aaa155d61d170dbc Mon Sep 17 00:00:00 2001
From: Ievgen Vovk <vovk@icrr.u-tokyo.ac.jp>
Date: Thu, 19 Oct 2023 12:31:02 +0900
Subject: [PATCH 03/33] High-level MAGIC MC file class.

---
 src/srcsim/magic/mc/__init__.py |   3 +-
 src/srcsim/magic/mc/file.py     | 111 ++++++++++++++++++++++++++++++++
 2 files changed, 113 insertions(+), 1 deletion(-)
 create mode 100644 src/srcsim/magic/mc/file.py

diff --git a/src/srcsim/magic/mc/__init__.py b/src/srcsim/magic/mc/__init__.py
index 88d8473..147d797 100644
--- a/src/srcsim/magic/mc/__init__.py
+++ b/src/srcsim/magic/mc/__init__.py
@@ -1,2 +1,3 @@
 from .info import MagicMcInfo
-from .events import MagicStereoEvents, MagicMcOrigEvents
\ No newline at end of file
+from .events import MagicStereoEvents, MagicMcOrigEvents
+from .file import MagicMcFile
\ No newline at end of file
diff --git a/src/srcsim/magic/mc/file.py b/src/srcsim/magic/mc/file.py
new file mode 100644
index 0000000..48827a2
--- /dev/null
+++ b/src/srcsim/magic/mc/file.py
@@ -0,0 +1,111 @@
+import numpy
+import healpy
+from dataclasses import dataclass
+
+from .info import MagicMcInfo
+from .events import MagicStereoEvents, MagicMcOrigEvents
+
+
+@dataclass(frozen=True)
+class MagicMcFile:
+    file_name: str
+    meta: MagicMcInfo
+    events_triggered: MagicStereoEvents
+    events_simulated: MagicMcOrigEvents
+
+    @classmethod
+    def from_file(cls, file_name):
+        self = MagicMcFile(
+            file_name = file_name,
+            meta = MagicMcInfo.from_file(file_name),
+            events_triggered = MagicStereoEvents.from_file(file_name),
+            events_simulated = MagicMcOrigEvents.from_file(file_name)
+        )
+        return self
+
+    def write(self, file_name):
+        self.meta.to_df().to_hdf(
+            file_name,
+            key='/simulation/run_config',
+            mode='w'
+        )
+        self.events_triggered.to_df().to_hdf(
+            file_name,
+            key='/dl2/event/telescope/parameters/MAGIC_MAGICCam',
+            mode='a'
+        )
+
+    def healpy_split(self, nside):
+        npix = healpy.nside2npix(nside)
+
+        pix_ids_triggered = healpy.ang2pix(
+            nside,
+            numpy.pi/2 - self.events_triggered.alt_tel.to('rad').value,
+            self.events_triggered.az_tel.to('rad').value,
+        )
+
+        pix_ids_simulated = healpy.ang2pix(
+            nside,
+            numpy.pi/2 - self.events_simulated.alt_tel.to('rad').value,
+            self.events_simulated.az_tel.to('rad').value,
+        )
+
+        files = []
+
+        for pix_id in range(npix):
+            n_simulated = numpy.sum(pix_ids_simulated == pix_id)
+            meta = MagicMcInfo(
+                num_showers = n_simulated,
+                energy_range_min = self.meta.energy_range_min,
+                energy_range_max = self.meta.energy_range_max,
+                spectral_index = self.meta.spectral_index,
+                max_scatter_range = self.meta.max_scatter_range,
+                min_scatter_range = self.meta.min_scatter_range,
+                max_viewcone_radius = self.meta.max_viewcone_radius,
+                min_viewcone_radius = self.meta.min_viewcone_radius
+            )
+
+            data = self.events_simulated
+            selection = pix_ids_simulated == pix_id
+            events_simulated = MagicMcOrigEvents(
+                az_tel = data.az_tel[selection],
+                alt_tel = data.alt_tel[selection],
+                energy = data.energy[selection],
+                src_x = data.src_x[selection],
+                src_y = data.src_y[selection],
+                file_name = data.file_name
+            )
+
+            data = self.events_triggered
+            selection = pix_ids_triggered == pix_id
+            events_triggered = MagicStereoEvents(
+                reco_alt = data.reco_alt[selection],
+                reco_az = data.reco_az[selection],
+                reco_ra = data.reco_ra[selection],
+                reco_dec = data.reco_dec[selection],
+                reco_energy = data.reco_energy[selection],
+                src_x = data.src_x[selection],
+                src_y = data.src_y[selection],
+                ra_tel = data.ra_tel[selection],
+                dec_tel = data.dec_tel[selection],
+                az_tel = data.az_tel[selection],
+                alt_tel = data.alt_tel[selection],
+                delta_t = data.delta_t[selection],
+                gammaness = data.gammaness[selection],
+                mjd = data.mjd[selection] if len(data.mjd) else data.mjd,
+                mc_alt = data.mc_alt[selection] if len(data.mc_alt) else data.mc_alt,
+                mc_az = data.mc_az[selection] if len(data.mc_az) else data.mc_az,
+                mc_energy = data.mc_energy[selection] if len(data.mc_energy) else data.mc_energy,
+                file_name = self.events_triggered.file_name
+            )
+
+            files.append(
+                MagicMcFile(
+                    file_name = self.file_name,
+                    meta = meta,
+                    events_triggered = events_triggered,
+                    events_simulated = events_simulated,
+                )
+            )
+
+        return files
\ No newline at end of file

From 6c201b385a25dbe30a3ad9cbc8dd2e64e4219fdb Mon Sep 17 00:00:00 2001
From: Ievgen Vovk <vovk@icrr.u-tokyo.ac.jp>
Date: Thu, 19 Oct 2023 12:52:02 +0900
Subject: [PATCH 04/33] MagicMcFile: string representation.

---
 src/srcsim/magic/mc/file.py | 15 +++++++++++++++
 1 file changed, 15 insertions(+)

diff --git a/src/srcsim/magic/mc/file.py b/src/srcsim/magic/mc/file.py
index 48827a2..1caa435 100644
--- a/src/srcsim/magic/mc/file.py
+++ b/src/srcsim/magic/mc/file.py
@@ -13,6 +13,21 @@ class MagicMcFile:
     events_triggered: MagicStereoEvents
     events_simulated: MagicMcOrigEvents
 
+    def __str__(self) -> str:
+        summary = \
+f"""{type(self).__name__} instance
+    {'File name':.<20s}: {self.file_name}
+    {'Simulated events':.<20s}: {self.events_simulated.n_events}
+    {'Triggered events':.<20s}: {self.events_triggered.n_events}
+    {'Energy range':.<20s}: {self.meta.energy_range_min} - {self.meta.energy_range_max}
+    {'Spectral index':.<20s}: {self.meta.spectral_index}
+    {'Scatter range':.<20s}: {self.meta.min_scatter_range} - {self.meta.max_scatter_range}
+    {'Viewcone':.<20s}: {self.meta.min_viewcone_radius} - {self.meta.max_viewcone_radius}
+    {'Telescope azimuth':.<20s}: {self.events_simulated.az_tel.to('deg').min()} - {self.events_simulated.az_tel.to('deg').max()}
+    {'Telescope altitude':.<20s}: {self.events_simulated.alt_tel.to('deg').min()} - {self.events_simulated.alt_tel.to('deg').max()}
+"""
+        return summary
+
     @classmethod
     def from_file(cls, file_name):
         self = MagicMcFile(

From 5367c6a29195f2b6a85c98e1fa5f620e333ff516 Mon Sep 17 00:00:00 2001
From: Ievgen Vovk <vovk@icrr.u-tokyo.ac.jp>
Date: Thu, 19 Oct 2023 13:25:28 +0900
Subject: [PATCH 05/33] Bug fix: MagicMcInfo.num_showers is expected to be a
 quantity.

---
 src/srcsim/magic/mc/file.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/srcsim/magic/mc/file.py b/src/srcsim/magic/mc/file.py
index 1caa435..c29ff43 100644
--- a/src/srcsim/magic/mc/file.py
+++ b/src/srcsim/magic/mc/file.py
@@ -68,7 +68,7 @@ def healpy_split(self, nside):
         files = []
 
         for pix_id in range(npix):
-            n_simulated = numpy.sum(pix_ids_simulated == pix_id)
+            n_simulated = u.one * numpy.sum(pix_ids_simulated == pix_id)
             meta = MagicMcInfo(
                 num_showers = n_simulated,
                 energy_range_min = self.meta.energy_range_min,

From b213ab61be3f9544fed4ef1686b09a0c72337aac Mon Sep 17 00:00:00 2001
From: Ievgen Vovk <vovk@icrr.u-tokyo.ac.jp>
Date: Thu, 19 Oct 2023 13:27:56 +0900
Subject: [PATCH 06/33] Missing import.

---
 src/srcsim/magic/mc/file.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/srcsim/magic/mc/file.py b/src/srcsim/magic/mc/file.py
index c29ff43..1b662d7 100644
--- a/src/srcsim/magic/mc/file.py
+++ b/src/srcsim/magic/mc/file.py
@@ -1,5 +1,6 @@
 import numpy
 import healpy
+import astropy.units as u
 from dataclasses import dataclass
 
 from .info import MagicMcInfo

From d85cc922c85b90c708bff4c1775ef8dc57d7b09e Mon Sep 17 00:00:00 2001
From: Ievgen Vovk <vovk@icrr.u-tokyo.ac.jp>
Date: Thu, 19 Oct 2023 13:30:46 +0900
Subject: [PATCH 07/33] MAGIC to DL2 coverter.

---
 pyproject.toml                     |  1 +
 src/srcsim/scripts/magic_to_dl2.py | 90 ++++++++++++++++++++++++++++++
 2 files changed, 91 insertions(+)
 create mode 100644 src/srcsim/scripts/magic_to_dl2.py

diff --git a/pyproject.toml b/pyproject.toml
index b5440dc..44c91b9 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -35,3 +35,4 @@ dependencies = [
 [project.scripts]
 getruns = "srcsim.scripts.get_runs:main"
 simrun = "srcsim.scripts.sim_run:main"
+magic_to_dl2 = "srcsim.scripts.magic_to_dl2:main"
diff --git a/src/srcsim/scripts/magic_to_dl2.py b/src/srcsim/scripts/magic_to_dl2.py
new file mode 100644
index 0000000..1b3a569
--- /dev/null
+++ b/src/srcsim/scripts/magic_to_dl2.py
@@ -0,0 +1,90 @@
+import argparse
+import logging
+import healpy
+
+from srcsim.magic.mc import MagicMcFile
+
+
+def main():
+    placeholder = "[HEALPY]"
+
+    arg_parser = argparse.ArgumentParser(
+        description="""
+        Convert MAGIC Melibea/Superstar MC files to the DL2 (HDF5) format.
+        """
+    )
+
+    arg_parser.add_argument(
+        "-i",
+        "--input-file", 
+        help='MAGIC MC file in Root format'
+    )
+    arg_parser.add_argument(
+        "-n",
+        "--nside",
+        type=int,
+        default=0,
+        help='NSIDE to use in HEALPix split of the input file.' \
+            ' If nside=0, no spliting is applied. Defaults to 0.'
+    )
+    arg_parser.add_argument(
+        "-o",
+        "--output-file", 
+        help='Output HDF5 file name.' \
+            f' If nside > 0 is given, must containt a "{placeholder}" placeholder,' \
+            ' that will be replaced with "healpy-nside-pixid" when writing the file.'
+    )
+    arg_parser.add_argument(
+        '-v',
+        "--verbose",
+        action='store_true',
+        help='extra verbosity'
+    )
+    args = arg_parser.parse_args()
+
+    logging.basicConfig(
+        format='%(asctime)s %(name)-10s : %(levelname)-8s %(message)s',
+        datefmt='%Y-%m-%dT%H:%M:%S',
+    )
+
+    log = logging.getLogger(__name__)
+
+    if args.verbose:
+        log.setLevel(logging.DEBUG)
+    else:
+        log.setLevel(logging.INFO)
+
+    log.info(f'reading MC file')
+    mcfile = MagicMcFile.from_file(args.input_file)
+    log.debug(str(mcfile))
+
+    if args.nside > 0:
+        if not placeholder in args.output_file:
+            raise ValueError(
+                f'nside > 0 specified, but no "{placeholder}" given in output name "{args.output_file}"'
+            )
+        
+        log.info('running HEALPix spliting')
+        files = mcfile.healpy_split(args.nside)
+
+        npix = healpy.nside2npix(args.nside)
+        ids = range(npix)
+        with_events = tuple(
+            filter(lambda x: x[1].events_simulated.n_events > 0, zip(ids, files))
+        )
+        log.info(f'split complete: {len(with_events)} non-empty files (npix={npix})')
+
+        log.info(f'writing the DL2 files')
+        for pix_id, mfile in with_events:
+            output_name = args.output_file.replace(placeholder, f'healpy-{args.nside}-{pix_id}')
+            log.debug(f'writing {output_name}')
+            mfile.write(output_name)
+    else:
+        log.info(f'writing the DL2 file')
+        mcfile.write(args.output_file)
+
+    log.info('conversion complete')
+
+
+if __name__ == '__main__':
+    main()

From e6ca653f0b8a1e676746c97c25b596480aaf820b Mon Sep 17 00:00:00 2001
From: Ievgen Vovk <vovk@icrr.u-tokyo.ac.jp>
Date: Thu, 19 Oct 2023 13:37:45 +0900
Subject: [PATCH 08/33] MAGIC to DL2 converter: resolution notification.

---
 src/srcsim/scripts/magic_to_dl2.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/src/srcsim/scripts/magic_to_dl2.py b/src/srcsim/scripts/magic_to_dl2.py
index 1b3a569..5db1f6f 100644
--- a/src/srcsim/scripts/magic_to_dl2.py
+++ b/src/srcsim/scripts/magic_to_dl2.py
@@ -1,6 +1,7 @@
 import argparse
 import logging
 import healpy
+import numpy
 
 from srcsim.magic.mc import MagicMcFile
 
@@ -65,6 +66,10 @@ def main():
             )
         
         log.info('running HEALPix spliting')
+        log.info(
+            f'approximate pixel resolution: {numpy.degrees(healpy.nside2resol(args.nside)):.1f} deg' \
+            f' (nside={args.nside})'
+        )
         files = mcfile.healpy_split(args.nside)
 
         npix = healpy.nside2npix(args.nside)

From 2d3f6d6207e53e6f7def765a5c2361c50523c12e Mon Sep 17 00:00:00 2001
From: Ievgen Vovk <vovk@icrr.u-tokyo.ac.jp>
Date: Thu, 19 Oct 2023 14:12:45 +0900
Subject: [PATCH 09/33] Fixing the conversion to quantities, likely broken by
 tables / pandas update.

---
 src/srcsim/mc.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/srcsim/mc.py b/src/srcsim/mc.py
index c74eca3..f9b1dd9 100644
--- a/src/srcsim/mc.py
+++ b/src/srcsim/mc.py
@@ -55,7 +55,7 @@ def __init__(self, file_name=None, obs_id=None, data_table=None, config_table=No
         self.evt_energy = self.data_table['mc_energy'].to_numpy() * self.units['energy']
 
         # Filtering out events with excessive offsets (e.g. due to the simulation numerical accuracy)
-        offset_min, offset_max = self.config_table[['min_viewcone_radius', 'max_viewcone_radius']].iloc[0] * self.units['viewcone']
+        offset_min, offset_max = self.config_table[['min_viewcone_radius', 'max_viewcone_radius']].iloc[0].tolist() * self.units['viewcone']
         evt_offset = self.evt_coord.separation(self.tel_pos)
 
         in_fov = (evt_offset >= offset_min) & (evt_offset <= offset_max)
@@ -128,7 +128,7 @@ def dnde(self, energy):
         return power_law(energy, **self.spec_data)
     
     def dndo(self, coord):
-        offset_min, offset_max = self.config_table[['min_viewcone_radius', 'max_viewcone_radius']].iloc[0] * self.units['viewcone']
+        offset_min, offset_max = self.config_table[['min_viewcone_radius', 'max_viewcone_radius']].iloc[0].tolist() * self.units['viewcone']
         sky_area = 2 * np.pi * (np.cos(offset_min) - np.cos(offset_max)) * u.sr
         norm = 1 / sky_area
 

From 1639f722f40d48e9dd7599c6b010fcb603795773 Mon Sep 17 00:00:00 2001
From: Ievgen Vovk <vovk@icrr.u-tokyo.ac.jp>
Date: Thu, 19 Oct 2023 15:25:01 +0900
Subject: [PATCH 10/33] MagicStereoEvents: added overlooked reconstructed
 position x/y.

---
 src/srcsim/magic/mc/events.py | 16 ++++++++++++++++
 1 file changed, 16 insertions(+)

diff --git a/src/srcsim/magic/mc/events.py b/src/srcsim/magic/mc/events.py
index 7797302..9d18ab1 100644
--- a/src/srcsim/magic/mc/events.py
+++ b/src/srcsim/magic/mc/events.py
@@ -7,6 +7,8 @@
 
 @dataclass(frozen=True)
 class EventSample:
+    reco_src_x: u.Quantity
+    reco_src_y: u.Quantity
     reco_alt: u.Quantity
     reco_az: u.Quantity
     reco_ra: u.Quantity
@@ -65,6 +67,7 @@ def is_mc(self):
 
     @classmethod
     def from_file(cls, file_name):
+        cam2angle = 1.5 / 450 *u.Unit('deg/mm')
         north_offset = 7
 
         event_data = dict()
@@ -83,6 +86,8 @@ def from_file(cls, file_name):
             'MRawEvtHeader_1.fDAQEvtNumber',
             'MSrcPosCam_1.fX',
             'MSrcPosCam_1.fY',
+            'MStereoParDisp.fDirectionX',
+            'MStereoParDisp.fDirectionY',
             'MStereoParDisp.fDirectionRA',
             'MStereoParDisp.fDirectionDec',
             'MStereoParDisp.fDirectionAz',
@@ -108,6 +113,8 @@ def from_file(cls, file_name):
             'MPointingPos_1.fDec': 'dec_tel',
             'MSrcPosCam_1.fX': 'src_x',
             'MSrcPosCam_1.fY': 'src_y',
+            'MStereoParDisp.fDirectionX': 'reco_src_x',
+            'MStereoParDisp.fDirectionY': 'reco_src_y',
             'MStereoParDisp.fDirectionRA': 'reco_ra',
             'MStereoParDisp.fDirectionDec': 'reco_dec',
             'MStereoParDisp.fDirectionAz': 'reco_az',
@@ -123,6 +130,8 @@ def from_file(cls, file_name):
             'delta_t': u.s,
             'src_x': u.mm,
             'src_y': u.mm,
+            'reco_src_x': u.deg,
+            'reco_src_y': u.deg,
             'reco_alt': u.deg,
             'reco_az': u.deg,
             'reco_zd': u.deg,
@@ -186,6 +195,9 @@ def from_file(cls, file_name):
         event_data['alt_tel'] = 90 - event_data['zd_tel']
         event_data['reco_alt'] = 90 - event_data['reco_zd']
         event_data['gammaness'] = 1 - event_data['hadronness']
+
+        event_data['reco_src_x'] = event_data['reco_src_x'] / cam2angle
+        event_data['reco_src_y'] = event_data['reco_src_y'] / cam2angle
         
         for key in event_data:
             event_data[key] = event_data[key].to_numpy()
@@ -198,6 +210,8 @@ def from_file(cls, file_name):
 
         if is_mc:
             events = MagicStereoEvents(
+                reco_src_x = event_data['reco_src_x'],
+                reco_src_y = event_data['reco_src_y'],
                 reco_alt = event_data['reco_alt'],
                 reco_az = event_data['reco_az'],
                 reco_ra = event_data['reco_ra'],
@@ -218,6 +232,8 @@ def from_file(cls, file_name):
             )
         else:
             events = MagicStereoEvents(
+                reco_src_x = event_data['reco_src_x'],
+                reco_src_y = event_data['reco_src_y'],
                 reco_alt = event_data['reco_alt'],
                 reco_az = event_data['reco_az'],
                 reco_ra = event_data['reco_ra'],

From 890de42fb8bf2449f78139b9dcc946c57a3f7882 Mon Sep 17 00:00:00 2001
From: Ievgen Vovk <vovk@icrr.u-tokyo.ac.jp>
Date: Thu, 19 Oct 2023 15:43:46 +0900
Subject: [PATCH 11/33] MagicStereoEvents: fixing reco_src_x/y conversion.

---
 src/srcsim/magic/mc/events.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/srcsim/magic/mc/events.py b/src/srcsim/magic/mc/events.py
index 9d18ab1..e2c5d42 100644
--- a/src/srcsim/magic/mc/events.py
+++ b/src/srcsim/magic/mc/events.py
@@ -195,9 +195,6 @@ def from_file(cls, file_name):
         event_data['alt_tel'] = 90 - event_data['zd_tel']
         event_data['reco_alt'] = 90 - event_data['reco_zd']
         event_data['gammaness'] = 1 - event_data['hadronness']
-
-        event_data['reco_src_x'] = event_data['reco_src_x'] / cam2angle
-        event_data['reco_src_y'] = event_data['reco_src_y'] / cam2angle
         
         for key in event_data:
             event_data[key] = event_data[key].to_numpy()
@@ -206,6 +203,9 @@ def from_file(cls, file_name):
             if key in data_units:
                 event_data[key] = event_data[key] * data_units[key]
 
+        event_data['reco_src_x'] /= cam2angle
+        event_data['reco_src_y'] /= cam2angle
+
         is_mc = 'mc_energy' in event_data
 
         if is_mc:

From 85636bda2d07e775c881dc4ecd6246cc21c99073 Mon Sep 17 00:00:00 2001
From: Ievgen Vovk <vovk@icrr.u-tokyo.ac.jp>
Date: Thu, 19 Oct 2023 15:45:09 +0900
Subject: [PATCH 12/33] MagicStereoEvents: optional MC telescope coordinates.

---
 src/srcsim/magic/mc/events.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/src/srcsim/magic/mc/events.py b/src/srcsim/magic/mc/events.py
index e2c5d42..49ce029 100644
--- a/src/srcsim/magic/mc/events.py
+++ b/src/srcsim/magic/mc/events.py
@@ -26,6 +26,8 @@ class EventSample:
     mc_alt: u.Quantity = u.Quantity([], unit='deg')
     mc_az: u.Quantity = u.Quantity([], unit='deg')
     mc_energy: u.Quantity = u.Quantity([], unit='TeV')
+    mc_alt_tel: u.Quantity = u.Quantity([], unit='deg')
+    mc_az_tel: u.Quantity = u.Quantity([], unit='deg')
 
 
 @dataclass(frozen=True)
@@ -228,6 +230,8 @@ def from_file(cls, file_name):
                 mc_alt = event_data['mc_alt'],
                 mc_az = event_data['mc_az'],
                 mc_energy = event_data['mc_energy'],
+                mc_alt_tel = event_data['alt_tel'],
+                mc_az_tel = event_data['az_tel'],
                 file_name = file_name
             )
         else:

From 2d3bcc47b14e4b9c6cebfbc0e5292f29d9859863 Mon Sep 17 00:00:00 2001
From: Ievgen Vovk <vovk@icrr.u-tokyo.ac.jp>
Date: Thu, 19 Oct 2023 15:45:39 +0900
Subject: [PATCH 13/33] MagicStereoEvents: adding new fields to the export data
 frame.

---
 src/srcsim/magic/mc/events.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/src/srcsim/magic/mc/events.py b/src/srcsim/magic/mc/events.py
index 49ce029..354d7f6 100644
--- a/src/srcsim/magic/mc/events.py
+++ b/src/srcsim/magic/mc/events.py
@@ -262,6 +262,8 @@ def to_df(self):
             # delta_t = u.s,
             src_x = u.mm,
             src_y = u.mm,
+            reco_src_x = u.mm,
+            reco_src_y = u.mm,
             reco_alt = u.deg,
             reco_az = u.deg,
             reco_ra = u.deg,
@@ -275,6 +277,8 @@ def to_df(self):
             dec_tel =u.deg,
             mc_alt =u.deg,
             mc_az =u.deg,
+            mc_alt_tel = u.deg,
+            mc_az_tel = u.deg,
         )
         data = {
             key: self.__getattribute__(key).to(units[key]).value

From b0bab4491018ab5235f5b9c7f347ace2c7edbdb9 Mon Sep 17 00:00:00 2001
From: Ievgen Vovk <vovk@icrr.u-tokyo.ac.jp>
Date: Thu, 19 Oct 2023 15:50:33 +0900
Subject: [PATCH 14/33] MagicStereoEvents: aligning the export units with
 expectation of the MCSample class.

---
 src/srcsim/magic/mc/events.py | 34 +++++++++++++++++-----------------
 1 file changed, 17 insertions(+), 17 deletions(-)

diff --git a/src/srcsim/magic/mc/events.py b/src/srcsim/magic/mc/events.py
index 354d7f6..a4e20f7 100644
--- a/src/srcsim/magic/mc/events.py
+++ b/src/srcsim/magic/mc/events.py
@@ -260,25 +260,25 @@ def from_file(cls, file_name):
     def to_df(self):
         units = dict(
             # delta_t = u.s,
-            src_x = u.mm,
-            src_y = u.mm,
-            reco_src_x = u.mm,
-            reco_src_y = u.mm,
-            reco_alt = u.deg,
-            reco_az = u.deg,
-            reco_ra = u.deg,
-            reco_dec = u.deg,
-            reco_energy = u.GeV,
+            src_x = u.m,
+            src_y = u.m,
+            reco_src_x = u.m,
+            reco_src_y = u.m,
+            reco_alt = u.rad,
+            reco_az = u.rad,
+            reco_ra = u.rad,
+            reco_dec = u.rad,
+            reco_energy = u.TeV,
             gammaness = u.one,
             mjd = u.d,
-            alt_tel = u.deg,
-            az_tel = u.deg,
-            ra_tel = u.deg,
-            dec_tel =u.deg,
-            mc_alt =u.deg,
-            mc_az =u.deg,
-            mc_alt_tel = u.deg,
-            mc_az_tel = u.deg,
+            alt_tel = u.rad,
+            az_tel = u.rad,
+            ra_tel = u.rad,
+            dec_tel =u.rad,
+            mc_alt =u.rad,
+            mc_az = u.rad,
+            mc_alt_tel = u.rad,
+            mc_az_tel = u.rad,
         )
         data = {
             key: self.__getattribute__(key).to(units[key]).value

From 57070eaf0f4a4a5a9ae388eeea7ad3532175c6e4 Mon Sep 17 00:00:00 2001
From: Ievgen Vovk <vovk@icrr.u-tokyo.ac.jp>
Date: Thu, 19 Oct 2023 15:52:34 +0900
Subject: [PATCH 15/33] MagicStereoEvents: fixing MC energy export.

---
 src/srcsim/magic/mc/events.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/srcsim/magic/mc/events.py b/src/srcsim/magic/mc/events.py
index a4e20f7..9e3971d 100644
--- a/src/srcsim/magic/mc/events.py
+++ b/src/srcsim/magic/mc/events.py
@@ -151,6 +151,7 @@ def from_file(cls, file_name):
             'mc_alt':u.deg,
             'mc_az':u.deg,
             'mc_zd':u.deg,
+            'mc_energy': u.GeV,
         }
 
         with uproot.open(file_name) as input_file:
@@ -277,6 +278,7 @@ def to_df(self):
             dec_tel =u.rad,
             mc_alt =u.rad,
             mc_az = u.rad,
+            mc_energy = u.TeV,
             mc_alt_tel = u.rad,
             mc_az_tel = u.rad,
         )

From 9e95e88e6477ebd74596f5f5006ded5ec2feee2c Mon Sep 17 00:00:00 2001
From: Ievgen Vovk <vovk@icrr.u-tokyo.ac.jp>
Date: Thu, 19 Oct 2023 16:04:10 +0900
Subject: [PATCH 16/33] MagicStereoEvents: event_id as DAQ event number.

---
 src/srcsim/magic/mc/events.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/src/srcsim/magic/mc/events.py b/src/srcsim/magic/mc/events.py
index 9e3971d..41146b6 100644
--- a/src/srcsim/magic/mc/events.py
+++ b/src/srcsim/magic/mc/events.py
@@ -7,6 +7,7 @@
 
 @dataclass(frozen=True)
 class EventSample:
+    event_id: u.Quantity
     reco_src_x: u.Quantity
     reco_src_y: u.Quantity
     reco_alt: u.Quantity
@@ -129,6 +130,7 @@ def from_file(cls, file_name):
         }
         
         data_units = {
+            'daq_event_number': u.one,
             'delta_t': u.s,
             'src_x': u.mm,
             'src_y': u.mm,
@@ -213,6 +215,7 @@ def from_file(cls, file_name):
 
         if is_mc:
             events = MagicStereoEvents(
+                event_id = event_data['daq_event_number'],
                 reco_src_x = event_data['reco_src_x'],
                 reco_src_y = event_data['reco_src_y'],
                 reco_alt = event_data['reco_alt'],
@@ -261,6 +264,7 @@ def from_file(cls, file_name):
     def to_df(self):
         units = dict(
             # delta_t = u.s,
+            event_id = u.one,
             src_x = u.m,
             src_y = u.m,
             reco_src_x = u.m,

From 90d71fb7fba1e7caa646fbd0972802ffb2a0a13d Mon Sep 17 00:00:00 2001
From: Ievgen Vovk <vovk@icrr.u-tokyo.ac.jp>
Date: Thu, 19 Oct 2023 16:13:29 +0900
Subject: [PATCH 17/33] MagicMcInfo: reading run number.

---
 src/srcsim/magic/mc/info.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/src/srcsim/magic/mc/info.py b/src/srcsim/magic/mc/info.py
index 82f24b7..29ba8f1 100644
--- a/src/srcsim/magic/mc/info.py
+++ b/src/srcsim/magic/mc/info.py
@@ -6,6 +6,7 @@
 
 @dataclass(frozen=True)
 class McInfo:
+    run_number: u.Quantity
     num_showers: u.Quantity
     energy_range_min: u.Quantity
     energy_range_max: u.Quantity
@@ -41,6 +42,9 @@ def __repr__(self):
     def from_file(cls, file_name):
         meta = dict()
         with uproot.open(file_name) as input_file:
+            meta['run_number'] = u.one * int(
+                input_file['RunHeaders']['MMcRunHeader_1./MMcRunHeader_1.fMcRunNumber'].array()[0]
+            )
             meta['num_showers'] = u.one * int(
                 input_file['RunHeaders']['MMcRunHeader_1./MMcRunHeader_1.fNumSimulatedShowers'].array()[0]
             )
@@ -53,6 +57,7 @@ def from_file(cls, file_name):
             meta['min_viewcone_radius'] = u.deg * 0
 
         info = MagicMcInfo(
+            run_number = meta['run_number'],
             num_showers = meta['num_showers'],
             energy_range_min = meta['energy_range_min'],
             energy_range_max = meta['energy_range_max'],

From a9b8e216c6bd977c743009a404e05b10114b50a8 Mon Sep 17 00:00:00 2001
From: Ievgen Vovk <vovk@icrr.u-tokyo.ac.jp>
Date: Thu, 19 Oct 2023 16:17:13 +0900
Subject: [PATCH 18/33] MagicMcInfo: using CTA naming style.

---
 src/srcsim/magic/mc/info.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/srcsim/magic/mc/info.py b/src/srcsim/magic/mc/info.py
index 29ba8f1..34b0b93 100644
--- a/src/srcsim/magic/mc/info.py
+++ b/src/srcsim/magic/mc/info.py
@@ -6,7 +6,7 @@
 
 @dataclass(frozen=True)
 class McInfo:
-    run_number: u.Quantity
+    obs_id: u.Quantity
     num_showers: u.Quantity
     energy_range_min: u.Quantity
     energy_range_max: u.Quantity
@@ -42,7 +42,7 @@ def __repr__(self):
     def from_file(cls, file_name):
         meta = dict()
         with uproot.open(file_name) as input_file:
-            meta['run_number'] = u.one * int(
+            meta['obs_id'] = u.one * int(
                 input_file['RunHeaders']['MMcRunHeader_1./MMcRunHeader_1.fMcRunNumber'].array()[0]
             )
             meta['num_showers'] = u.one * int(
@@ -57,7 +57,7 @@ def from_file(cls, file_name):
             meta['min_viewcone_radius'] = u.deg * 0
 
         info = MagicMcInfo(
-            run_number = meta['run_number'],
+            obs_id = meta['obs_id'],
             num_showers = meta['num_showers'],
             energy_range_min = meta['energy_range_min'],
             energy_range_max = meta['energy_range_max'],

From d01d1fc362bd83097e71d9cf92d9673239443f69 Mon Sep 17 00:00:00 2001
From: Ievgen Vovk <vovk@icrr.u-tokyo.ac.jp>
Date: Thu, 19 Oct 2023 16:25:15 +0900
Subject: [PATCH 19/33] MagicStereoEvents: adding obs_id.

---
 src/srcsim/magic/mc/events.py | 12 ++++++++++++
 1 file changed, 12 insertions(+)

diff --git a/src/srcsim/magic/mc/events.py b/src/srcsim/magic/mc/events.py
index 41146b6..44d1b66 100644
--- a/src/srcsim/magic/mc/events.py
+++ b/src/srcsim/magic/mc/events.py
@@ -7,6 +7,7 @@
 
 @dataclass(frozen=True)
 class EventSample:
+    obs_id: u.Quantity
     event_id: u.Quantity
     reco_src_x: u.Quantity
     reco_src_y: u.Quantity
@@ -130,6 +131,7 @@ def from_file(cls, file_name):
         }
         
         data_units = {
+            'run_number': u.one,
             'daq_event_number': u.one,
             'delta_t': u.s,
             'src_x': u.mm,
@@ -190,6 +192,13 @@ def from_file(cls, file_name):
 
                     event_data['mjd'] = mjd + (millisec / 1e3 + nanosec / 1e9) / 86400.0
 
+                event_data['run_number'] = numpy.repeat(
+                    int(
+                        input_file['RunHeaders']['MRawRunHeader_1./MRawRunHeader_1.fRunNumber'].array()[0]
+                    ),
+                    len(event_data['alt_tel'])
+                )
+
             else:
                 # The file is likely corrupted, so return empty arrays
                 for key in names_mapping:
@@ -215,6 +224,7 @@ def from_file(cls, file_name):
 
         if is_mc:
             events = MagicStereoEvents(
+                obs_id = event_data['run_number'],
                 event_id = event_data['daq_event_number'],
                 reco_src_x = event_data['reco_src_x'],
                 reco_src_y = event_data['reco_src_y'],
@@ -240,6 +250,7 @@ def from_file(cls, file_name):
             )
         else:
             events = MagicStereoEvents(
+                obs_id = event_data['run_number'],
                 reco_src_x = event_data['reco_src_x'],
                 reco_src_y = event_data['reco_src_y'],
                 reco_alt = event_data['reco_alt'],
@@ -264,6 +275,7 @@ def from_file(cls, file_name):
     def to_df(self):
         units = dict(
             # delta_t = u.s,
+            obs_id = u.one,
             event_id = u.one,
             src_x = u.m,
             src_y = u.m,

From 95437213d6133d98d0cfa847e017928b926534aa Mon Sep 17 00:00:00 2001
From: Ievgen Vovk <vovk@icrr.u-tokyo.ac.jp>
Date: Thu, 19 Oct 2023 16:25:52 +0900
Subject: [PATCH 20/33] Overlooked event_id.

---
 src/srcsim/magic/mc/events.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/srcsim/magic/mc/events.py b/src/srcsim/magic/mc/events.py
index 44d1b66..14b98d1 100644
--- a/src/srcsim/magic/mc/events.py
+++ b/src/srcsim/magic/mc/events.py
@@ -251,6 +251,7 @@ def from_file(cls, file_name):
         else:
             events = MagicStereoEvents(
                 obs_id = event_data['run_number'],
+                event_id = event_data['daq_event_number'],
                 reco_src_x = event_data['reco_src_x'],
                 reco_src_y = event_data['reco_src_y'],
                 reco_alt = event_data['reco_alt'],

From 49fd09b162a08d44ab5fb73909c0d247ec4e2f05 Mon Sep 17 00:00:00 2001
From: Ievgen Vovk <vovk@icrr.u-tokyo.ac.jp>
Date: Thu, 19 Oct 2023 16:35:23 +0900
Subject: [PATCH 21/33] Bug fix in run number reading (could not be cast to
 numpy with "to_numpy()" otherwise).

---
 src/srcsim/magic/mc/events.py | 12 +++++-------
 1 file changed, 5 insertions(+), 7 deletions(-)

diff --git a/src/srcsim/magic/mc/events.py b/src/srcsim/magic/mc/events.py
index 14b98d1..3cec07a 100644
--- a/src/srcsim/magic/mc/events.py
+++ b/src/srcsim/magic/mc/events.py
@@ -192,18 +192,14 @@ def from_file(cls, file_name):
 
                     event_data['mjd'] = mjd + (millisec / 1e3 + nanosec / 1e9) / 86400.0
 
-                event_data['run_number'] = numpy.repeat(
-                    int(
-                        input_file['RunHeaders']['MRawRunHeader_1./MRawRunHeader_1.fRunNumber'].array()[0]
-                    ),
-                    len(event_data['alt_tel'])
-                )
+                run_number = input_file['RunHeaders']['MRawRunHeader_1./MRawRunHeader_1.fRunNumber'].array()[0]
 
             else:
                 # The file is likely corrupted, so return empty arrays
                 for key in names_mapping:
                     name = names_mapping[key]
                     event_data[name] = numpy.zeros(0)
+                run_number = 0
                     
         event_data['mc_alt'] = 90 - event_data['mc_zd']
         event_data['alt_tel'] = 90 - event_data['zd_tel']
@@ -212,7 +208,9 @@ def from_file(cls, file_name):
         
         for key in event_data:
             event_data[key] = event_data[key].to_numpy()
-            
+
+        event_data['run_number'] = numpy.repeat(run_number, len(event_data['az_tel']))
+
         for key in event_data:
             if key in data_units:
                 event_data[key] = event_data[key] * data_units[key]

From 9449122044f125f0899796492a68f5df93de02f3 Mon Sep 17 00:00:00 2001
From: Ievgen Vovk <vovk@icrr.u-tokyo.ac.jp>
Date: Thu, 19 Oct 2023 16:40:17 +0900
Subject: [PATCH 22/33] Adding new fields to the files following the HEALPix
 split.

---
 src/srcsim/magic/mc/file.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/src/srcsim/magic/mc/file.py b/src/srcsim/magic/mc/file.py
index 1b662d7..70d4a2f 100644
--- a/src/srcsim/magic/mc/file.py
+++ b/src/srcsim/magic/mc/file.py
@@ -71,6 +71,7 @@ def healpy_split(self, nside):
         for pix_id in range(npix):
             n_simulated = u.one * numpy.sum(pix_ids_simulated == pix_id)
             meta = MagicMcInfo(
+                obs_id = self.meta.obs_id,
                 num_showers = n_simulated,
                 energy_range_min = self.meta.energy_range_min,
                 energy_range_max = self.meta.energy_range_max,
@@ -95,6 +96,10 @@ def healpy_split(self, nside):
             data = self.events_triggered
             selection = pix_ids_triggered == pix_id
             events_triggered = MagicStereoEvents(
+                obs_id = data.obs_id[selection],
+                event_id = data.event_id[selection],
+                reco_src_x = data.reco_src_x[selection],
+                reco_src_y = data.reco_src_y[selection],
                 reco_alt = data.reco_alt[selection],
                 reco_az = data.reco_az[selection],
                 reco_ra = data.reco_ra[selection],

From efff8e506fd1f664e25f9d7fb83805c8b27d59a8 Mon Sep 17 00:00:00 2001
From: Ievgen Vovk <vovk@icrr.u-tokyo.ac.jp>
Date: Thu, 19 Oct 2023 18:31:54 +0900
Subject: [PATCH 23/33] Effective focal length based on tel_id.

---
 src/srcsim/mc.py  | 20 +++++++++++++++-----
 src/srcsim/off.py | 16 +++++++++++++++-
 src/srcsim/run.py |  4 ++--
 3 files changed, 32 insertions(+), 8 deletions(-)

diff --git a/src/srcsim/mc.py b/src/srcsim/mc.py
index f9b1dd9..1a8636d 100644
--- a/src/srcsim/mc.py
+++ b/src/srcsim/mc.py
@@ -11,6 +11,12 @@ def power_law(e, e0, norm, index):
 
 
 class MCSample:
+    # Effective focal lengths
+    camscale = {
+        1: 1 / 29.30565 * u.Unit('rad/m'),
+        5: 1 / 18.2121 * u.Unit('rad/m'),
+    }
+
     def __init__(self, file_name=None, obs_id=None, data_table=None, config_table=None):
         self.units = dict(
             energy = u.TeV,
@@ -18,10 +24,6 @@ def __init__(self, file_name=None, obs_id=None, data_table=None, config_table=No
             distance = u.m,
             viewcone = u.deg
         )
-        
-        # TODO: refine this value
-        lst_focal_length = 28.01 * u.m
-        self.cam2angle = 1 * u.rad / lst_focal_length
 
         # TODO: refine the logic below / implement nicer
         if data_table is not None and config_table is not None:
@@ -49,7 +51,7 @@ def __init__(self, file_name=None, obs_id=None, data_table=None, config_table=No
         self.spec_data = self.get_spec_data(nevents, emin, emax, index)
         self.spec_data['norm'] /= ground_area
 
-        cam_x, cam_y = self.data_table[['src_x', 'src_y']].to_numpy().transpose() * self.units['distance'] * self.cam2angle
+        cam_x, cam_y = self.data_table[['src_x', 'src_y']].to_numpy().transpose() * self.units['distance'] * self.cam2angle(self.data_table['tel_id'].tolist())
         self.evt_coord = SkyCoord(cam_x, cam_y, frame=self.tel_pos.skyoffset_frame())
 
         self.evt_energy = self.data_table['mc_energy'].to_numpy() * self.units['energy']
@@ -106,6 +108,14 @@ def read_config(cls, file_name, obs_id):
 
             return pd.DataFrame(data=data)
 
+    def cam2angle(self, tel_id):
+        cam2angle = u.Quantity(
+            list(
+                map(lambda idx: self.camscale[idx], tel_id)
+            )
+        )
+        return cam2angle
+
     def get_spec_data(self, n_events, emin, emax, index=-1):
         e0 = (emin * emax)**0.5
 
diff --git a/src/srcsim/off.py b/src/srcsim/off.py
index 6ab4533..cc5d357 100644
--- a/src/srcsim/off.py
+++ b/src/srcsim/off.py
@@ -7,6 +7,12 @@
 
 
 class OffSample:
+    # Effective focal lengths
+    camscale = {
+        1: 1 / 29.30565 * u.Unit('rad/m'),
+        5: 1 / 18.2121 * u.Unit('rad/m'),
+    }
+
     def __init__(self, file_name=None, obs_id=None, data_table=None):
         self.units = dict(
             energy = u.TeV,
@@ -43,7 +49,7 @@ def __init__(self, file_name=None, obs_id=None, data_table=None):
             unit=self.units['angle'],
             frame='altaz'
         )
-        cam_x, cam_y = self.data_table[['reco_src_x', 'reco_src_y']].to_numpy().transpose() * self.units['distance'] * self.cam2angle
+        cam_x, cam_y = self.data_table[['reco_src_x', 'reco_src_y']].to_numpy().transpose() * self.units['distance'] * self.cam2angle(self.data_table['tel_id'].tolist())
         self.evt_coord = SkyCoord(cam_x, cam_y, frame=tel_pos.skyoffset_frame())
         self.evt_energy = self.data_table['reco_energy'].to_numpy() * self.units['energy']
         
@@ -72,6 +78,14 @@ def calc_obs_duration(self, data_table):
 
         return t_elapsed
 
+    def cam2angle(self, tel_id):
+        cam2angle = u.Quantity(
+            list(
+                map(lambda idx: self.camscale[idx], tel_id)
+            )
+        )
+        return cam2angle
+
     def dndedo(self, energy, coord):
         dummy = 1 + 0 * energy.value
         val = dummy * 1 / (1 / self.obs_duration * u.Unit('1/(s TeV sr)'))
diff --git a/src/srcsim/run.py b/src/srcsim/run.py
index 46b0abe..45416f2 100644
--- a/src/srcsim/run.py
+++ b/src/srcsim/run.py
@@ -209,8 +209,8 @@ def predict(self, mccollections, source, tel_pos_tolerance=None, time_step=1*u.m
 
                     # Reconstructed events coordinates
                     reco_coords = SkyCoord(
-                        evt['reco_src_x'].to_numpy() * sample.units['distance'] * sample.cam2angle,
-                        evt['reco_src_y'].to_numpy() * sample.units['distance'] * sample.cam2angle,
+                        evt['reco_src_x'].to_numpy() * sample.units['distance'] * sample.cam2angle(evt['tel_id'].to_numpy()),
+                        evt['reco_src_y'].to_numpy() * sample.units['distance'] * sample.cam2angle(evt['tel_id'].to_numpy()),
                         frame=offset_frame
                     )
                     evt = evt.assign(

From d20b2e51d83f89319e864046ca3d59f56b31cb0d Mon Sep 17 00:00:00 2001
From: Ievgen Vovk <vovk@icrr.u-tokyo.ac.jp>
Date: Thu, 19 Oct 2023 18:53:52 +0900
Subject: [PATCH 24/33] Fixing the conversion to quantities, likely broken by
 tables / pandas update.

---
 src/srcsim/mc.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/srcsim/mc.py b/src/srcsim/mc.py
index 1a8636d..9b57038 100644
--- a/src/srcsim/mc.py
+++ b/src/srcsim/mc.py
@@ -42,7 +42,7 @@ def __init__(self, file_name=None, obs_id=None, data_table=None, config_table=No
         self.tel_pos = SkyCoord(pointing_data['mc_az_tel'], pointing_data['mc_alt_tel'], unit=self.units['angle'], frame='altaz')
         
         # Working out the simulation spectrum
-        rmin, rmax = self.config_table[['min_scatter_range', 'max_scatter_range']].iloc[0] * self.units['distance']
+        rmin, rmax = self.config_table[['min_scatter_range', 'max_scatter_range']].iloc[0].tolist() * self.units['distance']
         ground_area = np.pi * (rmax**2 - rmin**2)
         nevents = self.config_table['num_showers'].iloc[0] * self.config_table['shower_reuse'].iloc[0]
         emin = self.config_table['energy_range_min'].iloc[0] * self.units['energy']

From ca6b96829ead712f68464428bc41e833297d2add Mon Sep 17 00:00:00 2001
From: Ievgen Vovk <vovk@icrr.u-tokyo.ac.jp>
Date: Thu, 19 Oct 2023 19:11:18 +0900
Subject: [PATCH 25/33] Added missing fields when splitting an MC file.

---
 src/srcsim/magic/mc/file.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/srcsim/magic/mc/file.py b/src/srcsim/magic/mc/file.py
index 70d4a2f..09866a0 100644
--- a/src/srcsim/magic/mc/file.py
+++ b/src/srcsim/magic/mc/file.py
@@ -117,6 +117,8 @@ def healpy_split(self, nside):
                 mc_alt = data.mc_alt[selection] if len(data.mc_alt) else data.mc_alt,
                 mc_az = data.mc_az[selection] if len(data.mc_az) else data.mc_az,
                 mc_energy = data.mc_energy[selection] if len(data.mc_energy) else data.mc_energy,
+                mc_alt_tel = data.mc_alt_tel[selection] if len(data.mc_alt_tel) else data.mc_alt_tel,
+                mc_az_tel = data.mc_az_tel[selection] if len(data.mc_az_tel) else data.mc_az_tel,
                 file_name = self.events_triggered.file_name
             )
 

From 1704aba35f43c96963c000b63542767defef95de Mon Sep 17 00:00:00 2001
From: Ievgen Vovk <vovk@icrr.u-tokyo.ac.jp>
Date: Tue, 24 Oct 2023 13:12:17 +0900
Subject: [PATCH 26/33] MCSample: when computing event coordinates, using
 per-event telescope position instead of the averaged one. Required for cases
 when telescope postion in MC is not fixed (e.g. for MAGIC).

---
 src/srcsim/mc.py | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/src/srcsim/mc.py b/src/srcsim/mc.py
index 9b57038..06dd558 100644
--- a/src/srcsim/mc.py
+++ b/src/srcsim/mc.py
@@ -38,8 +38,10 @@ def __init__(self, file_name=None, obs_id=None, data_table=None, config_table=No
             self.data_table = pd.read_hdf(file_name, 'dl2/event/telescope/parameters/LST_LSTCam').query(f'obs_id == {obs_id}')
 
         # Getting the telescope pointing
-        pointing_data = self.data_table[['mc_az_tel', 'mc_alt_tel']].mean()
-        self.tel_pos = SkyCoord(pointing_data['mc_az_tel'], pointing_data['mc_alt_tel'], unit=self.units['angle'], frame='altaz')
+        pointing_data = self.data_table[['mc_az_tel', 'mc_alt_tel']]
+        tel_pos = SkyCoord(pointing_data['mc_az_tel'], pointing_data['mc_alt_tel'], unit=self.units['angle'], frame='altaz')
+        # Average poiting
+        self.tel_pos = SkyCoord(tel_pos.az.mean(), tel_pos.alt.mean(), frame='altaz')
         
         # Working out the simulation spectrum
         rmin, rmax = self.config_table[['min_scatter_range', 'max_scatter_range']].iloc[0].tolist() * self.units['distance']
@@ -52,13 +54,13 @@ def __init__(self, file_name=None, obs_id=None, data_table=None, config_table=No
         self.spec_data['norm'] /= ground_area
 
         cam_x, cam_y = self.data_table[['src_x', 'src_y']].to_numpy().transpose() * self.units['distance'] * self.cam2angle(self.data_table['tel_id'].tolist())
-        self.evt_coord = SkyCoord(cam_x, cam_y, frame=self.tel_pos.skyoffset_frame())
+        self.evt_coord = SkyCoord(cam_x, cam_y, frame=tel_pos.skyoffset_frame())
 
         self.evt_energy = self.data_table['mc_energy'].to_numpy() * self.units['energy']
 
         # Filtering out events with excessive offsets (e.g. due to the simulation numerical accuracy)
         offset_min, offset_max = self.config_table[['min_viewcone_radius', 'max_viewcone_radius']].iloc[0].tolist() * self.units['viewcone']
-        evt_offset = self.evt_coord.separation(self.tel_pos)
+        evt_offset = self.evt_coord.separation(tel_pos)
 
         in_fov = (evt_offset >= offset_min) & (evt_offset <= offset_max)
         self.data_table = self.data_table[in_fov]

From 5f351cfaadeeaf420162d0958dea9c9196a7ee4c Mon Sep 17 00:00:00 2001
From: Ievgen Vovk <vovk@icrr.u-tokyo.ac.jp>
Date: Tue, 24 Oct 2023 13:36:07 +0900
Subject: [PATCH 27/33] MagicStereoEvents: using expected MAGIC tel_id.

---
 src/srcsim/magic/mc/events.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/src/srcsim/magic/mc/events.py b/src/srcsim/magic/mc/events.py
index 3cec07a..1ca03b8 100644
--- a/src/srcsim/magic/mc/events.py
+++ b/src/srcsim/magic/mc/events.py
@@ -7,6 +7,7 @@
 
 @dataclass(frozen=True)
 class EventSample:
+    tel_id: u.Quantity
     obs_id: u.Quantity
     event_id: u.Quantity
     reco_src_x: u.Quantity
@@ -73,6 +74,7 @@ def is_mc(self):
     def from_file(cls, file_name):
         cam2angle = 1.5 / 450 *u.Unit('deg/mm')
         north_offset = 7
+        magic_tel_id = 5 * u.one  # tel_id expected from the MAGIC+LST simulations
 
         event_data = dict()
 
@@ -222,6 +224,7 @@ def from_file(cls, file_name):
 
         if is_mc:
             events = MagicStereoEvents(
+                tel_id = numpy.repeat(magic_tel_id, len(event_data['run_number'])),
                 obs_id = event_data['run_number'],
                 event_id = event_data['daq_event_number'],
                 reco_src_x = event_data['reco_src_x'],
@@ -248,6 +251,7 @@ def from_file(cls, file_name):
             )
         else:
             events = MagicStereoEvents(
+                tel_id = numpy.repeat(magic_tel_id, len(event_data['run_number'])),
                 obs_id = event_data['run_number'],
                 event_id = event_data['daq_event_number'],
                 reco_src_x = event_data['reco_src_x'],
@@ -274,6 +278,7 @@ def from_file(cls, file_name):
     def to_df(self):
         units = dict(
             # delta_t = u.s,
+            tel_id = u.one,
             obs_id = u.one,
             event_id = u.one,
             src_x = u.m,

From 70d9e1b6566d36d8997ebea935a279c8c632e464 Mon Sep 17 00:00:00 2001
From: Ievgen Vovk <vovk@icrr.u-tokyo.ac.jp>
Date: Tue, 24 Oct 2023 13:39:58 +0900
Subject: [PATCH 28/33] MagicMcFile: special treatment for run_config, that
 appears to be too complex for Pandas.

---
 src/srcsim/magic/mc/file.py | 49 +++++++++++++++++++++++++++++++++----
 1 file changed, 44 insertions(+), 5 deletions(-)

diff --git a/src/srcsim/magic/mc/file.py b/src/srcsim/magic/mc/file.py
index 09866a0..922aabf 100644
--- a/src/srcsim/magic/mc/file.py
+++ b/src/srcsim/magic/mc/file.py
@@ -3,10 +3,27 @@
 import astropy.units as u
 from dataclasses import dataclass
 
+from tables import open_file
+from tables import IsDescription
+from tables import UInt32Col, UInt64Col, Float32Col
+
 from .info import MagicMcInfo
 from .events import MagicStereoEvents, MagicMcOrigEvents
 
 
+class TablesMcInfo(IsDescription):
+    obs_id = UInt64Col()
+    num_showers = UInt64Col()
+    shower_reuse = UInt32Col()
+    energy_range_min = Float32Col()
+    energy_range_max = Float32Col()
+    spectral_index = Float32Col()
+    max_scatter_range = Float32Col()
+    min_scatter_range = Float32Col()
+    max_viewcone_radius = Float32Col()
+    min_viewcone_radius = Float32Col()
+
+
 @dataclass(frozen=True)
 class MagicMcFile:
     file_name: str
@@ -40,12 +57,34 @@ def from_file(cls, file_name):
         return self
 
     def write(self, file_name):
-        self.meta.to_df().to_hdf(
-            file_name,
-            key='/simulation/run_config',
-            mode='w'
+        units = dict(
+            energy = u.TeV,
+            angle = u.rad,
+            distance = u.m,
+            viewcone = u.deg
         )
-        self.events_triggered.to_df().to_hdf(
+
+        # Special treatment for run config as it's structure may be too complex for Pandas
+        with open_file(file_name, mode="w", title="MAGIC MC file") as output:
+            group = output.create_group("/", 'simulation', 'Simulation information')
+            table = output.create_table(group, 'run_config', TablesMcInfo, 'Simulation run configuration')
+            
+            entry = table.row
+            entry['obs_id'] = self.meta.obs_id.value
+            entry['num_showers'] = self.meta.num_showers.value
+            entry['shower_reuse'] = 1
+            entry['spectral_index'] = self.meta.spectral_index.value
+            entry['energy_range_min'] = self.meta.energy_range_min.to(units['energy']).value
+            entry['energy_range_max'] = self.meta.energy_range_max.to(units['energy']).value
+            entry['max_scatter_range'] = self.meta.max_scatter_range.to(units['distance']).value
+            entry['min_scatter_range'] = self.meta.min_scatter_range.to(units['distance']).value
+            entry['max_viewcone_radius'] = self.meta.max_viewcone_radius.to(units['viewcone']).value
+            entry['min_viewcone_radius'] = self.meta.min_viewcone_radius.to(units['viewcone']).value
+            entry.append()
+            table.flush()
+
+        df = self.events_triggered.to_df()
+        df.to_hdf(
             file_name,
             key='/dl2/event/telescope/parameters/MAGIC_MAGICCam',
             mode='a'

From db97ba55ea5324f9641cb877a132681dd33a0c2e Mon Sep 17 00:00:00 2001
From: Ievgen Vovk <vovk@icrr.u-tokyo.ac.jp>
Date: Tue, 24 Oct 2023 13:41:14 +0900
Subject: [PATCH 29/33] MagicMcFile: overwriting obs_id with the correct value.

---
 src/srcsim/magic/mc/file.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/srcsim/magic/mc/file.py b/src/srcsim/magic/mc/file.py
index 922aabf..2d9ffe8 100644
--- a/src/srcsim/magic/mc/file.py
+++ b/src/srcsim/magic/mc/file.py
@@ -84,6 +84,8 @@ def write(self, file_name):
             table.flush()
 
         df = self.events_triggered.to_df()
+        # Overwriting obs_id as it is not properly stored for events in the MC files
+        df['obs_id'] = self.meta.obs_id.value
         df.to_hdf(
             file_name,
             key='/dl2/event/telescope/parameters/MAGIC_MAGICCam',

From 6a608a0b2bc4c8104d3fcd8344587bbeb66463ce Mon Sep 17 00:00:00 2001
From: Ievgen Vovk <vovk@icrr.u-tokyo.ac.jp>
Date: Tue, 24 Oct 2023 13:42:52 +0900
Subject: [PATCH 30/33] MagicMcFile: pretending camera to be that of LST (to
 enable loading with MCSample; possibly needs to be revised)

---
 src/srcsim/magic/mc/file.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/srcsim/magic/mc/file.py b/src/srcsim/magic/mc/file.py
index 2d9ffe8..b28c19e 100644
--- a/src/srcsim/magic/mc/file.py
+++ b/src/srcsim/magic/mc/file.py
@@ -86,9 +86,10 @@ def write(self, file_name):
         df = self.events_triggered.to_df()
         # Overwriting obs_id as it is not properly stored for events in the MC files
         df['obs_id'] = self.meta.obs_id.value
+        # TODO: check if one can use "MAGIC_MAGICCam" instead
         df.to_hdf(
             file_name,
-            key='/dl2/event/telescope/parameters/MAGIC_MAGICCam',
+            key='/dl2/event/telescope/parameters/LST_LSTCam',
             mode='a'
         )
 

From 9cb962a25d01689480e6553fd988797de72c69c2 Mon Sep 17 00:00:00 2001
From: Ievgen Vovk <vovk@icrr.u-tokyo.ac.jp>
Date: Tue, 24 Oct 2023 13:48:17 +0900
Subject: [PATCH 31/33] Incrementing code version due to the added new
 functionality (MAGIC simulation loading).

---
 pyproject.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyproject.toml b/pyproject.toml
index 44c91b9..6a1fbc5 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "srcsim"
-version = "1.0.0"
+version = "1.1.0"
 authors = [
   { name="Ievgen Vovk", email="vovk@icrr.u-tokyo.ac.jp" },
   { name="Marcel Strzys", email="strzys@icrr.u-tokyo.ac.jp" },

From 9846c09492e9250c8354b860601c0d72c13ac8c1 Mon Sep 17 00:00:00 2001
From: Ievgen Vovk <vovk@icrr.u-tokyo.ac.jp>
Date: Tue, 24 Oct 2023 13:50:05 +0900
Subject: [PATCH 32/33] Relaxing package dependencies.

---
 pyproject.toml | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/pyproject.toml b/pyproject.toml
index 6a1fbc5..fd2ded1 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -20,12 +20,12 @@ classifiers = [
     "Operating System :: OS Independent",
 ]
 dependencies = [
-    "astropy==5.1",
+    "astropy>=5.1",
     "numpy>=1.21.0",
-    "pandas==1.4.1",
+    "pandas>=1.4.1",
     "scipy>=1.5.0",
-    "tables==3.7.0",
-    "PyYAML==5.3.1"
+    "tables>=3.7.0",
+    "PyYAML>=5.3.1"
 ]
 
 [project.urls]

From 2874eb2e0110744952eb5cdf9e6ed1ef9e355690 Mon Sep 17 00:00:00 2001
From: Ievgen Vovk <vovk@icrr.u-tokyo.ac.jp>
Date: Tue, 24 Oct 2023 14:09:34 +0900
Subject: [PATCH 33/33] Revert "MCSample: when computing event coordinates,
 using per-event telescope position instead of the averaged one."

This reverts commit 1704aba35f43c96963c000b63542767defef95de.
Reason: commit breaks the MCSample.dndo() logic, which needs to be
substantially revised before implementing it.
---
 src/srcsim/mc.py | 10 ++++------
 1 file changed, 4 insertions(+), 6 deletions(-)

diff --git a/src/srcsim/mc.py b/src/srcsim/mc.py
index 06dd558..9b57038 100644
--- a/src/srcsim/mc.py
+++ b/src/srcsim/mc.py
@@ -38,10 +38,8 @@ def __init__(self, file_name=None, obs_id=None, data_table=None, config_table=No
             self.data_table = pd.read_hdf(file_name, 'dl2/event/telescope/parameters/LST_LSTCam').query(f'obs_id == {obs_id}')
 
         # Getting the telescope pointing
-        pointing_data = self.data_table[['mc_az_tel', 'mc_alt_tel']]
-        tel_pos = SkyCoord(pointing_data['mc_az_tel'], pointing_data['mc_alt_tel'], unit=self.units['angle'], frame='altaz')
-        # Average poiting
-        self.tel_pos = SkyCoord(tel_pos.az.mean(), tel_pos.alt.mean(), frame='altaz')
+        pointing_data = self.data_table[['mc_az_tel', 'mc_alt_tel']].mean()
+        self.tel_pos = SkyCoord(pointing_data['mc_az_tel'], pointing_data['mc_alt_tel'], unit=self.units['angle'], frame='altaz')
         
         # Working out the simulation spectrum
         rmin, rmax = self.config_table[['min_scatter_range', 'max_scatter_range']].iloc[0].tolist() * self.units['distance']
@@ -54,13 +52,13 @@ def __init__(self, file_name=None, obs_id=None, data_table=None, config_table=No
         self.spec_data['norm'] /= ground_area
 
         cam_x, cam_y = self.data_table[['src_x', 'src_y']].to_numpy().transpose() * self.units['distance'] * self.cam2angle(self.data_table['tel_id'].tolist())
-        self.evt_coord = SkyCoord(cam_x, cam_y, frame=tel_pos.skyoffset_frame())
+        self.evt_coord = SkyCoord(cam_x, cam_y, frame=self.tel_pos.skyoffset_frame())
 
         self.evt_energy = self.data_table['mc_energy'].to_numpy() * self.units['energy']
 
         # Filtering out events with excessive offsets (e.g. due to the simulation numerical accuracy)
         offset_min, offset_max = self.config_table[['min_viewcone_radius', 'max_viewcone_radius']].iloc[0].tolist() * self.units['viewcone']
-        evt_offset = self.evt_coord.separation(tel_pos)
+        evt_offset = self.evt_coord.separation(self.tel_pos)
 
         in_fov = (evt_offset >= offset_min) & (evt_offset <= offset_max)
         self.data_table = self.data_table[in_fov]