Merge pull request #84 from cta-observatory/laser_test

Add reading of LIDAR reports from ROOT files.

.codacy.yml

@@ -6,3 +6,5 @@ engines:
  pyflakes:
    disable:
      - F999
+ assert_used: skips: ['*/*/test_*.py', '*/test_*.py', 'test_*.py']
+

ctapipe_io_magic/__init__.py

@@ -10,14 +10,15 @@
 import uproot
 import logging
 import numpy as np
+from typing import List, Any
 from pathlib import Path
 from decimal import Decimal
 from astropy import units as u
 from astropy.time import Time
 from pkg_resources import resource_filename
 
 from ctapipe.io import EventSource, DataLevel
-from ctapipe.core import Provenance
+from ctapipe.core import Provenance, Container, Field
 from ctapipe.core.traits import Bool, UseEnum
 from ctapipe.coordinates import CameraFrame
 
@@ -92,6 +93,68 @@
     DATA_MAGIC_LST_TRIGGER: EventType.SUBARRAY,
 }
 
+class ReportLaserContainer(Container):
+    """ Container for Magic laser parameters """
+    UniqueID = Field(Any, 'No.')
+    Bits = Field(Any, 'ID')
+    MJD = Field(np.float64, 'Modified Julian Date')
+    BadReport = Field(Any, 'Bad Report')
+    State = Field(Any, 'State')
+    IsOffsetCorrection = Field(Any, 'Is Offset Correction')
+    IsOffsetFitted = Field(Any, 'Is Offset Fitted')
+    IsBGCorrection = Field(Any, 'Is BG Correction')
+    IsT0ShiftFitted = Field(Any, 'Is T0 Shift Fitted')
+    IsUseGDAS = Field(Any, 'Is Use GDAS')
+    IsUpwardMoving = Field(Any, 'Is Upward Moving')
+    OverShoot = Field(Any, 'Over Shoot')
+    UnderShoot = Field(Any, 'Under Shoot')
+    BGSamples = Field(Any, 'BG Samples')
+    Transmission3km = Field(Any, 'Transmission at 3 km')
+    Transmission6km = Field(Any, 'Transmission at 6 km')
+    Transmission9km = Field(Any, 'Transmission at 9 km')
+    Transmission12km = Field(Any, 'Transmission at 12 km')
+    Zenith = Field(Any, 'Zenith angle', unit=u.deg)
+    Azimuth = Field(Any, 'Azimuth angle', unit=u.deg)
+    FullOverlap = Field(Any, 'Full Overlap')
+    EndGroundLayer = Field(Any, 'End Ground Layer')
+    GroundLayerTrans = Field(Any, 'Ground Layer Trans')
+    Calimaness = Field(Any, 'Calimaness')
+    CloudLayerAlt = Field(Any, 'Altitude of cloud layer')
+    CloudLayerDens = Field(Any, 'Density of cloud layer')
+    Klett_k = Field(Any, 'Klett k')
+    PheCounts = Field(List[int], 'Phe Counts')
+    Offset = Field(Any, 'Offset')
+    Offset_Calculated = Field(Any, 'Offset calculated')
+    Offset_Fitted = Field(Any, 'Offset fitted')
+    Offset2 = Field(Any, 'Offset 2')
+    Offset3 = Field(Any, 'Offset 3')
+    Background1 = Field(Any, 'Background 1')
+    Background2 = Field(Any, 'Background 2')
+    BackgroundErr1 = Field(Any, 'Background error 1')
+    BackgroundErr2 = Field(Any, 'Background error 2')
+    RangeMax = Field(Any, 'Range max')
+    RangeMax_Clouds = Field(Any, 'Range max clouds')
+    ErrorCode = Field(Any, 'Error code')
+    ScaleHeight_fit = Field(Any, 'Scale Height fit')
+    Alpha_fit = Field(Any, 'Alpha fit')
+    Chi2Alpha_fit = Field(Any, 'Chi2 Alpha fit')
+    Alpha_firstY = Field(Any, 'Alpha first Y')
+    Alpha_Junge = Field(Any, 'Alpha Junge')
+    PBLHeight = Field(Any, 'PBL Height')
+    Chi2Full_fit = Field(Any, 'Chi2 Full fit')
+    SignalSamples = Field(Any, 'Signal Samples')
+    HWSwitch = Field(Any, 'HW Switch')
+    HWSwitchMaxOffset = Field(Any, 'HW Switch Max Offset')
+    NCollapse = Field(Any, 'N Collapse')
+    Shots = Field(Any, 'Shots')
+    T0Shift = Field(Any, 'T0 Shift')
+    Interval_0 = Field(Any, 'Interval 0')
+    RCS_min_perfect = Field(Any, 'RCS min perfect')
+    RCS_min_clouds = Field(Any, 'RCS min cloud')
+    RCS_min_mol = Field(Any, 'RCS min mol')
+    LIDAR_ratio = Field(Any, 'LIDAR ratio')
+    LIDAR_ratio_Cloud = Field(Any, 'LIDAR ratio cloud')
+    LIDAR_ratio_Junge = Field(Any, 'LIDAR ratio Junge')
 
 def load_camera_geometry():
     """Load camera geometry from bundled resources of this repo"""
@@ -223,6 +286,9 @@ def __init__(self, input_url=None, config=None, parent=None, **kwargs):
 
         self.metadata = self.parse_metadata_info()
 
+        if not self.is_simulation:
+            self.laser = self.parse_laser_info()
+
         # Retrieving the data level (so far HARDCODED Sorcerer)
         self.datalevel = DataLevel.DL0
 
@@ -381,13 +447,13 @@ def get_run_info_from_name(file_name):
         elif re.match(mask_data_superstar, file_name) is not None:
             parsed_info = re.match(mask_data_superstar, file_name)
             telescope = None
-            run_number = int(parsed_info.grou(1))
+            run_number = int(parsed_info.group(1))
             datalevel = MARSDataLevel.SUPERSTAR
             is_mc = False
         elif re.match(mask_data_melibea, file_name) is not None:
             parsed_info = re.match(mask_data_melibea, file_name)
             telescope = None
-            run_number = int(parsed_info.grou(1))
+            run_number = int(parsed_info.group(1))
             datalevel = MARSDataLevel.MELIBEA
             is_mc = False
         elif re.match(mask_mc_calibrated, file_name) is not None:
@@ -1035,6 +1101,150 @@ def parse_metadata_info(self):
 
         return metadata
 
+    def parse_laser_info(self):
+        laser_info_array_list_runh = [
+            'MReportLaser.MReport.fUniqueID',
+            'MReportLaser.MReport.fBits',
+            'MTimeLaser.fMjd',
+            'MTimeLaser.fTime.fMilliSec',
+            'MReportLaser.MReport.fBadReport',
+            'MReportLaser.MReport.fState',
+            'MReportLaser.fIsOffsetCorrection',
+            'MReportLaser.fIsOffsetFitted',
+            'MReportLaser.fIsBGCorrection',
+            'MReportLaser.fIsT0ShiftFitted',
+            'MReportLaser.fIsUseGDAS',
+            'MReportLaser.fIsUpwardMoving',
+            'MReportLaser.fOverShoot',
+            'MReportLaser.fUnderShoot',
+            'MReportLaser.fBGSamples',
+            'MReportLaser.fTransmission3km',
+            'MReportLaser.fTransmission6km',
+            'MReportLaser.fTransmission9km',
+            'MReportLaser.fTransmission12km',
+            'MReportLaser.fZenith',
+            'MReportLaser.fAzimuth',
+            'MReportLaser.fFullOverlap',
+            'MReportLaser.fEndGroundLayer',
+            'MReportLaser.fGroundLayerTrans',
+            'MReportLaser.fKlett_k',
+            'MReportLaser.fPheCounts',
+            'MReportLaser.fCalimaness',
+            'MReportLaser.fCloudLayerAlt',
+            'MReportLaser.fCloudLayerDens',
+            'MReportLaser.fOffset',
+            'MReportLaser.fOffset_Calculated',
+            'MReportLaser.fOffset_Fitted',
+            'MReportLaser.fOffset2',
+            'MReportLaser.fOffset3',
+            'MReportLaser.fBackground1',
+            'MReportLaser.fBackground2',
+            'MReportLaser.fBackgroundErr1',
+            'MReportLaser.fBackgroundErr2',
+            'MReportLaser.fRangeMax',
+            'MReportLaser.fRangeMax_Clouds',
+            'MReportLaser.fErrorCode',
+            'MReportLaser.fScaleHeight_fit',
+            'MReportLaser.fChi2Alpha_fit',
+            'MReportLaser.fChi2Alpha_fit',
+            'MReportLaser.fAlpha_firstY',
+            'MReportLaser.fAlpha_Junge',
+            'MReportLaser.fPBLHeight',
+            'MReportLaser.fChi2Full_fit',
+            'MReportLaser.fHWSwitchMaxOffset',
+            'MReportLaser.fNCollapse',
+            'MReportLaser.fShots',
+            'MReportLaser.fT0Shift',
+            'MReportLaser.fInterval_0',
+            'MReportLaser.fRCS_min_perfect',
+            'MReportLaser.fRCS_min_clouds',
+            'MReportLaser.fRCS_min_mol',
+            'MReportLaser.fLIDAR_ratio',
+            'MReportLaser.fLIDAR_ratio_Cloud',
+            'MReportLaser.fLIDAR_ratio_Junge',
+        ]
+
+        old_mjd = None
+        unique_reports = []
+        for rootf in self.files_:
+            try:
+                laser_info_runh = rootf['Laser'].arrays(
+                laser_info_array_list_runh, library="np")
+                mjd_value = laser_info_runh['MTimeLaser.fMjd']
+                millisec_value = laser_info_runh['MTimeLaser.fTime.fMilliSec']
+            except KeyError as e:
+                print(f"Required key not found in the file {rootf}: {e}")
+                continue
+            millisec_seconds = millisec_value * msec2sec
+            combined_mjd_value = mjd_value + millisec_seconds / 86400
+            for index in range(len(mjd_value)):
+                if combined_mjd_value[index] != old_mjd:
+                    laser = ReportLaserContainer()
+                    laser.MJD = combined_mjd_value[index]
+                    laser.UniqueID = laser_info_runh['MReportLaser.MReport.fUniqueID'][index]
+                    laser.Bits = laser_info_runh['MReportLaser.MReport.fBits'][index]
+                    laser.BadReport = laser_info_runh['MReportLaser.MReport.fBadReport'][index]
+                    laser.State = laser_info_runh['MReportLaser.MReport.fState'][index]
+                    laser.IsOffsetCorrection = laser_info_runh['MReportLaser.fIsOffsetCorrection'][index]
+                    laser.IsOffsetFitted = laser_info_runh['MReportLaser.fIsOffsetFitted'][index]
+                    laser.IsBGCorrection = laser_info_runh['MReportLaser.fIsBGCorrection'][index]
+                    laser.IsT0ShiftFitted = laser_info_runh['MReportLaser.fIsT0ShiftFitted'][index]
+                    laser.IsUseGDAS = laser_info_runh['MReportLaser.fIsUseGDAS'][index]
+                    laser.IsUpwardMoving = laser_info_runh['MReportLaser.fIsUpwardMoving'][index]
+                    laser.OverShoot = laser_info_runh['MReportLaser.fOverShoot'][index]
+                    laser.UnderShoot = laser_info_runh['MReportLaser.fUnderShoot'][index]
+                    laser.BGSamples = laser_info_runh['MReportLaser.fBGSamples'][index]
+                    laser.Transmission3km = laser_info_runh['MReportLaser.fTransmission3km'][index]
+                    laser.Transmission6km = laser_info_runh['MReportLaser.fTransmission6km'][index]
+                    laser.Transmission9km = laser_info_runh['MReportLaser.fTransmission9km'][index]
+                    laser.Transmission12km = laser_info_runh['MReportLaser.fTransmission12km'][index]
+                    laser.Zenith = laser_info_runh['MReportLaser.fZenith'][index] * u.deg
+                    laser.Azimuth = laser_info_runh['MReportLaser.fAzimuth'][index] * u.deg
+                    laser.FullOverlap = laser_info_runh['MReportLaser.fFullOverlap'][index]
+                    laser.EndGroundLayer = laser_info_runh['MReportLaser.fEndGroundLayer'][index]
+                    laser.GroundLayerTrans = laser_info_runh['MReportLaser.fGroundLayerTrans'][index]
+                    laser.Klett_k = laser_info_runh['MReportLaser.fKlett_k'][index]
+                    laser.PheCounts = laser_info_runh['MReportLaser.fPheCounts'][index]
+                    laser.Calimaness = laser_info_runh['MReportLaser.fCalimaness'][index]
+                    laser.CloudLayerAlt = laser_info_runh['MReportLaser.fCloudLayerAlt'][index]
+                    laser.CloudLayerDens = laser_info_runh['MReportLaser.fCloudLayerDens'][index]
+                    laser.Offset = laser_info_runh['MReportLaser.fOffset'][index]
+                    laser.Offset_Calculated = laser_info_runh['MReportLaser.fOffset_Calculated'][index]
+                    laser.Offset_Fitted = laser_info_runh['MReportLaser.fOffset_Fitted'][index]
+                    laser.Offset2 = laser_info_runh['MReportLaser.fOffset2'][index]
+                    laser.Offset3 = laser_info_runh['MReportLaser.fOffset3'][index]
+                    laser.Background1 = laser_info_runh['MReportLaser.fBackground1'][index]
+                    laser.Background2 = laser_info_runh['MReportLaser.fBackground2'][index]
+                    laser.BackgroundErr1 = laser_info_runh['MReportLaser.fBackgroundErr1'][index]
+                    laser.BackgroundErr2 = laser_info_runh['MReportLaser.fBackgroundErr2'][index]
+                    laser.RangeMax = laser_info_runh['MReportLaser.fRangeMax'][index]
+                    laser.RangeMax_Clouds = laser_info_runh['MReportLaser.fRangeMax_Clouds'][index]
+                    laser.ErrorCode = laser_info_runh['MReportLaser.fErrorCode'][index]
+                    laser.ScaleHeight_fit = laser_info_runh['MReportLaser.fScaleHeight_fit'][index]
+                    laser.Alpha_fit = laser_info_runh['MReportLaser.fChi2Alpha_fit'][index]
+                    laser.Chi2Alpha_fit = laser_info_runh['MReportLaser.fChi2Alpha_fit'][index]
+                    laser.Alpha_firstY = laser_info_runh['MReportLaser.fAlpha_firstY'][index]
+                    laser.Alpha_Junge = laser_info_runh['MReportLaser.fAlpha_Junge'][index]
+                    laser.PBLHeight = laser_info_runh['MReportLaser.fPBLHeight'][index]
+                    laser.Chi2Full_fit = laser_info_runh['MReportLaser.fChi2Full_fit'][index]
+                    laser.HWSwitchMaxOffset = laser_info_runh['MReportLaser.fHWSwitchMaxOffset'][index]
+                    laser.NCollapse = laser_info_runh['MReportLaser.fNCollapse'][index]
+                    laser.Shots = laser_info_runh['MReportLaser.fShots'][index]
+                    laser.T0Shift = laser_info_runh['MReportLaser.fT0Shift'][index]
+                    laser.Interval_0 = laser_info_runh['MReportLaser.fInterval_0'][index]
+                    laser.RCS_min_perfect = laser_info_runh['MReportLaser.fRCS_min_perfect'][index]
+                    laser.RCS_min_clouds = laser_info_runh['MReportLaser.fRCS_min_clouds'][index]
+                    laser.RCS_min_mol = laser_info_runh['MReportLaser.fRCS_min_mol'][index]
+                    laser.LIDAR_ratio = laser_info_runh['MReportLaser.fLIDAR_ratio'][index]
+                    laser.LIDAR_ratio_Cloud = laser_info_runh['MReportLaser.fLIDAR_ratio_Cloud'][index]
+                    laser.LIDAR_ratio_Junge = laser_info_runh['MReportLaser.fLIDAR_ratio_Junge'][index]
+
+                    unique_reports.append(laser)
+                    print(mjd_value)
+                old_mjd = combined_mjd_value[index]
+
+        return unique_reports
+
     def parse_simulation_header(self):
         """
         Parse the simulation information from the RunHeaders tree.

ctapipe_io_magic/tests/test_lidar.py

@@ -0,0 +1,29 @@
+import os
+from pathlib import Path
+
+import pytest
+
+test_data = Path(os.getenv("MAGIC_TEST_DATA", "test_data")).absolute()
+test_calibrated_real_dir = test_data / "real/calibrated"
+test_calibrated_real = test_calibrated_real_dir / "20210314_M1_05095172.002_Y_CrabNebula-W0.40+035.root"
+
+@pytest.fixture
+def sample_report_laser():
+    from ctapipe_io_magic import MAGICEventSource
+    with MAGICEventSource(input_url=test_calibrated_real,process_run=True) as source:
+        laser = source.laser
+        return laser
+
+def test_lidar_parameters(sample_report_laser):
+    '''Comparing the read Lidar report parameters with hardcoded values'''
+    laser = sample_report_laser
+    assert len(laser) == 1
+    for report in laser:
+        phe_counts = [425, 764, 683, 711, 746, 721, 658, 731, 721, 677, 707, 754, 684, 728, 751, 711, 690, 695, 738, 687, 706, 650, 719, 724, 670, 710, 737, 742, 708, 705, 759, 179623, 2464563, 1758611, 1280821, 986393, 812893, 722703, 743947, 1691760, 3368201, 4376413, 4746189, 4724701, 4540566, 4299863, 4011618, 3763189, 3499434, 3253066, 3037985, 2831435, 2643358, 2474402, 2323471, 2178678, 2045924, 1925799, 1810810, 1706918, 1607277, 1515936, 1436227, 1356119, 1280398, 1218400, 1147896, 1086930, 1033397, 979359, 930597, 887483, 843525, 799898, 765021, 726615, 696215, 661314, 634164, 607068, 581747, 558592, 535797, 516721, 495480, 476815, 459176, 442030, 424442, 410585, 397295, 383475, 371108, 357603, 347352, 335689, 326221, 314597, 303917, 295578, 288181, 277774, 270177, 261659, 254427, 248027, 238576, 231383, 225801, 220153, 213184, 208052, 201773, 195992, 191629, 187563, 181116, 177353, 171939, 168135, 165829, 159561, 157485, 152367, 149774, 145431, 142500, 138534, 267828, 256294, 245388, 233980, 225530, 214241, 205567, 197794, 189473, 182911, 175534, 167783, 162659, 156169, 151238, 144940, 139989, 136526, 130620, 126002, 122111, 118354, 114672, 111281, 106334, 103802, 99964, 97115, 93948, 91785, 88386, 86373, 84068, 81654, 78832, 77558, 75214, 72567, 70234, 68192, 66504, 65121, 62752, 61142, 60156, 57868, 56230, 54986, 53344, 53607, 52631, 49575, 49334, 47419, 46824, 45539, 43945, 42663, 42277, 42386, 40259, 39347, 38583, 38002, 37322, 36039, 35217, 34848, 34428, 33457, 32524, 31969, 31401, 30554, 30320, 29564, 28803, 28435, 27901, 27275, 26462, 26659, 26215, 25136, 24792, 24609, 23965, 23048, 23350, 23042, 22148, 21619, 21450, 21027, 21040, 20045, 19757, 19730, 19343, 19072, 18729, 19274, 17827, 17871, 17421, 17081, 16846, 16688, 16752, 16059, 15695, 15649, 15363, 15133, 15618, 15505, 15868, 18740, 20848, 19722, 18449, 15898, 14839, 13871, 13201, 13033, 13031, 12772, 23447, 22891, 21864, 21490, 20727, 20305, 19659, 18968, 18497, 18345, 17716, 17234, 16854, 16291, 15821, 15401, 15487, 15011, 14650, 14303, 13954, 13555, 13169, 13125, 12681, 12523, 12188, 11897, 11597, 11074, 11024, 10835, 10722, 10498, 10181, 10086, 9662, 9426, 9273, 9142, 8951, 8862, 8628, 8530, 8224, 8264, 8024, 7831, 7822, 7726, 7589, 7503, 7280, 7255, 7056, 6966, 6956, 6777, 6780, 6572, 6482, 6485, 6429, 6262, 6201, 6088, 6046, 5915, 5924, 5769, 5687, 5625, 5524, 5506, 5426, 5379, 5379, 5212, 5215, 5132, 5266, 5106, 5101, 5003, 4974, 4811, 4942, 4740, 4822, 4714, 4626, 4698, 4632, 4568, 4574, 4336, 4467, 4400, 4299, 4255, 4295, 4255, 4122, 4339, 4209, 4181, 4100, 4080, 4115, 4121, 4102, 4002, 3896, 3971, 3894, 3909, 3898, 3835, 3782, 3776, 3699, 3779, 3660, 3696, 3650, 3667, 3685, 3681, 7131, 7254, 7312, 7134, 6978, 7045, 6834, 6834, 6788, 6687, 6639, 6634, 6511, 6787, 6603, 6464, 6388, 6437, 6419, 6323, 6360, 6418, 6294, 6370, 6247, 6239, 6241, 6178, 6056, 6059, 5981, 6160, 6118, 5884, 5885, 5893, 5958, 5988, 5869, 6088, 5947, 5961, 5964, 5878, 5866, 5879, 5751, 5918, 5819, 5843, 5750, 5869, 5920, 5854, 5863, 5797, 5856, 5688, 5788, 5754, 5690, 5799, 5664, 5863, 5601, 5660, 5693, 5790, 5886, 5787, 5792, 5792, 5665, 5651, 5771, 5831, 5640, 5804, 5748, 5696, 5758, 5704, 5567, 5712, 5693, 5839, 5748, 5725, 5830, 5675, 5586, 5837, 5627, 5727, 5710, 5603, 5709, 5637, 5752, 5886, 5757, 5639, 5726, 5677, 5650, 5696, 5654, 5832, 5692, 5530, 5694, 5706, 5675, 5608, 5683, 5595, 5657, 5745, 5596, 5759, 5737, 5731, 5630, 5671, 5683, 5632, 5695, 5705]
+        assert report.Transmission12km == pytest.approx(0.89, abs=8.9e-07)
+        assert report.IsUseGDAS == [False]
+        azimuth_degrees = report.Azimuth.value
+        assert azimuth_degrees == pytest.approx(276.63)
+        assert list(report.PheCounts) == phe_counts
+        assert report.Bits == 50331656
+        assert report.MJD == pytest.approx(59286.91349633, rel = 1e-10, abs = 1e-6)