From 3bc1bf2cfcd186bee73cb5ff1ebf01c0b9735145 Mon Sep 17 00:00:00 2001 From: nzywucka Date: Thu, 23 Nov 2023 15:06:49 +0100 Subject: [PATCH 01/38] Reading of uncalibrated reports from Magic lidar and implementation of container to keep the parameters. --- ctapipe_io_magic/__init__.py | 87 +++++++++++++++++++++++++++++++++++- 1 file changed, 86 insertions(+), 1 deletion(-) diff --git a/ctapipe_io_magic/__init__.py b/ctapipe_io_magic/__init__.py index 846c928..b2d69f0 100644 --- a/ctapipe_io_magic/__init__.py +++ b/ctapipe_io_magic/__init__.py @@ -10,14 +10,16 @@ import uproot import logging import numpy as np +from numpy import nan from pathlib import Path from decimal import Decimal from astropy import units as u from astropy.time import Time +from astropy.coordinates import Angle from pkg_resources import resource_filename from ctapipe.io import EventSource, DataLevel -from ctapipe.core import Provenance +from ctapipe.core import Provenance, Container, Field, Map from ctapipe.core.traits import Bool, UseEnum from ctapipe.coordinates import CameraFrame @@ -80,6 +82,29 @@ DATA_MAGIC_LST_TRIGGER: EventType.SUBARRAY, } +class ReportLaserContainer(Container): + """ Container for Magic laser parameters """ + + OverShoot = Field(nan, 'Over Shoot') + UnderShoot = Field(nan, 'Under Shoot') + Transmission3km = Field(np.float32(np.nan), 'Transmission at 3 km') + Transmission6km = Field(np.float32(np.nan), 'Transmission at 6 km') + Transmission9km = Field(np.float32(np.nan), 'Transmission at 9 km') + Transmission12km = Field(np.float32(np.nan), 'Transmission at 12 km') + Zenith = Field(Angle(np.nan, u.deg), 'Zenith angle', unit=u.deg) + Azimuth = Field(Angle(np.nan, u.deg), 'Azimuth angle', unit=u.deg) + CloudFWHM = Field(np.float32(np.nan), 'Cloud FWHM') + CloudBase = Field(np.float32(np.nan), 'Cloud Base') + CloudTop = Field(np.float32(np.nan), 'Cloud Top') + CloudTrans = Field(np.float32(np.nan), 'Cloud Trans') + CloudHM = Field(np.float32(np.nan), 'Cloud HM') + CloudHStd = Field(np.float32(np.nan), 'Cloud HStd') + CloudLR = Field(np.float32(np.nan), 'Cloud LR') + FullOverlap = Field(np.float32(np.nan), 'Full Overlap') + EndGroundLayer = Field(np.float32(np.nan), 'End Ground Layer') + GroundLayerTrans = Field(np.float32(np.nan), 'Ground Layer Trans') + Klett_k = Field(np.float32(np.nan), 'Klett k') + PheCounts = Field(np.float32(np.nan), 'Phe Counts') def load_camera_geometry(): """Load camera geometry from bundled resources of this repo""" @@ -205,6 +230,8 @@ def __init__(self, input_url=None, config=None, parent=None, **kwargs): self.mars_datalevel = run_info[3][0] self.metadata = self.parse_metadata_info() + + self.laser = self.parse_laser_info() # Retrieving the data level (so far HARDCODED Sorcerer) self.datalevel = DataLevel.DL0 @@ -888,6 +915,64 @@ def parse_metadata_info(self): ) return metadata + + def parse_laser_info(self): + + laser_info_array_list_runh = [ + 'MReportLaser.fOverShoot', + 'MReportLaser.fUnderShoot', + 'MReportLaser.fTransmission3km', + 'MReportLaser.fTransmission6km', + 'MReportLaser.fTransmission9km', + 'MReportLaser.fTransmission12km', + 'MReportLaser.fZenith', + 'MReportLaser.fAzimuth', + 'MReportLaser.fCloudFWHM[10]', + 'MReportLaser.fCloudBase[10]', + 'MReportLaser.fCloudTop[10]', + 'MReportLaser.fCloudTrans[10]', + 'MReportLaser.fCloudHM[10]', + 'MReportLaser.fCloudHStd[10]', + 'MReportLaser.fCloudLR[10]', + 'MReportLaser.fFullOverlap', + 'MReportLaser.fEndGroundLayer', + 'MReportLaser.fGroundLayerTrans', + 'MReportLaser.fKlett_k', + 'MReportLaser.fPheCounts', + ] + + laser = ReportLaserContainer() # Create an instance of ReportLaserContainer + + for rootf in self.files_: + laser_info_runh = rootf['Laser'].arrays( + laser_info_array_list_runh, library="np" + ) + + # Populate the attributes of the laser object with values from report_laser + laser.OverShoot = int(laser_info_runh['MReportLaser.fOverShoot'][0]) + laser.UnderShoot = int(laser_info_runh['MReportLaser.fUnderShoot'][0]) + laser.Transmission3km = laser_info_runh['MReportLaser.fTransmission3km'][0] + laser.Transmission6km = laser_info_runh['MReportLaser.fTransmission6km'][0] + laser.Transmission9km = laser_info_runh['MReportLaser.fTransmission9km'][0] + laser.Transmission12km = laser_info_runh['MReportLaser.fTransmission12km'][0] + laser.Zenith = laser_info_runh['MReportLaser.fZenith'][0] + laser.Azimuth = laser_info_runh['MReportLaser.fAzimuth'][0] + laser.CloudFWHM = laser_info_runh['MReportLaser.fCloudFWHM[10]'][0] + laser.CloudBase = laser_info_runh['MReportLaser.fCloudBase[10]'][0] + laser.CloudTop = laser_info_runh['MReportLaser.fCloudTop[10]'][0] + laser.CloudTrans = laser_info_runh['MReportLaser.fCloudTrans[10]'][0] + laser.CloudHM = laser_info_runh['MReportLaser.fCloudHM[10]'][0] + laser.CloudHStd = laser_info_runh['MReportLaser.fCloudHStd[10]'][0] + laser.CloudLR = laser_info_runh['MReportLaser.fCloudLR[10]'][0] + laser.FullOverlap = laser_info_runh['MReportLaser.fFullOverlap'][0] + laser.EndGroundLayer = laser_info_runh['MReportLaser.fEndGroundLayer'][0] + laser.GroundLayerTrans = laser_info_runh['MReportLaser.fGroundLayerTrans'][0] + laser.Klett_k = laser_info_runh['MReportLaser.fKlett_k'][0] + laser.PheCounts = laser_info_runh['MReportLaser.fPheCounts'][0] + + print("laser init = ", laser) + return laser + def parse_simulation_header(self): """ From fbd80b1f081328f735f92c827d32425b1169c83f Mon Sep 17 00:00:00 2001 From: nzywucka Date: Thu, 23 Nov 2023 15:18:37 +0100 Subject: [PATCH 02/38] Removed print. --- ctapipe_io_magic/__init__.py | 1 - 1 file changed, 1 deletion(-) diff --git a/ctapipe_io_magic/__init__.py b/ctapipe_io_magic/__init__.py index b2d69f0..90bd619 100644 --- a/ctapipe_io_magic/__init__.py +++ b/ctapipe_io_magic/__init__.py @@ -970,7 +970,6 @@ def parse_laser_info(self): laser.Klett_k = laser_info_runh['MReportLaser.fKlett_k'][0] laser.PheCounts = laser_info_runh['MReportLaser.fPheCounts'][0] - print("laser init = ", laser) return laser From 7a975652dd4bf7d4bdc6bf0009a3cbc41d16ef00 Mon Sep 17 00:00:00 2001 From: nzywucka Date: Wed, 29 Nov 2023 21:53:03 +0100 Subject: [PATCH 03/38] In __init__.py, more parameters were added to ReportLaserContainer() and to parse_laser_info(). --- ctapipe_io_magic/__init__.py | 155 +++++++++++++++++++++++++++++++++-- 1 file changed, 149 insertions(+), 6 deletions(-) diff --git a/ctapipe_io_magic/__init__.py b/ctapipe_io_magic/__init__.py index 90bd619..60cc559 100644 --- a/ctapipe_io_magic/__init__.py +++ b/ctapipe_io_magic/__init__.py @@ -84,7 +84,20 @@ class ReportLaserContainer(Container): """ Container for Magic laser parameters """ - + UniqueID = Field(None, 'No.') + Bits = Field(None, 'ID') + MJD = Field(None, 'Modified Julian Date') + MiliSec = Field(None, 'Mili sec') + NanoSec = Field(None, 'Nano sec') + BadReport = Field(np.bool, 'Bad Report') + Telescope = Field(None, 'Telescope') + State = Field(None, 'State') + IsOffsetCorrection = Field(np.bool, 'Is Offset Correction') + IsOffsetFitted = Field(np.bool, 'Is Offset Fitted') + IsBGCorrection = Field(np.bool, 'Is BG Correction') + IsT0ShiftFitted = Field(np.bool, 'Is T0 Shift Fitted') + IsUseGDAS = Field(np.bool, 'Is Use GDAS') + IsUpwardMoving = Field(np.bool, 'Is Upward Moving') OverShoot = Field(nan, 'Over Shoot') UnderShoot = Field(nan, 'Under Shoot') Transmission3km = Field(np.float32(np.nan), 'Transmission at 3 km') @@ -103,8 +116,44 @@ class ReportLaserContainer(Container): FullOverlap = Field(np.float32(np.nan), 'Full Overlap') EndGroundLayer = Field(np.float32(np.nan), 'End Ground Layer') GroundLayerTrans = Field(np.float32(np.nan), 'Ground Layer Trans') + Calimaness = Field(np.float32(np.nan), 'Calimaness') + CloudLayerAlt = Field(np.float32(np.nan), 'Altitude of cloud layer') + CloudLayerDens = Field(np.float32(np.nan), 'Density of cloud layer') Klett_k = Field(np.float32(np.nan), 'Klett k') - PheCounts = Field(np.float32(np.nan), 'Phe Counts') + PheCounts = Field([], 'Phe Counts') + Offset = Field(np.float32(np.nan), 'Offset') + Offset_Calculated = Field(np.float32(np.nan), 'Offset calculated') + Offset_Fitted = Field(np.float32(np.nan), 'Offset fitted') + Offset2 = Field(np.float32(np.nan), 'Offset 2') + Offset3 = Field(np.float32(np.nan), 'Offset 3') + Background1 = Field(np.float32(np.nan), 'Background 1') + Background2 = Field(np.float32(np.nan), 'Background 2') + BackgroundErr1 = Field(np.float32(np.nan), 'Background error 1') + BackgroundErr2 = Field(np.float32(np.nan), 'Background error 2') + RangeMax = Field(np.float32(np.nan), 'Range max') + RangeMax_Clouds = Field(np.float32(np.nan), 'Range max clouds') + ErrorCode = Field(None, 'Error code') + ScaleHeight_fit = Field(np.float32(np.nan), 'Scale Height fit') + Alpha_fit = Field(np.float32(np.nan), 'Alpha fit') + Chi2Alpha_fit = Field(np.float32(np.nan), 'Chi2 Alpha fit') + Alpha_firstY = Field(np.float32(np.nan), 'Alpha first Y') + Alpha_Junge = Field(np.float32(np.nan), 'Alpha Junge') + PBLHeight = Field(np.float32(np.nan), 'PBL Height') + Chi2Full_fit = Field(np.float32(np.nan), 'Chi2 Full fit') + BGSamples = Field(np.float32(np.nan), 'BG Samples') + SignalSamples = Field(np.float32(np.nan), 'Signal Samples') + HWSwitch = Field(np.float32(np.nan), 'HW Switch') + HWSwitchMaxOffset = Field(np.float32(np.nan), 'HW Switch Max Offset') + NCollapse = Field(np.float32(np.nan), 'N Collapse') + Shots = Field(np.float32(np.nan), 'Shots') + T0Shift = Field(np.float32(np.nan), 'T0 Shift') + Interval_0 = Field(np.float32(np.nan), 'Interval 0') + RCS_min_perfect = Field(np.float32(np.nan), 'RCS min perfect') + RCS_min_clouds = Field(np.float32(np.nan), 'RCS min cloud') + RCS_min_mol = Field(np.float32(np.nan), 'RCS min mol') + LIDAR_ratio = Field(np.float32(np.nan), 'LIDAR ratio') + LIDAR_ratio_Cloud = Field(np.float32(np.nan), 'LIDAR ratio cloud') + LIDAR_ratio_Junge = Field(np.float32(np.nan), 'LIDAR ratio Junge') def load_camera_geometry(): """Load camera geometry from bundled resources of this repo""" @@ -919,6 +968,21 @@ def parse_metadata_info(self): def parse_laser_info(self): laser_info_array_list_runh = [ + 'MReportLaser.MReport.fUniqueID', + 'MReportLaser.MReport.fBits', + 'MTimeLaser.fMjd', + #'MTimeLaser.fTime', + 'MTimeLaser.fTime.fMilliSec', + 'MTimeLaser.fNanoSec', + 'MReportLaser.MReport.fBadReport', + 'MReportLaser.MReport.fTelescope', + 'MReportLaser.MReport.fState', + 'MReportLaser.fIsOffsetCorrection', + 'MReportLaser.fIsOffsetFitted', + 'MReportLaser.fIsBGCorrection', + 'MReportLaser.fIsT0ShiftFitted', + 'MReportLaser.fIsUseGDAS', + 'MReportLaser.fIsUpwardMoving', 'MReportLaser.fOverShoot', 'MReportLaser.fUnderShoot', 'MReportLaser.fTransmission3km', @@ -939,6 +1003,39 @@ def parse_laser_info(self): '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', ] laser = ReportLaserContainer() # Create an instance of ReportLaserContainer @@ -949,14 +1046,28 @@ def parse_laser_info(self): ) # Populate the attributes of the laser object with values from report_laser + laser.UniqueID = laser_info_runh['MReportLaser.MReport.fUniqueID'][0] + laser.Bits = laser_info_runh['MReportLaser.MReport.fBits'][0] + laser.MJD = laser_info_runh['MTimeLaser.fMjd'][0] + laser.MiliSec = laser_info_runh['MTimeLaser.fTime.fMilliSec'][0] + laser.NanoSec = laser_info_runh['MTimeLaser.fNanoSec'][0] + laser.BadReport = bool(laser_info_runh['MReportLaser.MReport.fBadReport'][0]) + laser.Telescope = laser_info_runh['MReportLaser.MReport.fTelescope'][0] + laser.State = laser_info_runh['MReportLaser.MReport.fState'][0] + laser.IsOffsetCorrection = bool(laser_info_runh['MReportLaser.fIsOffsetCorrection'][0]) + laser.IsOffsetFitted = bool(laser_info_runh['MReportLaser.fIsOffsetFitted'][0]) + laser.IsBGCorrection = bool(laser_info_runh['MReportLaser.fIsBGCorrection'][0]) + laser.IsT0ShiftFitted = bool(laser_info_runh['MReportLaser.fIsT0ShiftFitted'][0]) + laser.IsUseGDAS = bool(laser_info_runh['MReportLaser.fIsUseGDAS'][0]) + laser.IsUpwardMoving = bool(laser_info_runh['MReportLaser.fIsUpwardMoving'][0]) laser.OverShoot = int(laser_info_runh['MReportLaser.fOverShoot'][0]) laser.UnderShoot = int(laser_info_runh['MReportLaser.fUnderShoot'][0]) laser.Transmission3km = laser_info_runh['MReportLaser.fTransmission3km'][0] laser.Transmission6km = laser_info_runh['MReportLaser.fTransmission6km'][0] laser.Transmission9km = laser_info_runh['MReportLaser.fTransmission9km'][0] laser.Transmission12km = laser_info_runh['MReportLaser.fTransmission12km'][0] - laser.Zenith = laser_info_runh['MReportLaser.fZenith'][0] - laser.Azimuth = laser_info_runh['MReportLaser.fAzimuth'][0] + laser.Zenith = laser_info_runh['MReportLaser.fZenith'][0]* u.deg + laser.Azimuth = laser_info_runh['MReportLaser.fAzimuth'][0]* u.deg laser.CloudFWHM = laser_info_runh['MReportLaser.fCloudFWHM[10]'][0] laser.CloudBase = laser_info_runh['MReportLaser.fCloudBase[10]'][0] laser.CloudTop = laser_info_runh['MReportLaser.fCloudTop[10]'][0] @@ -969,9 +1080,41 @@ def parse_laser_info(self): laser.GroundLayerTrans = laser_info_runh['MReportLaser.fGroundLayerTrans'][0] laser.Klett_k = laser_info_runh['MReportLaser.fKlett_k'][0] laser.PheCounts = laser_info_runh['MReportLaser.fPheCounts'][0] - + laser.Calimaness = laser_info_runh['MReportLaser.fCalimaness'][0] + laser.CloudLayerAlt = laser_info_runh['MReportLaser.fCloudLayerAlt'][0] + laser.CloudLayerDens = laser_info_runh['MReportLaser.fCloudLayerDens'][0] + laser.Offset = laser_info_runh['MReportLaser.fOffset'][0] + laser.Offset_Calculated = laser_info_runh['MReportLaser.fOffset_Calculated'][0] + laser.Offset_Fitted = laser_info_runh['MReportLaser.fOffset_Fitted'][0] + laser.Offset2 = laser_info_runh['MReportLaser.fOffset2'][0] + laser.Offset3 = laser_info_runh['MReportLaser.fOffset3'][0] + laser.Background1 = laser_info_runh['MReportLaser.fBackground1'][0] + laser.Background2 = laser_info_runh['MReportLaser.fBackground2'][0] + laser.BackgroundErr1 = laser_info_runh['MReportLaser.fBackgroundErr1'][0] + laser.BackgroundErr2 = laser_info_runh['MReportLaser.fBackgroundErr2'][0] + laser.RangeMax = laser_info_runh['MReportLaser.fRangeMax'][0] + laser.RangeMax_Clouds = laser_info_runh['MReportLaser.fRangeMax_Clouds'][0] + laser.ErrorCode = laser_info_runh['MReportLaser.fErrorCode'][0] + laser.ScaleHeight_fit = laser_info_runh['MReportLaser.fScaleHeight_fit'][0] + laser.Alpha_fit = laser_info_runh['MReportLaser.fChi2Alpha_fit'][0] + laser.Chi2Alpha_fit = laser_info_runh['MReportLaser.fChi2Alpha_fit'][0] + laser.Alpha_firstY = laser_info_runh['MReportLaser.fAlpha_firstY'][0] + laser.Alpha_Junge = laser_info_runh['MReportLaser.fAlpha_Junge'][0] + laser.PBLHeight = laser_info_runh['MReportLaser.fPBLHeight'][0] + laser.Chi2Full_fit = laser_info_runh['MReportLaser.fChi2Full_fit'][0] + laser.HWSwitchMaxOffset = laser_info_runh['MReportLaser.fHWSwitchMaxOffset'][0] + laser.NCollapse = laser_info_runh['MReportLaser.fNCollapse'][0] + laser.Shots = laser_info_runh['MReportLaser.fShots'][0] + laser.T0Shift = laser_info_runh['MReportLaser.fT0Shift'][0] + laser.Interval_0 = laser_info_runh['MReportLaser.fInterval_0'][0] + laser.RCS_min_perfect = laser_info_runh['MReportLaser.fRCS_min_perfect'][0] + laser.RCS_min_clouds = laser_info_runh['MReportLaser.fRCS_min_clouds'][0] + laser.RCS_min_mol = laser_info_runh['MReportLaser.fRCS_min_mol'][0] + laser.LIDAR_ratio = laser_info_runh['MReportLaser.fLIDAR_ratio'][0] + laser.LIDAR_ratio_Cloud = laser_info_runh['MReportLaser.fLIDAR_ratio_Cloud'][0] + laser.LIDAR_ratio_Junge = laser_info_runh['MReportLaser.fLIDAR_ratio_Junge'][0] + return laser - def parse_simulation_header(self): """ From 3398c848213fc069f0ad53390eb998847659419f Mon Sep 17 00:00:00 2001 From: nzywucka Date: Mon, 4 Dec 2023 11:22:01 +0100 Subject: [PATCH 04/38] __init__.py changed regarding Laser test #84 comments in ReportLaserContainer() class and parse_laser_info() function. --- ctapipe_io_magic/__init__.py | 188 +++++++++++++++++------------------ 1 file changed, 93 insertions(+), 95 deletions(-) diff --git a/ctapipe_io_magic/__init__.py b/ctapipe_io_magic/__init__.py index 60cc559..955c2e9 100644 --- a/ctapipe_io_magic/__init__.py +++ b/ctapipe_io_magic/__init__.py @@ -19,7 +19,7 @@ from pkg_resources import resource_filename from ctapipe.io import EventSource, DataLevel -from ctapipe.core import Provenance, Container, Field, Map +from ctapipe.core import Provenance, Container, Field from ctapipe.core.traits import Bool, UseEnum from ctapipe.coordinates import CameraFrame @@ -87,19 +87,16 @@ class ReportLaserContainer(Container): UniqueID = Field(None, 'No.') Bits = Field(None, 'ID') MJD = Field(None, 'Modified Julian Date') - MiliSec = Field(None, 'Mili sec') - NanoSec = Field(None, 'Nano sec') - BadReport = Field(np.bool, 'Bad Report') - Telescope = Field(None, 'Telescope') + BadReport = Field(bool, 'Bad Report') State = Field(None, 'State') - IsOffsetCorrection = Field(np.bool, 'Is Offset Correction') - IsOffsetFitted = Field(np.bool, 'Is Offset Fitted') - IsBGCorrection = Field(np.bool, 'Is BG Correction') - IsT0ShiftFitted = Field(np.bool, 'Is T0 Shift Fitted') - IsUseGDAS = Field(np.bool, 'Is Use GDAS') - IsUpwardMoving = Field(np.bool, 'Is Upward Moving') - OverShoot = Field(nan, 'Over Shoot') - UnderShoot = Field(nan, 'Under Shoot') + IsOffsetCorrection = Field(bool, 'Is Offset Correction') + IsOffsetFitted = Field(bool, 'Is Offset Fitted') + IsBGCorrection = Field(bool, 'Is BG Correction') + IsT0ShiftFitted = Field(bool, 'Is T0 Shift Fitted') + IsUseGDAS = Field(bool, 'Is Use GDAS') + IsUpwardMoving = Field(bool, 'Is Upward Moving') + OverShoot = Field(nan, 'Over Shoot') + UnderShoot = Field(nan, 'Under Shoot') Transmission3km = Field(np.float32(np.nan), 'Transmission at 3 km') Transmission6km = Field(np.float32(np.nan), 'Transmission at 6 km') Transmission9km = Field(np.float32(np.nan), 'Transmission at 9 km') @@ -140,11 +137,10 @@ class ReportLaserContainer(Container): Alpha_Junge = Field(np.float32(np.nan), 'Alpha Junge') PBLHeight = Field(np.float32(np.nan), 'PBL Height') Chi2Full_fit = Field(np.float32(np.nan), 'Chi2 Full fit') - BGSamples = Field(np.float32(np.nan), 'BG Samples') - SignalSamples = Field(np.float32(np.nan), 'Signal Samples') - HWSwitch = Field(np.float32(np.nan), 'HW Switch') - HWSwitchMaxOffset = Field(np.float32(np.nan), 'HW Switch Max Offset') - NCollapse = Field(np.float32(np.nan), 'N Collapse') + SignalSamples = Field(np.float32(np.nan), 'Signal Samples') + HWSwitch = Field(np.float32(np.nan), 'HW Switch') + HWSwitchMaxOffset = Field(np.float32(np.nan), 'HW Switch Max Offset') + NCollapse = Field(np.float32(np.nan), 'N Collapse') Shots = Field(np.float32(np.nan), 'Shots') T0Shift = Field(np.float32(np.nan), 'T0 Shift') Interval_0 = Field(np.float32(np.nan), 'Interval 0') @@ -153,7 +149,7 @@ class ReportLaserContainer(Container): RCS_min_mol = Field(np.float32(np.nan), 'RCS min mol') LIDAR_ratio = Field(np.float32(np.nan), 'LIDAR ratio') LIDAR_ratio_Cloud = Field(np.float32(np.nan), 'LIDAR ratio cloud') - LIDAR_ratio_Junge = Field(np.float32(np.nan), 'LIDAR ratio Junge') + LIDAR_ratio_Junge = Field(np.float32(np.nan), 'LIDAR ratio Junge') def load_camera_geometry(): """Load camera geometry from bundled resources of this repo""" @@ -971,11 +967,8 @@ def parse_laser_info(self): 'MReportLaser.MReport.fUniqueID', 'MReportLaser.MReport.fBits', 'MTimeLaser.fMjd', - #'MTimeLaser.fTime', 'MTimeLaser.fTime.fMilliSec', - 'MTimeLaser.fNanoSec', 'MReportLaser.MReport.fBadReport', - 'MReportLaser.MReport.fTelescope', 'MReportLaser.MReport.fState', 'MReportLaser.fIsOffsetCorrection', 'MReportLaser.fIsOffsetFitted', @@ -1038,82 +1031,87 @@ def parse_laser_info(self): 'MReportLaser.fLIDAR_ratio_Junge', ] - laser = ReportLaserContainer() # Create an instance of ReportLaserContainer + laser = ReportLaserContainer() for rootf in self.files_: - laser_info_runh = rootf['Laser'].arrays( - laser_info_array_list_runh, library="np" - ) - - # Populate the attributes of the laser object with values from report_laser - laser.UniqueID = laser_info_runh['MReportLaser.MReport.fUniqueID'][0] - laser.Bits = laser_info_runh['MReportLaser.MReport.fBits'][0] - laser.MJD = laser_info_runh['MTimeLaser.fMjd'][0] - laser.MiliSec = laser_info_runh['MTimeLaser.fTime.fMilliSec'][0] - laser.NanoSec = laser_info_runh['MTimeLaser.fNanoSec'][0] - laser.BadReport = bool(laser_info_runh['MReportLaser.MReport.fBadReport'][0]) - laser.Telescope = laser_info_runh['MReportLaser.MReport.fTelescope'][0] - laser.State = laser_info_runh['MReportLaser.MReport.fState'][0] - laser.IsOffsetCorrection = bool(laser_info_runh['MReportLaser.fIsOffsetCorrection'][0]) - laser.IsOffsetFitted = bool(laser_info_runh['MReportLaser.fIsOffsetFitted'][0]) - laser.IsBGCorrection = bool(laser_info_runh['MReportLaser.fIsBGCorrection'][0]) - laser.IsT0ShiftFitted = bool(laser_info_runh['MReportLaser.fIsT0ShiftFitted'][0]) - laser.IsUseGDAS = bool(laser_info_runh['MReportLaser.fIsUseGDAS'][0]) - laser.IsUpwardMoving = bool(laser_info_runh['MReportLaser.fIsUpwardMoving'][0]) - laser.OverShoot = int(laser_info_runh['MReportLaser.fOverShoot'][0]) - laser.UnderShoot = int(laser_info_runh['MReportLaser.fUnderShoot'][0]) - laser.Transmission3km = laser_info_runh['MReportLaser.fTransmission3km'][0] - laser.Transmission6km = laser_info_runh['MReportLaser.fTransmission6km'][0] - laser.Transmission9km = laser_info_runh['MReportLaser.fTransmission9km'][0] - laser.Transmission12km = laser_info_runh['MReportLaser.fTransmission12km'][0] - laser.Zenith = laser_info_runh['MReportLaser.fZenith'][0]* u.deg - laser.Azimuth = laser_info_runh['MReportLaser.fAzimuth'][0]* u.deg - laser.CloudFWHM = laser_info_runh['MReportLaser.fCloudFWHM[10]'][0] - laser.CloudBase = laser_info_runh['MReportLaser.fCloudBase[10]'][0] - laser.CloudTop = laser_info_runh['MReportLaser.fCloudTop[10]'][0] - laser.CloudTrans = laser_info_runh['MReportLaser.fCloudTrans[10]'][0] - laser.CloudHM = laser_info_runh['MReportLaser.fCloudHM[10]'][0] - laser.CloudHStd = laser_info_runh['MReportLaser.fCloudHStd[10]'][0] - laser.CloudLR = laser_info_runh['MReportLaser.fCloudLR[10]'][0] - laser.FullOverlap = laser_info_runh['MReportLaser.fFullOverlap'][0] - laser.EndGroundLayer = laser_info_runh['MReportLaser.fEndGroundLayer'][0] - laser.GroundLayerTrans = laser_info_runh['MReportLaser.fGroundLayerTrans'][0] - laser.Klett_k = laser_info_runh['MReportLaser.fKlett_k'][0] - laser.PheCounts = laser_info_runh['MReportLaser.fPheCounts'][0] - laser.Calimaness = laser_info_runh['MReportLaser.fCalimaness'][0] - laser.CloudLayerAlt = laser_info_runh['MReportLaser.fCloudLayerAlt'][0] - laser.CloudLayerDens = laser_info_runh['MReportLaser.fCloudLayerDens'][0] - laser.Offset = laser_info_runh['MReportLaser.fOffset'][0] - laser.Offset_Calculated = laser_info_runh['MReportLaser.fOffset_Calculated'][0] - laser.Offset_Fitted = laser_info_runh['MReportLaser.fOffset_Fitted'][0] - laser.Offset2 = laser_info_runh['MReportLaser.fOffset2'][0] - laser.Offset3 = laser_info_runh['MReportLaser.fOffset3'][0] - laser.Background1 = laser_info_runh['MReportLaser.fBackground1'][0] - laser.Background2 = laser_info_runh['MReportLaser.fBackground2'][0] - laser.BackgroundErr1 = laser_info_runh['MReportLaser.fBackgroundErr1'][0] - laser.BackgroundErr2 = laser_info_runh['MReportLaser.fBackgroundErr2'][0] - laser.RangeMax = laser_info_runh['MReportLaser.fRangeMax'][0] - laser.RangeMax_Clouds = laser_info_runh['MReportLaser.fRangeMax_Clouds'][0] - laser.ErrorCode = laser_info_runh['MReportLaser.fErrorCode'][0] - laser.ScaleHeight_fit = laser_info_runh['MReportLaser.fScaleHeight_fit'][0] - laser.Alpha_fit = laser_info_runh['MReportLaser.fChi2Alpha_fit'][0] - laser.Chi2Alpha_fit = laser_info_runh['MReportLaser.fChi2Alpha_fit'][0] - laser.Alpha_firstY = laser_info_runh['MReportLaser.fAlpha_firstY'][0] - laser.Alpha_Junge = laser_info_runh['MReportLaser.fAlpha_Junge'][0] - laser.PBLHeight = laser_info_runh['MReportLaser.fPBLHeight'][0] - laser.Chi2Full_fit = laser_info_runh['MReportLaser.fChi2Full_fit'][0] - laser.HWSwitchMaxOffset = laser_info_runh['MReportLaser.fHWSwitchMaxOffset'][0] - laser.NCollapse = laser_info_runh['MReportLaser.fNCollapse'][0] - laser.Shots = laser_info_runh['MReportLaser.fShots'][0] - laser.T0Shift = laser_info_runh['MReportLaser.fT0Shift'][0] - laser.Interval_0 = laser_info_runh['MReportLaser.fInterval_0'][0] - laser.RCS_min_perfect = laser_info_runh['MReportLaser.fRCS_min_perfect'][0] - laser.RCS_min_clouds = laser_info_runh['MReportLaser.fRCS_min_clouds'][0] - laser.RCS_min_mol = laser_info_runh['MReportLaser.fRCS_min_mol'][0] - laser.LIDAR_ratio = laser_info_runh['MReportLaser.fLIDAR_ratio'][0] - laser.LIDAR_ratio_Cloud = laser_info_runh['MReportLaser.fLIDAR_ratio_Cloud'][0] - laser.LIDAR_ratio_Junge = laser_info_runh['MReportLaser.fLIDAR_ratio_Junge'][0] + try: + laser_info_runh = rootf['Laser'].arrays( + laser_info_array_list_runh, library="np" + ) + laser.UniqueID = laser_info_runh['MReportLaser.MReport.fUniqueID'][0] + laser.Bits = laser_info_runh['MReportLaser.MReport.fBits'][0] + mjd_value = laser_info_runh['MTimeLaser.fMjd'][0] + millisec_value = laser_info_runh['MTimeLaser.fTime.fMilliSec'][0] + laser.BadReport = bool(laser_info_runh['MReportLaser.MReport.fBadReport'][0]) + laser.State = laser_info_runh['MReportLaser.MReport.fState'][0] + laser.IsOffsetCorrection = bool(laser_info_runh['MReportLaser.fIsOffsetCorrection'][0]) + laser.IsOffsetFitted = bool(laser_info_runh['MReportLaser.fIsOffsetFitted'][0]) + laser.IsBGCorrection = bool(laser_info_runh['MReportLaser.fIsBGCorrection'][0]) + laser.IsT0ShiftFitted = bool(laser_info_runh['MReportLaser.fIsT0ShiftFitted'][0]) + laser.IsUseGDAS = bool(laser_info_runh['MReportLaser.fIsUseGDAS'][0]) + laser.IsUpwardMoving = bool(laser_info_runh['MReportLaser.fIsUpwardMoving'][0]) + laser.OverShoot = int(laser_info_runh['MReportLaser.fOverShoot'][0]) + laser.UnderShoot = int(laser_info_runh['MReportLaser.fUnderShoot'][0]) + laser.Transmission3km = laser_info_runh['MReportLaser.fTransmission3km'][0] + laser.Transmission6km = laser_info_runh['MReportLaser.fTransmission6km'][0] + laser.Transmission9km = laser_info_runh['MReportLaser.fTransmission9km'][0] + laser.Transmission12km = laser_info_runh['MReportLaser.fTransmission12km'][0] + laser.Zenith = laser_info_runh['MReportLaser.fZenith'][0]* u.deg + laser.Azimuth = laser_info_runh['MReportLaser.fAzimuth'][0]* u.deg + laser.CloudFWHM = laser_info_runh['MReportLaser.fCloudFWHM[10]'][0] + laser.CloudBase = laser_info_runh['MReportLaser.fCloudBase[10]'][0] + laser.CloudTop = laser_info_runh['MReportLaser.fCloudTop[10]'][0] + laser.CloudTrans = laser_info_runh['MReportLaser.fCloudTrans[10]'][0] + laser.CloudHM = laser_info_runh['MReportLaser.fCloudHM[10]'][0] + laser.CloudHStd = laser_info_runh['MReportLaser.fCloudHStd[10]'][0] + laser.CloudLR = laser_info_runh['MReportLaser.fCloudLR[10]'][0] + laser.FullOverlap = laser_info_runh['MReportLaser.fFullOverlap'][0] + laser.EndGroundLayer = laser_info_runh['MReportLaser.fEndGroundLayer'][0] + laser.GroundLayerTrans = laser_info_runh['MReportLaser.fGroundLayerTrans'][0] + laser.Klett_k = laser_info_runh['MReportLaser.fKlett_k'][0] + laser.PheCounts = laser_info_runh['MReportLaser.fPheCounts'][0] + laser.Calimaness = laser_info_runh['MReportLaser.fCalimaness'][0] + laser.CloudLayerAlt = laser_info_runh['MReportLaser.fCloudLayerAlt'][0] + laser.CloudLayerDens = laser_info_runh['MReportLaser.fCloudLayerDens'][0] + laser.Offset = laser_info_runh['MReportLaser.fOffset'][0] + laser.Offset_Calculated = laser_info_runh['MReportLaser.fOffset_Calculated'][0] + laser.Offset_Fitted = laser_info_runh['MReportLaser.fOffset_Fitted'][0] + laser.Offset2 = laser_info_runh['MReportLaser.fOffset2'][0] + laser.Offset3 = laser_info_runh['MReportLaser.fOffset3'][0] + laser.Background1 = laser_info_runh['MReportLaser.fBackground1'][0] + laser.Background2 = laser_info_runh['MReportLaser.fBackground2'][0] + laser.BackgroundErr1 = laser_info_runh['MReportLaser.fBackgroundErr1'][0] + laser.BackgroundErr2 = laser_info_runh['MReportLaser.fBackgroundErr2'][0] + laser.RangeMax = laser_info_runh['MReportLaser.fRangeMax'][0] + laser.RangeMax_Clouds = laser_info_runh['MReportLaser.fRangeMax_Clouds'][0] + laser.ErrorCode = laser_info_runh['MReportLaser.fErrorCode'][0] + laser.ScaleHeight_fit = laser_info_runh['MReportLaser.fScaleHeight_fit'][0] + laser.Alpha_fit = laser_info_runh['MReportLaser.fChi2Alpha_fit'][0] + laser.Chi2Alpha_fit = laser_info_runh['MReportLaser.fChi2Alpha_fit'][0] + laser.Alpha_firstY = laser_info_runh['MReportLaser.fAlpha_firstY'][0] + laser.Alpha_Junge = laser_info_runh['MReportLaser.fAlpha_Junge'][0] + laser.PBLHeight = laser_info_runh['MReportLaser.fPBLHeight'][0] + laser.Chi2Full_fit = laser_info_runh['MReportLaser.fChi2Full_fit'][0] + laser.HWSwitchMaxOffset = laser_info_runh['MReportLaser.fHWSwitchMaxOffset'][0] + laser.NCollapse = laser_info_runh['MReportLaser.fNCollapse'][0] + laser.Shots = laser_info_runh['MReportLaser.fShots'][0] + laser.T0Shift = laser_info_runh['MReportLaser.fT0Shift'][0] + laser.Interval_0 = laser_info_runh['MReportLaser.fInterval_0'][0] + laser.RCS_min_perfect = laser_info_runh['MReportLaser.fRCS_min_perfect'][0] + laser.RCS_min_clouds = laser_info_runh['MReportLaser.fRCS_min_clouds'][0] + laser.RCS_min_mol = laser_info_runh['MReportLaser.fRCS_min_mol'][0] + laser.LIDAR_ratio = laser_info_runh['MReportLaser.fLIDAR_ratio'][0] + laser.LIDAR_ratio_Cloud = laser_info_runh['MReportLaser.fLIDAR_ratio_Cloud'][0] + laser.LIDAR_ratio_Junge = laser_info_runh['MReportLaser.fLIDAR_ratio_Junge'][0] + + millisec_seconds = millisec_value * 1e-3 + combined_mjd_value = mjd_value + millisec_seconds / 86400 + laser.MJD = combined_mjd_value + except KeyError as e: + print(f"Required key not found in the file {rootf}: {e}") + continue + return laser def parse_simulation_header(self): From 867aa7313daaf8cf88f9096e95bd9ab2c8f50f27 Mon Sep 17 00:00:00 2001 From: nzywucka Date: Mon, 26 Feb 2024 12:17:46 +0100 Subject: [PATCH 05/38] Unnecessary spaces are removed from the file. --- ctapipe_io_magic/__init__.py | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/ctapipe_io_magic/__init__.py b/ctapipe_io_magic/__init__.py index 955c2e9..70df3c9 100644 --- a/ctapipe_io_magic/__init__.py +++ b/ctapipe_io_magic/__init__.py @@ -960,9 +960,8 @@ def parse_metadata_info(self): ) return metadata - - def parse_laser_info(self): + def parse_laser_info(self): laser_info_array_list_runh = [ 'MReportLaser.MReport.fUniqueID', 'MReportLaser.MReport.fBits', @@ -975,7 +974,7 @@ def parse_laser_info(self): 'MReportLaser.fIsBGCorrection', 'MReportLaser.fIsT0ShiftFitted', 'MReportLaser.fIsUseGDAS', - 'MReportLaser.fIsUpwardMoving', + 'MReportLaser.fIsUpwardMoving', 'MReportLaser.fOverShoot', 'MReportLaser.fUnderShoot', 'MReportLaser.fTransmission3km', @@ -984,10 +983,10 @@ def parse_laser_info(self): 'MReportLaser.fTransmission12km', 'MReportLaser.fZenith', 'MReportLaser.fAzimuth', - 'MReportLaser.fCloudFWHM[10]', - 'MReportLaser.fCloudBase[10]', - 'MReportLaser.fCloudTop[10]', - 'MReportLaser.fCloudTrans[10]', + 'MReportLaser.fCloudFWHM[10]', + 'MReportLaser.fCloudBase[10]', + 'MReportLaser.fCloudTop[10]', + 'MReportLaser.fCloudTrans[10]', 'MReportLaser.fCloudHM[10]', 'MReportLaser.fCloudHStd[10]', 'MReportLaser.fCloudLR[10]', @@ -1020,7 +1019,7 @@ def parse_laser_info(self): 'MReportLaser.fChi2Full_fit', 'MReportLaser.fHWSwitchMaxOffset', 'MReportLaser.fNCollapse', - 'MReportLaser.fShots', + 'MReportLaser.fShots', 'MReportLaser.fT0Shift', 'MReportLaser.fInterval_0', 'MReportLaser.fRCS_min_perfect', @@ -1103,7 +1102,7 @@ def parse_laser_info(self): laser.LIDAR_ratio = laser_info_runh['MReportLaser.fLIDAR_ratio'][0] laser.LIDAR_ratio_Cloud = laser_info_runh['MReportLaser.fLIDAR_ratio_Cloud'][0] laser.LIDAR_ratio_Junge = laser_info_runh['MReportLaser.fLIDAR_ratio_Junge'][0] - + millisec_seconds = millisec_value * 1e-3 combined_mjd_value = mjd_value + millisec_seconds / 86400 laser.MJD = combined_mjd_value @@ -1112,7 +1111,7 @@ def parse_laser_info(self): print(f"Required key not found in the file {rootf}: {e}") continue - return laser + return laser def parse_simulation_header(self): """ From aad168a1322b02369159384b6771a4c0dcee1368 Mon Sep 17 00:00:00 2001 From: nzywucka Date: Mon, 26 Feb 2024 13:48:27 +0100 Subject: [PATCH 06/38] Unnecessary spaces are removed --- ctapipe_io_magic/__init__.py | 1 - 1 file changed, 1 deletion(-) diff --git a/ctapipe_io_magic/__init__.py b/ctapipe_io_magic/__init__.py index 70df3c9..0ccb487 100644 --- a/ctapipe_io_magic/__init__.py +++ b/ctapipe_io_magic/__init__.py @@ -275,7 +275,6 @@ def __init__(self, input_url=None, config=None, parent=None, **kwargs): self.mars_datalevel = run_info[3][0] self.metadata = self.parse_metadata_info() - self.laser = self.parse_laser_info() # Retrieving the data level (so far HARDCODED Sorcerer) From 9cebd7983e1d5c28db04e51370aebb2351001a4c Mon Sep 17 00:00:00 2001 From: nzywucka Date: Mon, 26 Feb 2024 13:59:02 +0100 Subject: [PATCH 07/38] Unnecessary spaces removed --- ctapipe_io_magic/__init__.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ctapipe_io_magic/__init__.py b/ctapipe_io_magic/__init__.py index 0ccb487..abd67d5 100644 --- a/ctapipe_io_magic/__init__.py +++ b/ctapipe_io_magic/__init__.py @@ -1,4 +1,4 @@ -""" +=""" # Event source for MAGIC calibrated data files. # Requires uproot package (https://github.com/scikit-hep/uproot). """ @@ -1105,7 +1105,7 @@ def parse_laser_info(self): millisec_seconds = millisec_value * 1e-3 combined_mjd_value = mjd_value + millisec_seconds / 86400 laser.MJD = combined_mjd_value - + except KeyError as e: print(f"Required key not found in the file {rootf}: {e}") continue From 74d5c6230d7fd6eac53bb432421c51e2f1fba4b5 Mon Sep 17 00:00:00 2001 From: nzywucka Date: Mon, 26 Feb 2024 14:02:13 +0100 Subject: [PATCH 08/38] Typo removed --- ctapipe_io_magic/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ctapipe_io_magic/__init__.py b/ctapipe_io_magic/__init__.py index abd67d5..f5e917c 100644 --- a/ctapipe_io_magic/__init__.py +++ b/ctapipe_io_magic/__init__.py @@ -1,4 +1,4 @@ -=""" +""" # Event source for MAGIC calibrated data files. # Requires uproot package (https://github.com/scikit-hep/uproot). """ From 8bc648797bbeb827901eb0b45aa1cfae0632bf36 Mon Sep 17 00:00:00 2001 From: nzywucka Date: Mon, 26 Feb 2024 14:43:48 +0100 Subject: [PATCH 09/38] If statement for self.laser = self.parse_laser_info() added --- ctapipe_io_magic/__init__.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/ctapipe_io_magic/__init__.py b/ctapipe_io_magic/__init__.py index f5e917c..e904f4f 100644 --- a/ctapipe_io_magic/__init__.py +++ b/ctapipe_io_magic/__init__.py @@ -1,4 +1,4 @@ -""" +y""" # Event source for MAGIC calibrated data files. # Requires uproot package (https://github.com/scikit-hep/uproot). """ @@ -275,7 +275,9 @@ def __init__(self, input_url=None, config=None, parent=None, **kwargs): self.mars_datalevel = run_info[3][0] self.metadata = self.parse_metadata_info() - self.laser = self.parse_laser_info() + + if not self.is_simulation: + self.laser = self.parse_laser_info() # Retrieving the data level (so far HARDCODED Sorcerer) self.datalevel = DataLevel.DL0 From e0c1118107d02ec34e533f30c2e0ac7a005b80c6 Mon Sep 17 00:00:00 2001 From: nzywucka Date: Mon, 26 Feb 2024 14:48:15 +0100 Subject: [PATCH 10/38] Typo removed --- ctapipe_io_magic/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ctapipe_io_magic/__init__.py b/ctapipe_io_magic/__init__.py index e904f4f..df0d6ee 100644 --- a/ctapipe_io_magic/__init__.py +++ b/ctapipe_io_magic/__init__.py @@ -1,4 +1,4 @@ -y""" +""" # Event source for MAGIC calibrated data files. # Requires uproot package (https://github.com/scikit-hep/uproot). """ From 7e0252f8b3baccfc72103863f043d61b7786ab7a Mon Sep 17 00:00:00 2001 From: nzywucka Date: Mon, 26 Feb 2024 15:08:42 +0100 Subject: [PATCH 11/38] BGSamples parameter added --- ctapipe_io_magic/__init__.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/ctapipe_io_magic/__init__.py b/ctapipe_io_magic/__init__.py index df0d6ee..a0afaa9 100644 --- a/ctapipe_io_magic/__init__.py +++ b/ctapipe_io_magic/__init__.py @@ -97,6 +97,7 @@ class ReportLaserContainer(Container): IsUpwardMoving = Field(bool, 'Is Upward Moving') OverShoot = Field(nan, 'Over Shoot') UnderShoot = Field(nan, 'Under Shoot') + BGSamples = Field(np.float32(np.nan), 'BG Samples') Transmission3km = Field(np.float32(np.nan), 'Transmission at 3 km') Transmission6km = Field(np.float32(np.nan), 'Transmission at 6 km') Transmission9km = Field(np.float32(np.nan), 'Transmission at 9 km') @@ -978,6 +979,7 @@ def parse_laser_info(self): 'MReportLaser.fIsUpwardMoving', 'MReportLaser.fOverShoot', 'MReportLaser.fUnderShoot', + 'MReportLaser.fBGSamples', 'MReportLaser.fTransmission3km', 'MReportLaser.fTransmission6km', 'MReportLaser.fTransmission9km', @@ -1052,6 +1054,7 @@ def parse_laser_info(self): laser.IsUpwardMoving = bool(laser_info_runh['MReportLaser.fIsUpwardMoving'][0]) laser.OverShoot = int(laser_info_runh['MReportLaser.fOverShoot'][0]) laser.UnderShoot = int(laser_info_runh['MReportLaser.fUnderShoot'][0]) + laser.BGSamples = int(laser_info_runh['MReportLaser.fBGSamples'][0]) laser.Transmission3km = laser_info_runh['MReportLaser.fTransmission3km'][0] laser.Transmission6km = laser_info_runh['MReportLaser.fTransmission6km'][0] laser.Transmission9km = laser_info_runh['MReportLaser.fTransmission9km'][0] From be6a5366ab13eae06d7f0656d238f381159963f9 Mon Sep 17 00:00:00 2001 From: nzywucka Date: Tue, 12 Mar 2024 09:19:33 +0100 Subject: [PATCH 12/38] The [0]s have been removed from the arrays while reading the Lidar reports. Filtering for the unique reports has been added. --- ctapipe_io_magic/__init__.py | 148 ++++++++++++++++++----------------- 1 file changed, 78 insertions(+), 70 deletions(-) diff --git a/ctapipe_io_magic/__init__.py b/ctapipe_io_magic/__init__.py index a0afaa9..5787cb7 100644 --- a/ctapipe_io_magic/__init__.py +++ b/ctapipe_io_magic/__init__.py @@ -1040,76 +1040,84 @@ def parse_laser_info(self): laser_info_runh = rootf['Laser'].arrays( laser_info_array_list_runh, library="np" ) - laser.UniqueID = laser_info_runh['MReportLaser.MReport.fUniqueID'][0] - laser.Bits = laser_info_runh['MReportLaser.MReport.fBits'][0] - mjd_value = laser_info_runh['MTimeLaser.fMjd'][0] - millisec_value = laser_info_runh['MTimeLaser.fTime.fMilliSec'][0] - laser.BadReport = bool(laser_info_runh['MReportLaser.MReport.fBadReport'][0]) - laser.State = laser_info_runh['MReportLaser.MReport.fState'][0] - laser.IsOffsetCorrection = bool(laser_info_runh['MReportLaser.fIsOffsetCorrection'][0]) - laser.IsOffsetFitted = bool(laser_info_runh['MReportLaser.fIsOffsetFitted'][0]) - laser.IsBGCorrection = bool(laser_info_runh['MReportLaser.fIsBGCorrection'][0]) - laser.IsT0ShiftFitted = bool(laser_info_runh['MReportLaser.fIsT0ShiftFitted'][0]) - laser.IsUseGDAS = bool(laser_info_runh['MReportLaser.fIsUseGDAS'][0]) - laser.IsUpwardMoving = bool(laser_info_runh['MReportLaser.fIsUpwardMoving'][0]) - laser.OverShoot = int(laser_info_runh['MReportLaser.fOverShoot'][0]) - laser.UnderShoot = int(laser_info_runh['MReportLaser.fUnderShoot'][0]) - laser.BGSamples = int(laser_info_runh['MReportLaser.fBGSamples'][0]) - laser.Transmission3km = laser_info_runh['MReportLaser.fTransmission3km'][0] - laser.Transmission6km = laser_info_runh['MReportLaser.fTransmission6km'][0] - laser.Transmission9km = laser_info_runh['MReportLaser.fTransmission9km'][0] - laser.Transmission12km = laser_info_runh['MReportLaser.fTransmission12km'][0] - laser.Zenith = laser_info_runh['MReportLaser.fZenith'][0]* u.deg - laser.Azimuth = laser_info_runh['MReportLaser.fAzimuth'][0]* u.deg - laser.CloudFWHM = laser_info_runh['MReportLaser.fCloudFWHM[10]'][0] - laser.CloudBase = laser_info_runh['MReportLaser.fCloudBase[10]'][0] - laser.CloudTop = laser_info_runh['MReportLaser.fCloudTop[10]'][0] - laser.CloudTrans = laser_info_runh['MReportLaser.fCloudTrans[10]'][0] - laser.CloudHM = laser_info_runh['MReportLaser.fCloudHM[10]'][0] - laser.CloudHStd = laser_info_runh['MReportLaser.fCloudHStd[10]'][0] - laser.CloudLR = laser_info_runh['MReportLaser.fCloudLR[10]'][0] - laser.FullOverlap = laser_info_runh['MReportLaser.fFullOverlap'][0] - laser.EndGroundLayer = laser_info_runh['MReportLaser.fEndGroundLayer'][0] - laser.GroundLayerTrans = laser_info_runh['MReportLaser.fGroundLayerTrans'][0] - laser.Klett_k = laser_info_runh['MReportLaser.fKlett_k'][0] - laser.PheCounts = laser_info_runh['MReportLaser.fPheCounts'][0] - laser.Calimaness = laser_info_runh['MReportLaser.fCalimaness'][0] - laser.CloudLayerAlt = laser_info_runh['MReportLaser.fCloudLayerAlt'][0] - laser.CloudLayerDens = laser_info_runh['MReportLaser.fCloudLayerDens'][0] - laser.Offset = laser_info_runh['MReportLaser.fOffset'][0] - laser.Offset_Calculated = laser_info_runh['MReportLaser.fOffset_Calculated'][0] - laser.Offset_Fitted = laser_info_runh['MReportLaser.fOffset_Fitted'][0] - laser.Offset2 = laser_info_runh['MReportLaser.fOffset2'][0] - laser.Offset3 = laser_info_runh['MReportLaser.fOffset3'][0] - laser.Background1 = laser_info_runh['MReportLaser.fBackground1'][0] - laser.Background2 = laser_info_runh['MReportLaser.fBackground2'][0] - laser.BackgroundErr1 = laser_info_runh['MReportLaser.fBackgroundErr1'][0] - laser.BackgroundErr2 = laser_info_runh['MReportLaser.fBackgroundErr2'][0] - laser.RangeMax = laser_info_runh['MReportLaser.fRangeMax'][0] - laser.RangeMax_Clouds = laser_info_runh['MReportLaser.fRangeMax_Clouds'][0] - laser.ErrorCode = laser_info_runh['MReportLaser.fErrorCode'][0] - laser.ScaleHeight_fit = laser_info_runh['MReportLaser.fScaleHeight_fit'][0] - laser.Alpha_fit = laser_info_runh['MReportLaser.fChi2Alpha_fit'][0] - laser.Chi2Alpha_fit = laser_info_runh['MReportLaser.fChi2Alpha_fit'][0] - laser.Alpha_firstY = laser_info_runh['MReportLaser.fAlpha_firstY'][0] - laser.Alpha_Junge = laser_info_runh['MReportLaser.fAlpha_Junge'][0] - laser.PBLHeight = laser_info_runh['MReportLaser.fPBLHeight'][0] - laser.Chi2Full_fit = laser_info_runh['MReportLaser.fChi2Full_fit'][0] - laser.HWSwitchMaxOffset = laser_info_runh['MReportLaser.fHWSwitchMaxOffset'][0] - laser.NCollapse = laser_info_runh['MReportLaser.fNCollapse'][0] - laser.Shots = laser_info_runh['MReportLaser.fShots'][0] - laser.T0Shift = laser_info_runh['MReportLaser.fT0Shift'][0] - laser.Interval_0 = laser_info_runh['MReportLaser.fInterval_0'][0] - laser.RCS_min_perfect = laser_info_runh['MReportLaser.fRCS_min_perfect'][0] - laser.RCS_min_clouds = laser_info_runh['MReportLaser.fRCS_min_clouds'][0] - laser.RCS_min_mol = laser_info_runh['MReportLaser.fRCS_min_mol'][0] - laser.LIDAR_ratio = laser_info_runh['MReportLaser.fLIDAR_ratio'][0] - laser.LIDAR_ratio_Cloud = laser_info_runh['MReportLaser.fLIDAR_ratio_Cloud'][0] - laser.LIDAR_ratio_Junge = laser_info_runh['MReportLaser.fLIDAR_ratio_Junge'][0] - - millisec_seconds = millisec_value * 1e-3 - combined_mjd_value = mjd_value + millisec_seconds / 86400 - laser.MJD = combined_mjd_value + laser.UniqueID = laser_info_runh['MReportLaser.MReport.fUniqueID'] + laser.Bits = laser_info_runh['MReportLaser.MReport.fBits'] + mjd_value = laser_info_runh['MTimeLaser.fMjd'] + millisec_value = laser_info_runh['MTimeLaser.fTime.fMilliSec'] + laser.BadReport = bool(laser_info_runh['MReportLaser.MReport.fBadReport']) + laser.State = laser_info_runh['MReportLaser.MReport.fState'] + laser.IsOffsetCorrection = bool(laser_info_runh['MReportLaser.fIsOffsetCorrection']) + laser.IsOffsetFitted = bool(laser_info_runh['MReportLaser.fIsOffsetFitted']) + laser.IsBGCorrection = bool(laser_info_runh['MReportLaser.fIsBGCorrection']) + laser.IsT0ShiftFitted = bool(laser_info_runh['MReportLaser.fIsT0ShiftFitted']) + laser.IsUseGDAS = bool(laser_info_runh['MReportLaser.fIsUseGDAS']) + laser.IsUpwardMoving = bool(laser_info_runh['MReportLaser.fIsUpwardMoving']) + laser.OverShoot = int(laser_info_runh['MReportLaser.fOverShoot']) + laser.UnderShoot = int(laser_info_runh['MReportLaser.fUnderShoot']) + laser.BGSamples = int(laser_info_runh['MReportLaser.fBGSamples']) + laser.Transmission3km = laser_info_runh['MReportLaser.fTransmission3km'] + laser.Transmission6km = laser_info_runh['MReportLaser.fTransmission6km'] + laser.Transmission9km = laser_info_runh['MReportLaser.fTransmission9km'] + laser.Transmission12km = laser_info_runh['MReportLaser.fTransmission12km'] + laser.Zenith = laser_info_runh['MReportLaser.fZenith']* u.deg + laser.Azimuth = laser_info_runh['MReportLaser.fAzimuth']* u.deg + laser.CloudFWHM = laser_info_runh['MReportLaser.fCloudFWHM[10]'] + laser.CloudBase = laser_info_runh['MReportLaser.fCloudBase[10]'] + laser.CloudTop = laser_info_runh['MReportLaser.fCloudTop[10]'] + laser.CloudTrans = laser_info_runh['MReportLaser.fCloudTrans[10]'] + laser.CloudHM = laser_info_runh['MReportLaser.fCloudHM[10]'] + laser.CloudHStd = laser_info_runh['MReportLaser.fCloudHStd[10]'] + laser.CloudLR = laser_info_runh['MReportLaser.fCloudLR[10]'] + laser.FullOverlap = laser_info_runh['MReportLaser.fFullOverlap'] + laser.EndGroundLayer = laser_info_runh['MReportLaser.fEndGroundLayer'] + laser.GroundLayerTrans = laser_info_runh['MReportLaser.fGroundLayerTrans'] + laser.Klett_k = laser_info_runh['MReportLaser.fKlett_k'] + laser.PheCounts = laser_info_runh['MReportLaser.fPheCounts'] + laser.Calimaness = laser_info_runh['MReportLaser.fCalimaness'] + laser.CloudLayerAlt = laser_info_runh['MReportLaser.fCloudLayerAlt'] + laser.CloudLayerDens = laser_info_runh['MReportLaser.fCloudLayerDens'] + laser.Offset = laser_info_runh['MReportLaser.fOffset'] + laser.Offset_Calculated = laser_info_runh['MReportLaser.fOffset_Calculated'] + laser.Offset_Fitted = laser_info_runh['MReportLaser.fOffset_Fitted'] + laser.Offset2 = laser_info_runh['MReportLaser.fOffset2'] + laser.Offset3 = laser_info_runh['MReportLaser.fOffset3'] + laser.Background1 = laser_info_runh['MReportLaser.fBackground1'] + laser.Background2 = laser_info_runh['MReportLaser.fBackground2'] + laser.BackgroundErr1 = laser_info_runh['MReportLaser.fBackgroundErr1'] + laser.BackgroundErr2 = laser_info_runh['MReportLaser.fBackgroundErr2'] + laser.RangeMax = laser_info_runh['MReportLaser.fRangeMax'] + laser.RangeMax_Clouds = laser_info_runh['MReportLaser.fRangeMax_Clouds'] + laser.ErrorCode = laser_info_runh['MReportLaser.fErrorCode'] + laser.ScaleHeight_fit = laser_info_runh['MReportLaser.fScaleHeight_fit'] + laser.Alpha_fit = laser_info_runh['MReportLaser.fChi2Alpha_fit'] + laser.Chi2Alpha_fit = laser_info_runh['MReportLaser.fChi2Alpha_fit'] + laser.Alpha_firstY = laser_info_runh['MReportLaser.fAlpha_firstY'] + laser.Alpha_Junge = laser_info_runh['MReportLaser.fAlpha_Junge'] + laser.PBLHeight = laser_info_runh['MReportLaser.fPBLHeight'] + laser.Chi2Full_fit = laser_info_runh['MReportLaser.fChi2Full_fit'] + laser.HWSwitchMaxOffset = laser_info_runh['MReportLaser.fHWSwitchMaxOffset'] + laser.NCollapse = laser_info_runh['MReportLaser.fNCollapse'] + laser.Shots = laser_info_runh['MReportLaser.fShots'] + laser.T0Shift = laser_info_runh['MReportLaser.fT0Shift'] + laser.Interval_0 = laser_info_runh['MReportLaser.fInterval_0'] + laser.RCS_min_perfect = laser_info_runh['MReportLaser.fRCS_min_perfect'] + laser.RCS_min_clouds = laser_info_runh['MReportLaser.fRCS_min_clouds'] + laser.RCS_min_mol = laser_info_runh['MReportLaser.fRCS_min_mol'] + laser.LIDAR_ratio = laser_info_runh['MReportLaser.fLIDAR_ratio'] + laser.LIDAR_ratio_Cloud = laser_info_runh['MReportLaser.fLIDAR_ratio_Cloud'] + laser.LIDAR_ratio_Junge = laser_info_runh['MReportLaser.fLIDAR_ratio_Junge'] + + unique_reports = {} + + for report in reports: + mjd_value = report['fMjD'] + millisec_value = report['fMilliSec'] + millisec_seconds = millisec_value * 1e-3 + combined_mjd_value = mjd_value + millisec_seconds / 86400 + + if (mjd_value, millisec_value) not in unique_reports: + unique_reports[(mjd_value, millisec_value)] = combined_mjd_value + report['laser.MJD'] = combined_mjd_value except KeyError as e: print(f"Required key not found in the file {rootf}: {e}") From 6c424c83fd90d0647b14a7c4b6887a32e35902a6 Mon Sep 17 00:00:00 2001 From: nzywucka Date: Tue, 12 Mar 2024 09:29:41 +0100 Subject: [PATCH 13/38] Unindent that does not match any outer indentation level has been corrected. --- ctapipe_io_magic/__init__.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/ctapipe_io_magic/__init__.py b/ctapipe_io_magic/__init__.py index 5787cb7..310b77b 100644 --- a/ctapipe_io_magic/__init__.py +++ b/ctapipe_io_magic/__init__.py @@ -1115,9 +1115,9 @@ def parse_laser_info(self): millisec_seconds = millisec_value * 1e-3 combined_mjd_value = mjd_value + millisec_seconds / 86400 - if (mjd_value, millisec_value) not in unique_reports: - unique_reports[(mjd_value, millisec_value)] = combined_mjd_value - report['laser.MJD'] = combined_mjd_value + if (mjd_value, millisec_value) not in unique_reports: + unique_reports[(mjd_value, millisec_value)] = combined_mjd_value + report['laser.MJD'] = combined_mjd_value except KeyError as e: print(f"Required key not found in the file {rootf}: {e}") From f5e166931f2a230133412edcac04c7eadd1efc62 Mon Sep 17 00:00:00 2001 From: nzywucka Date: Tue, 12 Mar 2024 09:52:48 +0100 Subject: [PATCH 14/38] Filtering for the unique reports has been corrected. --- ctapipe_io_magic/__init__.py | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/ctapipe_io_magic/__init__.py b/ctapipe_io_magic/__init__.py index 310b77b..983654e 100644 --- a/ctapipe_io_magic/__init__.py +++ b/ctapipe_io_magic/__init__.py @@ -1034,7 +1034,7 @@ def parse_laser_info(self): ] laser = ReportLaserContainer() - + unique_reports = {} for rootf in self.files_: try: laser_info_runh = rootf['Laser'].arrays( @@ -1107,21 +1107,18 @@ def parse_laser_info(self): laser.LIDAR_ratio_Cloud = laser_info_runh['MReportLaser.fLIDAR_ratio_Cloud'] laser.LIDAR_ratio_Junge = laser_info_runh['MReportLaser.fLIDAR_ratio_Junge'] - unique_reports = {} - - for report in reports: - mjd_value = report['fMjD'] - millisec_value = report['fMilliSec'] + mjd_values = laser_info_runh['MTimeLaser.fMjd'] + millisec_values = laser_info_runh['MTimeLaser.fTime.fMilliSec'] + for mjd_value, millisec_value in zip(mjd_values, millisec_values): millisec_seconds = millisec_value * 1e-3 combined_mjd_value = mjd_value + millisec_seconds / 86400 - if (mjd_value, millisec_value) not in unique_reports: unique_reports[(mjd_value, millisec_value)] = combined_mjd_value - report['laser.MJD'] = combined_mjd_value - except KeyError as e: print(f"Required key not found in the file {rootf}: {e}") continue + for key, value in unique_reports.items(): + laser.MJD = value return laser From 27d8a8f6e2b12a2635408ca540906339281ed2d8 Mon Sep 17 00:00:00 2001 From: nzywucka Date: Tue, 12 Mar 2024 10:02:57 +0100 Subject: [PATCH 15/38] Unused variable has been removed. --- ctapipe_io_magic/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ctapipe_io_magic/__init__.py b/ctapipe_io_magic/__init__.py index 983654e..62f80dd 100644 --- a/ctapipe_io_magic/__init__.py +++ b/ctapipe_io_magic/__init__.py @@ -1117,7 +1117,7 @@ def parse_laser_info(self): except KeyError as e: print(f"Required key not found in the file {rootf}: {e}") continue - for key, value in unique_reports.items(): + for _ , value in unique_reports.items(): laser.MJD = value return laser From 4135a5ed3cc50c360c6041e10351a3a6eeb5389c Mon Sep 17 00:00:00 2001 From: nzywucka Date: Tue, 12 Mar 2024 11:28:58 +0100 Subject: [PATCH 16/38] Reading of the unnecesary parameters from the lidar reports are removed. Filtering the unique reports is improved. --- ctapipe_io_magic/__init__.py | 141 +++++++++++++++-------------------- 1 file changed, 61 insertions(+), 80 deletions(-) diff --git a/ctapipe_io_magic/__init__.py b/ctapipe_io_magic/__init__.py index 62f80dd..5cce5c5 100644 --- a/ctapipe_io_magic/__init__.py +++ b/ctapipe_io_magic/__init__.py @@ -11,6 +11,7 @@ import logging import numpy as np from numpy import nan +from typing import List from pathlib import Path from decimal import Decimal from astropy import units as u @@ -87,7 +88,7 @@ class ReportLaserContainer(Container): UniqueID = Field(None, 'No.') Bits = Field(None, 'ID') MJD = Field(None, 'Modified Julian Date') - BadReport = Field(bool, 'Bad Report') + BadReport = Field(bool, 'Bad Report') State = Field(None, 'State') IsOffsetCorrection = Field(bool, 'Is Offset Correction') IsOffsetFitted = Field(bool, 'Is Offset Fitted') @@ -95,62 +96,55 @@ class ReportLaserContainer(Container): IsT0ShiftFitted = Field(bool, 'Is T0 Shift Fitted') IsUseGDAS = Field(bool, 'Is Use GDAS') IsUpwardMoving = Field(bool, 'Is Upward Moving') - OverShoot = Field(nan, 'Over Shoot') - UnderShoot = Field(nan, 'Under Shoot') - BGSamples = Field(np.float32(np.nan), 'BG Samples') - Transmission3km = Field(np.float32(np.nan), 'Transmission at 3 km') - Transmission6km = Field(np.float32(np.nan), 'Transmission at 6 km') - Transmission9km = Field(np.float32(np.nan), 'Transmission at 9 km') - Transmission12km = Field(np.float32(np.nan), 'Transmission at 12 km') - Zenith = Field(Angle(np.nan, u.deg), 'Zenith angle', unit=u.deg) - Azimuth = Field(Angle(np.nan, u.deg), 'Azimuth angle', unit=u.deg) - CloudFWHM = Field(np.float32(np.nan), 'Cloud FWHM') - CloudBase = Field(np.float32(np.nan), 'Cloud Base') - CloudTop = Field(np.float32(np.nan), 'Cloud Top') - CloudTrans = Field(np.float32(np.nan), 'Cloud Trans') - CloudHM = Field(np.float32(np.nan), 'Cloud HM') - CloudHStd = Field(np.float32(np.nan), 'Cloud HStd') - CloudLR = Field(np.float32(np.nan), 'Cloud LR') - FullOverlap = Field(np.float32(np.nan), 'Full Overlap') - EndGroundLayer = Field(np.float32(np.nan), 'End Ground Layer') - GroundLayerTrans = Field(np.float32(np.nan), 'Ground Layer Trans') - Calimaness = Field(np.float32(np.nan), 'Calimaness') - CloudLayerAlt = Field(np.float32(np.nan), 'Altitude of cloud layer') - CloudLayerDens = Field(np.float32(np.nan), 'Density of cloud layer') - Klett_k = Field(np.float32(np.nan), 'Klett k') - PheCounts = Field([], 'Phe Counts') - Offset = Field(np.float32(np.nan), 'Offset') - Offset_Calculated = Field(np.float32(np.nan), 'Offset calculated') - Offset_Fitted = Field(np.float32(np.nan), 'Offset fitted') - Offset2 = Field(np.float32(np.nan), 'Offset 2') - Offset3 = Field(np.float32(np.nan), 'Offset 3') - Background1 = Field(np.float32(np.nan), 'Background 1') - Background2 = Field(np.float32(np.nan), 'Background 2') - BackgroundErr1 = Field(np.float32(np.nan), 'Background error 1') - BackgroundErr2 = Field(np.float32(np.nan), 'Background error 2') - RangeMax = Field(np.float32(np.nan), 'Range max') - RangeMax_Clouds = Field(np.float32(np.nan), 'Range max clouds') + OverShoot = Field(np.nan, 'Over Shoot') + UnderShoot = Field(np.nan, 'Under Shoot') + BGSamples = Field(List[np.float32], 'BG Samples') + Transmission3km = Field(List[np.float32], 'Transmission at 3 km') + Transmission6km = Field(List[np.float32], 'Transmission at 6 km') + Transmission9km = Field(List[np.float32], 'Transmission at 9 km') + Transmission12km = Field(List[np.float32], 'Transmission at 12 km') + Zenith = Field(List[Angle], 'Zenith angle', unit=u.deg) + Azimuth = Field(List[Angle], 'Azimuth angle', unit=u.deg) + FullOverlap = Field(List[np.float32], 'Full Overlap') + EndGroundLayer = Field(List[np.float32], 'End Ground Layer') + GroundLayerTrans = Field(List[np.float32], 'Ground Layer Trans') + Calimaness = Field(List[np.float32], 'Calimaness') + CloudLayerAlt = Field(List[np.float32], 'Altitude of cloud layer') + CloudLayerDens = Field(List[np.float32], 'Density of cloud layer') + Klett_k = Field(List[np.float32], 'Klett k') + PheCounts = Field(List[np.float32], 'Phe Counts') + Offset = Field(List[np.float32], 'Offset') + Offset_Calculated = Field(List[np.float32], 'Offset calculated') + Offset_Fitted = Field(List[np.float32], 'Offset fitted') + Offset2 = Field(List[np.float32], 'Offset 2') + Offset3 = Field(List[np.float32], 'Offset 3') + Background1 = Field(List[np.float32], 'Background 1') + Background2 = Field(List[np.float32], 'Background 2') + BackgroundErr1 = Field(List[np.float32], 'Background error 1') + BackgroundErr2 = Field(List[np.float32], 'Background error 2') + RangeMax = Field(List[np.float32], 'Range max') + RangeMax_Clouds = Field(List[np.float32], 'Range max clouds') ErrorCode = Field(None, 'Error code') - ScaleHeight_fit = Field(np.float32(np.nan), 'Scale Height fit') - Alpha_fit = Field(np.float32(np.nan), 'Alpha fit') - Chi2Alpha_fit = Field(np.float32(np.nan), 'Chi2 Alpha fit') - Alpha_firstY = Field(np.float32(np.nan), 'Alpha first Y') - Alpha_Junge = Field(np.float32(np.nan), 'Alpha Junge') - PBLHeight = Field(np.float32(np.nan), 'PBL Height') - Chi2Full_fit = Field(np.float32(np.nan), 'Chi2 Full fit') - SignalSamples = Field(np.float32(np.nan), 'Signal Samples') - HWSwitch = Field(np.float32(np.nan), 'HW Switch') - HWSwitchMaxOffset = Field(np.float32(np.nan), 'HW Switch Max Offset') - NCollapse = Field(np.float32(np.nan), 'N Collapse') - Shots = Field(np.float32(np.nan), 'Shots') - T0Shift = Field(np.float32(np.nan), 'T0 Shift') - Interval_0 = Field(np.float32(np.nan), 'Interval 0') - RCS_min_perfect = Field(np.float32(np.nan), 'RCS min perfect') - RCS_min_clouds = Field(np.float32(np.nan), 'RCS min cloud') - RCS_min_mol = Field(np.float32(np.nan), 'RCS min mol') - LIDAR_ratio = Field(np.float32(np.nan), 'LIDAR ratio') - LIDAR_ratio_Cloud = Field(np.float32(np.nan), 'LIDAR ratio cloud') - LIDAR_ratio_Junge = Field(np.float32(np.nan), 'LIDAR ratio Junge') + ScaleHeight_fit = Field(List[np.float32], 'Scale Height fit') + Alpha_fit = Field(List[np.float32], 'Alpha fit') + Chi2Alpha_fit = Field(List[np.float32], 'Chi2 Alpha fit') + Alpha_firstY = Field(List[np.float32], 'Alpha first Y') + Alpha_Junge = Field(List[np.float32], 'Alpha Junge') + PBLHeight = Field(List[np.float32], 'PBL Height') + Chi2Full_fit = Field(List[np.float32], 'Chi2 Full fit') + SignalSamples = Field(List[np.float32], 'Signal Samples') + HWSwitch = Field(List[np.float32], 'HW Switch') + HWSwitchMaxOffset = Field(List[np.float32], 'HW Switch Max Offset') + NCollapse = Field(List[np.float32], 'N Collapse') + Shots = Field(List[np.float32], 'Shots') + T0Shift = Field(List[np.float32], 'T0 Shift') + Interval_0 = Field(List[np.float32], 'Interval 0') + RCS_min_perfect = Field(List[np.float32], 'RCS min perfect') + RCS_min_clouds = Field(List[np.float32], 'RCS min cloud') + RCS_min_mol = Field(List[np.float32], 'RCS min mol') + LIDAR_ratio = Field(List[np.float32], 'LIDAR ratio') + LIDAR_ratio_Cloud = Field(List[np.float32], 'LIDAR ratio cloud') + LIDAR_ratio_Junge = Field(List[np.float32], 'LIDAR ratio Junge') def load_camera_geometry(): """Load camera geometry from bundled resources of this repo""" @@ -986,13 +980,6 @@ def parse_laser_info(self): 'MReportLaser.fTransmission12km', 'MReportLaser.fZenith', 'MReportLaser.fAzimuth', - 'MReportLaser.fCloudFWHM[10]', - 'MReportLaser.fCloudBase[10]', - 'MReportLaser.fCloudTop[10]', - 'MReportLaser.fCloudTrans[10]', - 'MReportLaser.fCloudHM[10]', - 'MReportLaser.fCloudHStd[10]', - 'MReportLaser.fCloudLR[10]', 'MReportLaser.fFullOverlap', 'MReportLaser.fEndGroundLayer', 'MReportLaser.fGroundLayerTrans', @@ -1061,13 +1048,6 @@ def parse_laser_info(self): laser.Transmission12km = laser_info_runh['MReportLaser.fTransmission12km'] laser.Zenith = laser_info_runh['MReportLaser.fZenith']* u.deg laser.Azimuth = laser_info_runh['MReportLaser.fAzimuth']* u.deg - laser.CloudFWHM = laser_info_runh['MReportLaser.fCloudFWHM[10]'] - laser.CloudBase = laser_info_runh['MReportLaser.fCloudBase[10]'] - laser.CloudTop = laser_info_runh['MReportLaser.fCloudTop[10]'] - laser.CloudTrans = laser_info_runh['MReportLaser.fCloudTrans[10]'] - laser.CloudHM = laser_info_runh['MReportLaser.fCloudHM[10]'] - laser.CloudHStd = laser_info_runh['MReportLaser.fCloudHStd[10]'] - laser.CloudLR = laser_info_runh['MReportLaser.fCloudLR[10]'] laser.FullOverlap = laser_info_runh['MReportLaser.fFullOverlap'] laser.EndGroundLayer = laser_info_runh['MReportLaser.fEndGroundLayer'] laser.GroundLayerTrans = laser_info_runh['MReportLaser.fGroundLayerTrans'] @@ -1107,19 +1087,20 @@ def parse_laser_info(self): laser.LIDAR_ratio_Cloud = laser_info_runh['MReportLaser.fLIDAR_ratio_Cloud'] laser.LIDAR_ratio_Junge = laser_info_runh['MReportLaser.fLIDAR_ratio_Junge'] - mjd_values = laser_info_runh['MTimeLaser.fMjd'] - millisec_values = laser_info_runh['MTimeLaser.fTime.fMilliSec'] - for mjd_value, millisec_value in zip(mjd_values, millisec_values): - millisec_seconds = millisec_value * 1e-3 - combined_mjd_value = mjd_value + millisec_seconds / 86400 - if (mjd_value, millisec_value) not in unique_reports: - unique_reports[(mjd_value, millisec_value)] = combined_mjd_value + for i in range(len(mjd_values)): + mjd_value = mjd_values[i] + millisec_value = millisec_values[i] + unique_key = (mjd_value, millisec_value) + if unique_key not in unique_reports: + unique_reports[unique_key] = True # Marking as seen + laser.UniqueID = laser_info_runh['MReportLaser.MReport.fUniqueID'][i] + laser.Bits = laser_info_runh['MReportLaser.MReport.fBits'][i] + millisec_seconds = millisec_value * 1e-3 + combined_mjd_value = mjd_value + millisec_seconds / 86400 + laser.MJD = combined_mjd_value except KeyError as e: print(f"Required key not found in the file {rootf}: {e}") continue - for _ , value in unique_reports.items(): - laser.MJD = value - return laser def parse_simulation_header(self): From 7d6e7cec0309bdcc34e7f019569f347b126c9a9d Mon Sep 17 00:00:00 2001 From: nzywucka Date: Tue, 12 Mar 2024 11:41:16 +0100 Subject: [PATCH 17/38] Unused import has been removed Reading the unique reports has been improved --- ctapipe_io_magic/__init__.py | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/ctapipe_io_magic/__init__.py b/ctapipe_io_magic/__init__.py index 5cce5c5..f8345d8 100644 --- a/ctapipe_io_magic/__init__.py +++ b/ctapipe_io_magic/__init__.py @@ -10,7 +10,6 @@ import uproot import logging import numpy as np -from numpy import nan from typing import List from pathlib import Path from decimal import Decimal @@ -1087,17 +1086,17 @@ def parse_laser_info(self): laser.LIDAR_ratio_Cloud = laser_info_runh['MReportLaser.fLIDAR_ratio_Cloud'] laser.LIDAR_ratio_Junge = laser_info_runh['MReportLaser.fLIDAR_ratio_Junge'] - for i in range(len(mjd_values)): - mjd_value = mjd_values[i] - millisec_value = millisec_values[i] - unique_key = (mjd_value, millisec_value) + for idx, (mjd_values, millisec_values) in enumerate(zip(mjd_value, millisec_value)): + unique_key = (mjd_values, millisec_values) + if unique_key not in unique_reports: - unique_reports[unique_key] = True # Marking as seen - laser.UniqueID = laser_info_runh['MReportLaser.MReport.fUniqueID'][i] - laser.Bits = laser_info_runh['MReportLaser.MReport.fBits'][i] - millisec_seconds = millisec_value * 1e-3 - combined_mjd_value = mjd_value + millisec_seconds / 86400 + unique_reports[unique_key] = True + laser.UniqueID = laser_info_runh['MReportLaser.MReport.fUniqueID'][idx] + laser.Bits = laser_info_runh['MReportLaser.MReport.fBits'][idx] + millisec_seconds = millisec_values * 1e-3 + combined_mjd_value = mjd_values + millisec_seconds / 86400 laser.MJD = combined_mjd_value + except KeyError as e: print(f"Required key not found in the file {rootf}: {e}") continue From ec4c4b7744da2d63901309e58a2924206e4ffc0c Mon Sep 17 00:00:00 2001 From: nzywucka Date: Tue, 12 Mar 2024 11:48:43 +0100 Subject: [PATCH 18/38] Trailing whitespace has been removed --- ctapipe_io_magic/__init__.py | 1 - 1 file changed, 1 deletion(-) diff --git a/ctapipe_io_magic/__init__.py b/ctapipe_io_magic/__init__.py index f8345d8..f574c86 100644 --- a/ctapipe_io_magic/__init__.py +++ b/ctapipe_io_magic/__init__.py @@ -1088,7 +1088,6 @@ def parse_laser_info(self): for idx, (mjd_values, millisec_values) in enumerate(zip(mjd_value, millisec_value)): unique_key = (mjd_values, millisec_values) - if unique_key not in unique_reports: unique_reports[unique_key] = True laser.UniqueID = laser_info_runh['MReportLaser.MReport.fUniqueID'][idx] From 2b2dd17a1e1fa2f809a03433b709a3ca4cb8b676 Mon Sep 17 00:00:00 2001 From: nzywucka Date: Tue, 12 Mar 2024 14:28:35 +0100 Subject: [PATCH 19/38] Reading the unique reports has been improved to keep all the unique parameters in the container. All the parameters in the container are stored in the list of parameters now. --- ctapipe_io_magic/__init__.py | 39 ++++++++++++++++++------------------ 1 file changed, 20 insertions(+), 19 deletions(-) diff --git a/ctapipe_io_magic/__init__.py b/ctapipe_io_magic/__init__.py index f574c86..4d4335d 100644 --- a/ctapipe_io_magic/__init__.py +++ b/ctapipe_io_magic/__init__.py @@ -84,19 +84,19 @@ class ReportLaserContainer(Container): """ Container for Magic laser parameters """ - UniqueID = Field(None, 'No.') - Bits = Field(None, 'ID') - MJD = Field(None, 'Modified Julian Date') - BadReport = Field(bool, 'Bad Report') - State = Field(None, 'State') - IsOffsetCorrection = Field(bool, 'Is Offset Correction') - IsOffsetFitted = Field(bool, 'Is Offset Fitted') - IsBGCorrection = Field(bool, 'Is BG Correction') - IsT0ShiftFitted = Field(bool, 'Is T0 Shift Fitted') - IsUseGDAS = Field(bool, 'Is Use GDAS') - IsUpwardMoving = Field(bool, 'Is Upward Moving') - OverShoot = Field(np.nan, 'Over Shoot') - UnderShoot = Field(np.nan, 'Under Shoot') + UniqueID = Field(List[], 'No.') + Bits = Field(List[], 'ID') + MJD = Field(List[], 'Modified Julian Date') + BadReport = Field(List[], 'Bad Report') + State = Field(List[], 'State') + IsOffsetCorrection = Field(List[], 'Is Offset Correction') + IsOffsetFitted = Field(List[], 'Is Offset Fitted') + IsBGCorrection = Field(List[], 'Is BG Correction') + IsT0ShiftFitted = Field(List[], 'Is T0 Shift Fitted') + IsUseGDAS = Field(List[], 'Is Use GDAS') + IsUpwardMoving = Field(List[], 'Is Upward Moving') + OverShoot = Field(List[], 'Over Shoot') + UnderShoot = Field(List[], 'Under Shoot') BGSamples = Field(List[np.float32], 'BG Samples') Transmission3km = Field(List[np.float32], 'Transmission at 3 km') Transmission6km = Field(List[np.float32], 'Transmission at 6 km') @@ -111,7 +111,7 @@ class ReportLaserContainer(Container): CloudLayerAlt = Field(List[np.float32], 'Altitude of cloud layer') CloudLayerDens = Field(List[np.float32], 'Density of cloud layer') Klett_k = Field(List[np.float32], 'Klett k') - PheCounts = Field(List[np.float32], 'Phe Counts') + PheCounts = Field(List[List[np.float32]], 'Phe Counts') Offset = Field(List[np.float32], 'Offset') Offset_Calculated = Field(List[np.float32], 'Offset calculated') Offset_Fitted = Field(List[np.float32], 'Offset fitted') @@ -123,7 +123,7 @@ class ReportLaserContainer(Container): BackgroundErr2 = Field(List[np.float32], 'Background error 2') RangeMax = Field(List[np.float32], 'Range max') RangeMax_Clouds = Field(List[np.float32], 'Range max clouds') - ErrorCode = Field(None, 'Error code') + ErrorCode = Field(List[], 'Error code') ScaleHeight_fit = Field(List[np.float32], 'Scale Height fit') Alpha_fit = Field(List[np.float32], 'Alpha fit') Chi2Alpha_fit = Field(List[np.float32], 'Chi2 Alpha fit') @@ -1090,12 +1090,13 @@ def parse_laser_info(self): unique_key = (mjd_values, millisec_values) if unique_key not in unique_reports: unique_reports[unique_key] = True - laser.UniqueID = laser_info_runh['MReportLaser.MReport.fUniqueID'][idx] - laser.Bits = laser_info_runh['MReportLaser.MReport.fBits'][idx] + unique_laser = ReportLaserContainer() + unique_laser.UniqueID = laser_info_runh['MReportLaser.MReport.fUniqueID'][idx] + unique_laser.Bits = laser_info_runh['MReportLaser.MReport.fBits'][idx] millisec_seconds = millisec_values * 1e-3 combined_mjd_value = mjd_values + millisec_seconds / 86400 - laser.MJD = combined_mjd_value - + unique_laser.MJD = combined_mjd_value + laser.append(unique_laser) except KeyError as e: print(f"Required key not found in the file {rootf}: {e}") continue From 08b5101a211c7ed4d96da6bbe71446c4fd2e3eee Mon Sep 17 00:00:00 2001 From: nzywucka Date: Tue, 12 Mar 2024 14:44:12 +0100 Subject: [PATCH 20/38] The definitions of the container parameters have been improved. --- ctapipe_io_magic/__init__.py | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/ctapipe_io_magic/__init__.py b/ctapipe_io_magic/__init__.py index 4d4335d..1fbc0ac 100644 --- a/ctapipe_io_magic/__init__.py +++ b/ctapipe_io_magic/__init__.py @@ -10,7 +10,7 @@ import uproot import logging import numpy as np -from typing import List +from typing import List, Any from pathlib import Path from decimal import Decimal from astropy import units as u @@ -84,19 +84,19 @@ class ReportLaserContainer(Container): """ Container for Magic laser parameters """ - UniqueID = Field(List[], 'No.') - Bits = Field(List[], 'ID') - MJD = Field(List[], 'Modified Julian Date') - BadReport = Field(List[], 'Bad Report') - State = Field(List[], 'State') - IsOffsetCorrection = Field(List[], 'Is Offset Correction') - IsOffsetFitted = Field(List[], 'Is Offset Fitted') - IsBGCorrection = Field(List[], 'Is BG Correction') - IsT0ShiftFitted = Field(List[], 'Is T0 Shift Fitted') - IsUseGDAS = Field(List[], 'Is Use GDAS') - IsUpwardMoving = Field(List[], 'Is Upward Moving') - OverShoot = Field(List[], 'Over Shoot') - UnderShoot = Field(List[], 'Under Shoot') + UniqueID = Field(List[Any], 'No.') + Bits = Field(List[Any], 'ID') + MJD = Field(List[Any], 'Modified Julian Date') + BadReport = Field(List[bool], 'Bad Report') + State = Field(List[Any], 'State') + IsOffsetCorrection = Field(List[bool], 'Is Offset Correction') + IsOffsetFitted = Field(List[bool], 'Is Offset Fitted') + IsBGCorrection = Field(List[bool], 'Is BG Correction') + IsT0ShiftFitted = Field(List[bool], 'Is T0 Shift Fitted') + IsUseGDAS = Field(List[bool], 'Is Use GDAS') + IsUpwardMoving = Field(List[bool], 'Is Upward Moving') + OverShoot = Field(List[Any], 'Over Shoot') + UnderShoot = Field(List[Any], 'Under Shoot') BGSamples = Field(List[np.float32], 'BG Samples') Transmission3km = Field(List[np.float32], 'Transmission at 3 km') Transmission6km = Field(List[np.float32], 'Transmission at 6 km') @@ -123,7 +123,7 @@ class ReportLaserContainer(Container): BackgroundErr2 = Field(List[np.float32], 'Background error 2') RangeMax = Field(List[np.float32], 'Range max') RangeMax_Clouds = Field(List[np.float32], 'Range max clouds') - ErrorCode = Field(List[], 'Error code') + ErrorCode = Field(List[Any], 'Error code') ScaleHeight_fit = Field(List[np.float32], 'Scale Height fit') Alpha_fit = Field(List[np.float32], 'Alpha fit') Chi2Alpha_fit = Field(List[np.float32], 'Chi2 Alpha fit') From f1024fae12ccc59d5504c2c755ae2eb73e5e66c0 Mon Sep 17 00:00:00 2001 From: nzywucka Date: Tue, 12 Mar 2024 15:08:08 +0100 Subject: [PATCH 21/38] Appending the unique parameters to the container was improved. --- ctapipe_io_magic/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ctapipe_io_magic/__init__.py b/ctapipe_io_magic/__init__.py index 1fbc0ac..abd11f3 100644 --- a/ctapipe_io_magic/__init__.py +++ b/ctapipe_io_magic/__init__.py @@ -1096,7 +1096,7 @@ def parse_laser_info(self): millisec_seconds = millisec_values * 1e-3 combined_mjd_value = mjd_values + millisec_seconds / 86400 unique_laser.MJD = combined_mjd_value - laser.append(unique_laser) + unique_reports.append(unique_laser) except KeyError as e: print(f"Required key not found in the file {rootf}: {e}") continue From 4625c5d5727f6c0f8b466abc681060f0ac9e3434 Mon Sep 17 00:00:00 2001 From: nzywucka Date: Wed, 13 Mar 2024 09:29:13 +0100 Subject: [PATCH 22/38] AttributeError has been fixed. --- ctapipe_io_magic/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ctapipe_io_magic/__init__.py b/ctapipe_io_magic/__init__.py index abd11f3..d35cab6 100644 --- a/ctapipe_io_magic/__init__.py +++ b/ctapipe_io_magic/__init__.py @@ -1096,7 +1096,7 @@ def parse_laser_info(self): millisec_seconds = millisec_values * 1e-3 combined_mjd_value = mjd_values + millisec_seconds / 86400 unique_laser.MJD = combined_mjd_value - unique_reports.append(unique_laser) + unique_reports[unique_key] = unique_laser except KeyError as e: print(f"Required key not found in the file {rootf}: {e}") continue From bb145237475a9b300bfb4eb497506cfd5aeb0586 Mon Sep 17 00:00:00 2001 From: nzywucka Date: Mon, 25 Mar 2024 19:05:01 +0100 Subject: [PATCH 23/38] Filtering of lidar parameters has been improved. --- ctapipe_io_magic/__init__.py | 152 +++++++++--------- .../tests/test_magic_event_source.py | 57 ++++++- 2 files changed, 132 insertions(+), 77 deletions(-) diff --git a/ctapipe_io_magic/__init__.py b/ctapipe_io_magic/__init__.py index 837e1c3..81eb900 100644 --- a/ctapipe_io_magic/__init__.py +++ b/ctapipe_io_magic/__init__.py @@ -99,14 +99,14 @@ class ReportLaserContainer(Container): UniqueID = Field(List[Any], 'No.') Bits = Field(List[Any], 'ID') MJD = Field(List[Any], 'Modified Julian Date') - BadReport = Field(List[bool], 'Bad Report') - State = Field(List[Any], 'State') - IsOffsetCorrection = Field(List[bool], 'Is Offset Correction') - IsOffsetFitted = Field(List[bool], 'Is Offset Fitted') - IsBGCorrection = Field(List[bool], 'Is BG Correction') - IsT0ShiftFitted = Field(List[bool], 'Is T0 Shift Fitted') - IsUseGDAS = Field(List[bool], 'Is Use GDAS') - IsUpwardMoving = Field(List[bool], 'Is Upward Moving') +# BadReport = Field(List[Any], 'Bad Report') # to be improved +# State = Field(List[Any], 'State') # to be improved + IsOffsetCorrection = Field(List[Any], 'Is Offset Correction') + IsOffsetFitted = Field(List[Any], 'Is Offset Fitted') + IsBGCorrection = Field(List[Any], 'Is BG Correction') + IsT0ShiftFitted = Field(List[Any], 'Is T0 Shift Fitted') + IsUseGDAS = Field(List[Any], 'Is Use GDAS') + IsUpwardMoving = Field(List[Any], 'Is Upward Moving') OverShoot = Field(List[Any], 'Over Shoot') UnderShoot = Field(List[Any], 'Under Shoot') BGSamples = Field(List[np.float32], 'BG Samples') @@ -448,13 +448,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: @@ -1108,8 +1108,8 @@ def parse_laser_info(self): 'MReportLaser.MReport.fBits', 'MTimeLaser.fMjd', 'MTimeLaser.fTime.fMilliSec', - 'MReportLaser.MReport.fBadReport', - 'MReportLaser.MReport.fState', + #'MReportLaser.MReport.fBadReport', # to be improved + #'MReportLaser.MReport.fState', # to be improved 'MReportLaser.fIsOffsetCorrection', 'MReportLaser.fIsOffsetFitted', 'MReportLaser.fIsBGCorrection', @@ -1165,84 +1165,84 @@ def parse_laser_info(self): 'MReportLaser.fLIDAR_ratio_Junge', ] - laser = ReportLaserContainer() unique_reports = {} for rootf in self.files_: try: laser_info_runh = rootf['Laser'].arrays( - laser_info_array_list_runh, library="np" + laser_info_array_list_runh, library="np" ) - laser.UniqueID = laser_info_runh['MReportLaser.MReport.fUniqueID'] - laser.Bits = laser_info_runh['MReportLaser.MReport.fBits'] mjd_value = laser_info_runh['MTimeLaser.fMjd'] millisec_value = laser_info_runh['MTimeLaser.fTime.fMilliSec'] - laser.BadReport = bool(laser_info_runh['MReportLaser.MReport.fBadReport']) - laser.State = laser_info_runh['MReportLaser.MReport.fState'] - laser.IsOffsetCorrection = bool(laser_info_runh['MReportLaser.fIsOffsetCorrection']) - laser.IsOffsetFitted = bool(laser_info_runh['MReportLaser.fIsOffsetFitted']) - laser.IsBGCorrection = bool(laser_info_runh['MReportLaser.fIsBGCorrection']) - laser.IsT0ShiftFitted = bool(laser_info_runh['MReportLaser.fIsT0ShiftFitted']) - laser.IsUseGDAS = bool(laser_info_runh['MReportLaser.fIsUseGDAS']) - laser.IsUpwardMoving = bool(laser_info_runh['MReportLaser.fIsUpwardMoving']) - laser.OverShoot = int(laser_info_runh['MReportLaser.fOverShoot']) - laser.UnderShoot = int(laser_info_runh['MReportLaser.fUnderShoot']) - laser.BGSamples = int(laser_info_runh['MReportLaser.fBGSamples']) - laser.Transmission3km = laser_info_runh['MReportLaser.fTransmission3km'] - laser.Transmission6km = laser_info_runh['MReportLaser.fTransmission6km'] - laser.Transmission9km = laser_info_runh['MReportLaser.fTransmission9km'] - laser.Transmission12km = laser_info_runh['MReportLaser.fTransmission12km'] - laser.Zenith = laser_info_runh['MReportLaser.fZenith']* u.deg - laser.Azimuth = laser_info_runh['MReportLaser.fAzimuth']* u.deg - laser.FullOverlap = laser_info_runh['MReportLaser.fFullOverlap'] - laser.EndGroundLayer = laser_info_runh['MReportLaser.fEndGroundLayer'] - laser.GroundLayerTrans = laser_info_runh['MReportLaser.fGroundLayerTrans'] - laser.Klett_k = laser_info_runh['MReportLaser.fKlett_k'] - laser.PheCounts = laser_info_runh['MReportLaser.fPheCounts'] - laser.Calimaness = laser_info_runh['MReportLaser.fCalimaness'] - laser.CloudLayerAlt = laser_info_runh['MReportLaser.fCloudLayerAlt'] - laser.CloudLayerDens = laser_info_runh['MReportLaser.fCloudLayerDens'] - laser.Offset = laser_info_runh['MReportLaser.fOffset'] - laser.Offset_Calculated = laser_info_runh['MReportLaser.fOffset_Calculated'] - laser.Offset_Fitted = laser_info_runh['MReportLaser.fOffset_Fitted'] - laser.Offset2 = laser_info_runh['MReportLaser.fOffset2'] - laser.Offset3 = laser_info_runh['MReportLaser.fOffset3'] - laser.Background1 = laser_info_runh['MReportLaser.fBackground1'] - laser.Background2 = laser_info_runh['MReportLaser.fBackground2'] - laser.BackgroundErr1 = laser_info_runh['MReportLaser.fBackgroundErr1'] - laser.BackgroundErr2 = laser_info_runh['MReportLaser.fBackgroundErr2'] - laser.RangeMax = laser_info_runh['MReportLaser.fRangeMax'] - laser.RangeMax_Clouds = laser_info_runh['MReportLaser.fRangeMax_Clouds'] - laser.ErrorCode = laser_info_runh['MReportLaser.fErrorCode'] - laser.ScaleHeight_fit = laser_info_runh['MReportLaser.fScaleHeight_fit'] - laser.Alpha_fit = laser_info_runh['MReportLaser.fChi2Alpha_fit'] - laser.Chi2Alpha_fit = laser_info_runh['MReportLaser.fChi2Alpha_fit'] - laser.Alpha_firstY = laser_info_runh['MReportLaser.fAlpha_firstY'] - laser.Alpha_Junge = laser_info_runh['MReportLaser.fAlpha_Junge'] - laser.PBLHeight = laser_info_runh['MReportLaser.fPBLHeight'] - laser.Chi2Full_fit = laser_info_runh['MReportLaser.fChi2Full_fit'] - laser.HWSwitchMaxOffset = laser_info_runh['MReportLaser.fHWSwitchMaxOffset'] - laser.NCollapse = laser_info_runh['MReportLaser.fNCollapse'] - laser.Shots = laser_info_runh['MReportLaser.fShots'] - laser.T0Shift = laser_info_runh['MReportLaser.fT0Shift'] - laser.Interval_0 = laser_info_runh['MReportLaser.fInterval_0'] - laser.RCS_min_perfect = laser_info_runh['MReportLaser.fRCS_min_perfect'] - laser.RCS_min_clouds = laser_info_runh['MReportLaser.fRCS_min_clouds'] - laser.RCS_min_mol = laser_info_runh['MReportLaser.fRCS_min_mol'] - laser.LIDAR_ratio = laser_info_runh['MReportLaser.fLIDAR_ratio'] - laser.LIDAR_ratio_Cloud = laser_info_runh['MReportLaser.fLIDAR_ratio_Cloud'] - laser.LIDAR_ratio_Junge = laser_info_runh['MReportLaser.fLIDAR_ratio_Junge'] for idx, (mjd_values, millisec_values) in enumerate(zip(mjd_value, millisec_value)): unique_key = (mjd_values, millisec_values) if unique_key not in unique_reports: unique_reports[unique_key] = True - unique_laser = ReportLaserContainer() - unique_laser.UniqueID = laser_info_runh['MReportLaser.MReport.fUniqueID'][idx] - unique_laser.Bits = laser_info_runh['MReportLaser.MReport.fBits'][idx] + laser = ReportLaserContainer() + laser.UniqueID = laser_info_runh['MReportLaser.MReport.fUniqueID'] + laser.Bits = laser_info_runh['MReportLaser.MReport.fBits'] + #mjd_value = laser_info_runh['MTimeLaser.fMjd'] + #millisec_value = laser_info_runh['MTimeLaser.fTime.fMilliSec'] + #laser.BadReport = laser_info_runh['MReportLaser.MReport.fBadReport'] # to be improved + #laser.State = laser_info_runh['MReportLaser.MReport.fState'] # to be improved + laser.IsOffsetCorrection = laser_info_runh['MReportLaser.fIsOffsetCorrection'] + laser.IsOffsetFitted = laser_info_runh['MReportLaser.fIsOffsetFitted'] + laser.IsBGCorrection = laser_info_runh['MReportLaser.fIsBGCorrection'] + laser.IsT0ShiftFitted = laser_info_runh['MReportLaser.fIsT0ShiftFitted'] + laser.IsUseGDAS = laser_info_runh['MReportLaser.fIsUseGDAS'] + laser.IsUpwardMoving = laser_info_runh['MReportLaser.fIsUpwardMoving'] + laser.OverShoot = int(laser_info_runh['MReportLaser.fOverShoot']) + laser.UnderShoot = int(laser_info_runh['MReportLaser.fUnderShoot']) + laser.BGSamples = int(laser_info_runh['MReportLaser.fBGSamples']) + laser.Transmission3km = laser_info_runh['MReportLaser.fTransmission3km'] + laser.Transmission6km = laser_info_runh['MReportLaser.fTransmission6km'] + laser.Transmission9km = laser_info_runh['MReportLaser.fTransmission9km'] + laser.Transmission12km = laser_info_runh['MReportLaser.fTransmission12km'] + laser.Zenith = laser_info_runh['MReportLaser.fZenith']* u.deg + laser.Azimuth = laser_info_runh['MReportLaser.fAzimuth']* u.deg + laser.FullOverlap = laser_info_runh['MReportLaser.fFullOverlap'] + laser.EndGroundLayer = laser_info_runh['MReportLaser.fEndGroundLayer'] + laser.GroundLayerTrans = laser_info_runh['MReportLaser.fGroundLayerTrans'] + laser.Klett_k = laser_info_runh['MReportLaser.fKlett_k'] + laser.PheCounts = laser_info_runh['MReportLaser.fPheCounts'] + laser.Calimaness = laser_info_runh['MReportLaser.fCalimaness'] + laser.CloudLayerAlt = laser_info_runh['MReportLaser.fCloudLayerAlt'] + laser.CloudLayerDens = laser_info_runh['MReportLaser.fCloudLayerDens'] + laser.Offset = laser_info_runh['MReportLaser.fOffset'] + laser.Offset_Calculated = laser_info_runh['MReportLaser.fOffset_Calculated'] + laser.Offset_Fitted = laser_info_runh['MReportLaser.fOffset_Fitted'] + laser.Offset2 = laser_info_runh['MReportLaser.fOffset2'] + laser.Offset3 = laser_info_runh['MReportLaser.fOffset3'] + laser.Background1 = laser_info_runh['MReportLaser.fBackground1'] + laser.Background2 = laser_info_runh['MReportLaser.fBackground2'] + laser.BackgroundErr1 = laser_info_runh['MReportLaser.fBackgroundErr1'] + laser.BackgroundErr2 = laser_info_runh['MReportLaser.fBackgroundErr2'] + laser.RangeMax = laser_info_runh['MReportLaser.fRangeMax'] + laser.RangeMax_Clouds = laser_info_runh['MReportLaser.fRangeMax_Clouds'] + laser.ErrorCode = laser_info_runh['MReportLaser.fErrorCode'] + laser.ScaleHeight_fit = laser_info_runh['MReportLaser.fScaleHeight_fit'] + laser.Alpha_fit = laser_info_runh['MReportLaser.fChi2Alpha_fit'] + laser.Chi2Alpha_fit = laser_info_runh['MReportLaser.fChi2Alpha_fit'] + laser.Alpha_firstY = laser_info_runh['MReportLaser.fAlpha_firstY'] + laser.Alpha_Junge = laser_info_runh['MReportLaser.fAlpha_Junge'] + laser.PBLHeight = laser_info_runh['MReportLaser.fPBLHeight'] + laser.Chi2Full_fit = laser_info_runh['MReportLaser.fChi2Full_fit'] + laser.HWSwitchMaxOffset = laser_info_runh['MReportLaser.fHWSwitchMaxOffset'] + laser.NCollapse = laser_info_runh['MReportLaser.fNCollapse'] + laser.Shots = laser_info_runh['MReportLaser.fShots'] + laser.T0Shift = laser_info_runh['MReportLaser.fT0Shift'] + laser.Interval_0 = laser_info_runh['MReportLaser.fInterval_0'] + laser.RCS_min_perfect = laser_info_runh['MReportLaser.fRCS_min_perfect'] + laser.RCS_min_clouds = laser_info_runh['MReportLaser.fRCS_min_clouds'] + laser.RCS_min_mol = laser_info_runh['MReportLaser.fRCS_min_mol'] + laser.LIDAR_ratio = laser_info_runh['MReportLaser.fLIDAR_ratio'] + laser.LIDAR_ratio_Cloud = laser_info_runh['MReportLaser.fLIDAR_ratio_Cloud'] + laser.LIDAR_ratio_Junge = laser_info_runh['MReportLaser.fLIDAR_ratio_Junge'] + millisec_seconds = millisec_values * 1e-3 combined_mjd_value = mjd_values + millisec_seconds / 86400 - unique_laser.MJD = combined_mjd_value - unique_reports[unique_key] = unique_laser + laser.MJD = combined_mjd_value + unique_reports[unique_key] = laser except KeyError as e: print(f"Required key not found in the file {rootf}: {e}") continue diff --git a/ctapipe_io_magic/tests/test_magic_event_source.py b/ctapipe_io_magic/tests/test_magic_event_source.py index fbac975..1b46c86 100644 --- a/ctapipe_io_magic/tests/test_magic_event_source.py +++ b/ctapipe_io_magic/tests/test_magic_event_source.py @@ -286,7 +286,7 @@ def test_allowed_tels(): dataset = ( test_calibrated_real_dir - / "20210314_M1_05095172.001_Y_CrabNebula-W0.40+035.root" + / "20210314_M1_05095172.002_Y_CrabNebula-W0.40+035.root" #"20210314_M1_05095172.001_Y_CrabNebula-W0.40+035.root" ) allowed_tels = {1} with MAGICEventSource( @@ -559,3 +559,58 @@ def test_broken_subruns_missing_arrays(): input_url=input_file, process_run=True, ) +# def test_eventseeker(): +# dataset = get_dataset_path("20131004_M1_05029747.003_Y_MagicCrab-W0.40+035.root") +# +# with MAGICEventSource(input_url=dataset) as source: +# seeker = EventSeeker(source) +# event = seeker.get_event_index(0) +# assert event.count == 0 +# assert event.index.event_id == 29795 +# +# event = seeker.get_event_index(2) +# assert event.count == 2 +# assert event.index.event_id == 29798 +# +# def test_eventcontent(): +# dataset = get_dataset_path("20131004_M1_05029747.003_Y_MagicCrab-W0.40+035.root") +# +# with MAGICEventSource(input_url=dataset) as source: +# seeker = EventSeeker(source) +# event = seeker.get_event_index(0) +# assert event.dl1.tel[1].image[0] == -0.53125 +# assert event.dl1.tel[1].peak_time[0] == 49.125 + + +#@pytest.mark.parametrize("dataset", test_calibrated_all) +#def test_lidar_parameters(dataset): +# try: +# with MAGICEventSource(input_url=dataset) as source: +# params = ReportLaserContainer() +# for item, event in params: +# #params.parse_laser_info() +# # assert params.laser.MJD == 59286.91349633102 +# assert params.Transmission12km == 0.89 +# # assert params.laser.Azimuth == u.Quantity(277, u.deg) +# #assert params.laser.BadReport == False +# except KeyInFileError: +# pytest.skip("Skipping test for file without 'Laser' key.") + +#dataset = ( +# test_calibrated_real_dir +# / "20210314_M1_05095172.002_Y_CrabNebula-W0.40+035.root" +#) + +@pytest.mark.parametrize("dataset", test_calibrated_all) +def test_lidar_parameters(dataset): + from ctapipe_io_magic import MAGICEventSource, ReportLaserContainer + + dataset = (test_calibrated_real_dir/"20210314_M1_05095172.002_Y_CrabNebula-W0.40+035.root") + + with MAGICEventSource(input_url=dataset) as source: + assert source.laser.Transmission12km == pytest.approx(0.89) + assert source.laser.Transmission3km == pytest.approx(0.96) + assert source.laser.Transmission6km == pytest.approx(0.93) + assert source.laser.Transmission9km == pytest.approx(0.89) + + From bcaaef2de6b3e2e0972490e68f932b94547585b5 Mon Sep 17 00:00:00 2001 From: nzywucka Date: Wed, 27 Mar 2024 10:51:21 +0100 Subject: [PATCH 24/38] Filtering of lidar parameters has been improved # Please enter the commit message for your changes. Lines starting --- ctapipe_io_magic/__init__.py | 134 +++++++++++++++++------------------ 1 file changed, 66 insertions(+), 68 deletions(-) diff --git a/ctapipe_io_magic/__init__.py b/ctapipe_io_magic/__init__.py index 81eb900..4d55c02 100644 --- a/ctapipe_io_magic/__init__.py +++ b/ctapipe_io_magic/__init__.py @@ -1174,79 +1174,77 @@ def parse_laser_info(self): mjd_value = laser_info_runh['MTimeLaser.fMjd'] millisec_value = laser_info_runh['MTimeLaser.fTime.fMilliSec'] - for idx, (mjd_values, millisec_values) in enumerate(zip(mjd_value, millisec_value)): + for mjd_values, millisec_values in zip(mjd_value, millisec_value): unique_key = (mjd_values, millisec_values) if unique_key not in unique_reports: - unique_reports[unique_key] = True - laser = ReportLaserContainer() - laser.UniqueID = laser_info_runh['MReportLaser.MReport.fUniqueID'] - laser.Bits = laser_info_runh['MReportLaser.MReport.fBits'] - #mjd_value = laser_info_runh['MTimeLaser.fMjd'] - #millisec_value = laser_info_runh['MTimeLaser.fTime.fMilliSec'] - #laser.BadReport = laser_info_runh['MReportLaser.MReport.fBadReport'] # to be improved - #laser.State = laser_info_runh['MReportLaser.MReport.fState'] # to be improved - laser.IsOffsetCorrection = laser_info_runh['MReportLaser.fIsOffsetCorrection'] - laser.IsOffsetFitted = laser_info_runh['MReportLaser.fIsOffsetFitted'] - laser.IsBGCorrection = laser_info_runh['MReportLaser.fIsBGCorrection'] - laser.IsT0ShiftFitted = laser_info_runh['MReportLaser.fIsT0ShiftFitted'] - laser.IsUseGDAS = laser_info_runh['MReportLaser.fIsUseGDAS'] - laser.IsUpwardMoving = laser_info_runh['MReportLaser.fIsUpwardMoving'] - laser.OverShoot = int(laser_info_runh['MReportLaser.fOverShoot']) - laser.UnderShoot = int(laser_info_runh['MReportLaser.fUnderShoot']) - laser.BGSamples = int(laser_info_runh['MReportLaser.fBGSamples']) - laser.Transmission3km = laser_info_runh['MReportLaser.fTransmission3km'] - laser.Transmission6km = laser_info_runh['MReportLaser.fTransmission6km'] - laser.Transmission9km = laser_info_runh['MReportLaser.fTransmission9km'] - laser.Transmission12km = laser_info_runh['MReportLaser.fTransmission12km'] - laser.Zenith = laser_info_runh['MReportLaser.fZenith']* u.deg - laser.Azimuth = laser_info_runh['MReportLaser.fAzimuth']* u.deg - laser.FullOverlap = laser_info_runh['MReportLaser.fFullOverlap'] - laser.EndGroundLayer = laser_info_runh['MReportLaser.fEndGroundLayer'] - laser.GroundLayerTrans = laser_info_runh['MReportLaser.fGroundLayerTrans'] - laser.Klett_k = laser_info_runh['MReportLaser.fKlett_k'] - laser.PheCounts = laser_info_runh['MReportLaser.fPheCounts'] - laser.Calimaness = laser_info_runh['MReportLaser.fCalimaness'] - laser.CloudLayerAlt = laser_info_runh['MReportLaser.fCloudLayerAlt'] - laser.CloudLayerDens = laser_info_runh['MReportLaser.fCloudLayerDens'] - laser.Offset = laser_info_runh['MReportLaser.fOffset'] - laser.Offset_Calculated = laser_info_runh['MReportLaser.fOffset_Calculated'] - laser.Offset_Fitted = laser_info_runh['MReportLaser.fOffset_Fitted'] - laser.Offset2 = laser_info_runh['MReportLaser.fOffset2'] - laser.Offset3 = laser_info_runh['MReportLaser.fOffset3'] - laser.Background1 = laser_info_runh['MReportLaser.fBackground1'] - laser.Background2 = laser_info_runh['MReportLaser.fBackground2'] - laser.BackgroundErr1 = laser_info_runh['MReportLaser.fBackgroundErr1'] - laser.BackgroundErr2 = laser_info_runh['MReportLaser.fBackgroundErr2'] - laser.RangeMax = laser_info_runh['MReportLaser.fRangeMax'] - laser.RangeMax_Clouds = laser_info_runh['MReportLaser.fRangeMax_Clouds'] - laser.ErrorCode = laser_info_runh['MReportLaser.fErrorCode'] - laser.ScaleHeight_fit = laser_info_runh['MReportLaser.fScaleHeight_fit'] - laser.Alpha_fit = laser_info_runh['MReportLaser.fChi2Alpha_fit'] - laser.Chi2Alpha_fit = laser_info_runh['MReportLaser.fChi2Alpha_fit'] - laser.Alpha_firstY = laser_info_runh['MReportLaser.fAlpha_firstY'] - laser.Alpha_Junge = laser_info_runh['MReportLaser.fAlpha_Junge'] - laser.PBLHeight = laser_info_runh['MReportLaser.fPBLHeight'] - laser.Chi2Full_fit = laser_info_runh['MReportLaser.fChi2Full_fit'] - laser.HWSwitchMaxOffset = laser_info_runh['MReportLaser.fHWSwitchMaxOffset'] - laser.NCollapse = laser_info_runh['MReportLaser.fNCollapse'] - laser.Shots = laser_info_runh['MReportLaser.fShots'] - laser.T0Shift = laser_info_runh['MReportLaser.fT0Shift'] - laser.Interval_0 = laser_info_runh['MReportLaser.fInterval_0'] - laser.RCS_min_perfect = laser_info_runh['MReportLaser.fRCS_min_perfect'] - laser.RCS_min_clouds = laser_info_runh['MReportLaser.fRCS_min_clouds'] - laser.RCS_min_mol = laser_info_runh['MReportLaser.fRCS_min_mol'] - laser.LIDAR_ratio = laser_info_runh['MReportLaser.fLIDAR_ratio'] - laser.LIDAR_ratio_Cloud = laser_info_runh['MReportLaser.fLIDAR_ratio_Cloud'] - laser.LIDAR_ratio_Junge = laser_info_runh['MReportLaser.fLIDAR_ratio_Junge'] - - millisec_seconds = millisec_values * 1e-3 - combined_mjd_value = mjd_values + millisec_seconds / 86400 - laser.MJD = combined_mjd_value - unique_reports[unique_key] = laser + unique_reports[unique_key] = [] + laser = ReportLaserContainer() + laser.UniqueID = laser_info_runh['MReportLaser.MReport.fUniqueID'] + laser.Bits = laser_info_runh['MReportLaser.MReport.fBits'] + #laser.BadReport = laser_info_runh['MReportLaser.MReport.fBadReport'] # to be improved + #laser.State = laser_info_runh['MReportLaser.MReport.fState'] # to be improved + laser.IsOffsetCorrection = laser_info_runh['MReportLaser.fIsOffsetCorrection'] + laser.IsOffsetFitted = laser_info_runh['MReportLaser.fIsOffsetFitted'] + laser.IsBGCorrection = laser_info_runh['MReportLaser.fIsBGCorrection'] + laser.IsT0ShiftFitted = laser_info_runh['MReportLaser.fIsT0ShiftFitted'] + laser.IsUseGDAS = laser_info_runh['MReportLaser.fIsUseGDAS'] + laser.IsUpwardMoving = laser_info_runh['MReportLaser.fIsUpwardMoving'] + laser.OverShoot = int(laser_info_runh['MReportLaser.fOverShoot']) + laser.UnderShoot = int(laser_info_runh['MReportLaser.fUnderShoot']) + laser.BGSamples = int(laser_info_runh['MReportLaser.fBGSamples']) + laser.Transmission3km = laser_info_runh['MReportLaser.fTransmission3km'] + laser.Transmission6km = laser_info_runh['MReportLaser.fTransmission6km'] + laser.Transmission9km = laser_info_runh['MReportLaser.fTransmission9km'] + laser.Transmission12km = laser_info_runh['MReportLaser.fTransmission12km'] + laser.Zenith = laser_info_runh['MReportLaser.fZenith']* u.deg + laser.Azimuth = laser_info_runh['MReportLaser.fAzimuth']* u.deg + laser.FullOverlap = laser_info_runh['MReportLaser.fFullOverlap'] + laser.EndGroundLayer = laser_info_runh['MReportLaser.fEndGroundLayer'] + laser.GroundLayerTrans = laser_info_runh['MReportLaser.fGroundLayerTrans'] + laser.Klett_k = laser_info_runh['MReportLaser.fKlett_k'] + laser.PheCounts = laser_info_runh['MReportLaser.fPheCounts'] + laser.Calimaness = laser_info_runh['MReportLaser.fCalimaness'] + laser.CloudLayerAlt = laser_info_runh['MReportLaser.fCloudLayerAlt'] + laser.CloudLayerDens = laser_info_runh['MReportLaser.fCloudLayerDens'] + laser.Offset = laser_info_runh['MReportLaser.fOffset'] + laser.Offset_Calculated = laser_info_runh['MReportLaser.fOffset_Calculated'] + laser.Offset_Fitted = laser_info_runh['MReportLaser.fOffset_Fitted'] + laser.Offset2 = laser_info_runh['MReportLaser.fOffset2'] + laser.Offset3 = laser_info_runh['MReportLaser.fOffset3'] + laser.Background1 = laser_info_runh['MReportLaser.fBackground1'] + laser.Background2 = laser_info_runh['MReportLaser.fBackground2'] + laser.BackgroundErr1 = laser_info_runh['MReportLaser.fBackgroundErr1'] + laser.BackgroundErr2 = laser_info_runh['MReportLaser.fBackgroundErr2'] + laser.RangeMax = laser_info_runh['MReportLaser.fRangeMax'] + laser.RangeMax_Clouds = laser_info_runh['MReportLaser.fRangeMax_Clouds'] + laser.ErrorCode = laser_info_runh['MReportLaser.fErrorCode'] + laser.ScaleHeight_fit = laser_info_runh['MReportLaser.fScaleHeight_fit'] + laser.Alpha_fit = laser_info_runh['MReportLaser.fChi2Alpha_fit'] + laser.Chi2Alpha_fit = laser_info_runh['MReportLaser.fChi2Alpha_fit'] + laser.Alpha_firstY = laser_info_runh['MReportLaser.fAlpha_firstY'] + laser.Alpha_Junge = laser_info_runh['MReportLaser.fAlpha_Junge'] + laser.PBLHeight = laser_info_runh['MReportLaser.fPBLHeight'] + laser.Chi2Full_fit = laser_info_runh['MReportLaser.fChi2Full_fit'] + laser.HWSwitchMaxOffset = laser_info_runh['MReportLaser.fHWSwitchMaxOffset'] + laser.NCollapse = laser_info_runh['MReportLaser.fNCollapse'] + laser.Shots = laser_info_runh['MReportLaser.fShots'] + laser.T0Shift = laser_info_runh['MReportLaser.fT0Shift'] + laser.Interval_0 = laser_info_runh['MReportLaser.fInterval_0'] + laser.RCS_min_perfect = laser_info_runh['MReportLaser.fRCS_min_perfect'] + laser.RCS_min_clouds = laser_info_runh['MReportLaser.fRCS_min_clouds'] + laser.RCS_min_mol = laser_info_runh['MReportLaser.fRCS_min_mol'] + laser.LIDAR_ratio = laser_info_runh['MReportLaser.fLIDAR_ratio'] + laser.LIDAR_ratio_Cloud = laser_info_runh['MReportLaser.fLIDAR_ratio_Cloud'] + laser.LIDAR_ratio_Junge = laser_info_runh['MReportLaser.fLIDAR_ratio_Junge'] + + millisec_seconds = millisec_values * 1e-3 + combined_mjd_value = mjd_values + millisec_seconds / 86400 + laser.MJD = combined_mjd_value + unique_reports[unique_key].append(laser) except KeyError as e: print(f"Required key not found in the file {rootf}: {e}") continue - return laser + return unique_reports def parse_simulation_header(self): """ From 1864e765f4c05011266432ce968d0363b887a75f Mon Sep 17 00:00:00 2001 From: nzywucka Date: Wed, 27 Mar 2024 11:06:29 +0100 Subject: [PATCH 25/38] test_lidar_parameters was removed --- ctapipe_io_magic/tests/test_magic_event_source.py | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/ctapipe_io_magic/tests/test_magic_event_source.py b/ctapipe_io_magic/tests/test_magic_event_source.py index 1b46c86..f57988f 100644 --- a/ctapipe_io_magic/tests/test_magic_event_source.py +++ b/ctapipe_io_magic/tests/test_magic_event_source.py @@ -601,16 +601,3 @@ def test_broken_subruns_missing_arrays(): # / "20210314_M1_05095172.002_Y_CrabNebula-W0.40+035.root" #) -@pytest.mark.parametrize("dataset", test_calibrated_all) -def test_lidar_parameters(dataset): - from ctapipe_io_magic import MAGICEventSource, ReportLaserContainer - - dataset = (test_calibrated_real_dir/"20210314_M1_05095172.002_Y_CrabNebula-W0.40+035.root") - - with MAGICEventSource(input_url=dataset) as source: - assert source.laser.Transmission12km == pytest.approx(0.89) - assert source.laser.Transmission3km == pytest.approx(0.96) - assert source.laser.Transmission6km == pytest.approx(0.93) - assert source.laser.Transmission9km == pytest.approx(0.89) - - From 61f4a7506f5427e12c960038131000d71c73cb1d Mon Sep 17 00:00:00 2001 From: nzywucka Date: Wed, 27 Mar 2024 11:14:24 +0100 Subject: [PATCH 26/38] Implementation of some tests of reading lidar parameters --- ctapipe_io_magic/tests/test_lidar.py | 49 ++++++++++++++++++++++++++++ 1 file changed, 49 insertions(+) create mode 100644 ctapipe_io_magic/tests/test_lidar.py diff --git a/ctapipe_io_magic/tests/test_lidar.py b/ctapipe_io_magic/tests/test_lidar.py new file mode 100644 index 0000000..57cfe01 --- /dev/null +++ b/ctapipe_io_magic/tests/test_lidar.py @@ -0,0 +1,49 @@ +import copy +import os +import numpy as np +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", +] + +data_dict = dict() + +data_dict["20210314_M1_05095172.002_Y_CrabNebula-W0.40+035.root"] = dict() + +data_dict["20210314_M1_05095172.002_Y_CrabNebula-W0.40+035.root"]["n_events_tot"] = 500 +data_dict["20210314_M1_05095172.002_Y_CrabNebula-W0.40+035.root"]["n_events_stereo"] = 452 +data_dict["20210314_M1_05095172.002_Y_CrabNebula-W0.40+035.root"]["n_events_pedestal"] = 48 +data_dict["20210314_M1_05095172.002_Y_CrabNebula-W0.40+035.root"]["n_events_mc_mono"] = 0 + +def test_sample_report_laser(sample_report_laser): + # Write your test code here + assert sample_report_laser is not None + +# Define a fixture to create an instance of ReportLaserContainer for testing +@pytest.fixture +def sample_report_laser(): + from ctapipe_io_magic import MAGICEventSource, ReportLaserContainer + dataset = (test_calibrated_real_dir/"20210314_M1_05095172.002_Y_CrabNebula-W0.40+035.root") + with MAGICEventSource(input_url=dataset) as source: + print(ReportLaserContainer()) + return ReportLaserContainer() + +@pytest.mark.parametrize("dataset", test_calibrated_real) +def test_lidar_parameters(dataset): + from ctapipe_io_magic import MAGICEventSource + + dataset = (test_calibrated_real_dir/"20210314_M1_05095172.002_Y_CrabNebula-W0.40+035.root") + with MAGICEventSource(input_url=dataset) as source: + assert source.laser.Transmission12km == pytest.approx(0.89) #float + assert source.laser.IsUseGDAS == pytest.approx([False]) # bool + azimuth_degrees = source.laser.Azimuth.value # Convert Angle to degrees + assert azimuth_degrees == pytest.approx(276.63) # Angle + flattened_array = np.concatenate(source.laser.PheCounts) + phe_counts = [int(str(value).replace('\n', '').strip()) for value in flattened_array] + assert 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] # array + assert source.laser.Bits == (50331656) # Any From 94efd5e5dea93b5124345d8714eb86cd1b2cd301 Mon Sep 17 00:00:00 2001 From: nzywucka Date: Wed, 27 Mar 2024 11:46:56 +0100 Subject: [PATCH 27/38] test_lidar_parameters test has been improved. --- ctapipe_io_magic/tests/test_lidar.py | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/ctapipe_io_magic/tests/test_lidar.py b/ctapipe_io_magic/tests/test_lidar.py index 57cfe01..c16e0d9 100644 --- a/ctapipe_io_magic/tests/test_lidar.py +++ b/ctapipe_io_magic/tests/test_lidar.py @@ -36,14 +36,18 @@ def sample_report_laser(): @pytest.mark.parametrize("dataset", test_calibrated_real) def test_lidar_parameters(dataset): from ctapipe_io_magic import MAGICEventSource - - dataset = (test_calibrated_real_dir/"20210314_M1_05095172.002_Y_CrabNebula-W0.40+035.root") + + # Ensure dataset is properly passed from parametrize with MAGICEventSource(input_url=dataset) as source: - assert source.laser.Transmission12km == pytest.approx(0.89) #float - assert source.laser.IsUseGDAS == pytest.approx([False]) # bool - azimuth_degrees = source.laser.Azimuth.value # Convert Angle to degrees - assert azimuth_degrees == pytest.approx(276.63) # Angle - flattened_array = np.concatenate(source.laser.PheCounts) - phe_counts = [int(str(value).replace('\n', '').strip()) for value in flattened_array] - assert 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] # array - assert source.laser.Bits == (50331656) # Any + laser = source.laser + for key, report_list in laser.items(): + for report in report_list: + assert report.Transmission12km == pytest.approx(0.89, abs=8.9e-07) + assert report.IsUseGDAS == pytest.approx([False]) # bool + azimuth_degrees = report.Azimuth.value # Convert Angle to degrees + assert azimuth_degrees == pytest.approx(276.63) # Angle + flattened_array = np.concatenate(report.PheCounts) + phe_counts = [int(str(value).replace('\n', '').strip()) for value in flattened_array] + assert 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.Bits == (50331656) # Any + From d46221fd5638cf2bc29b49f610c28b88fa8c81c1 Mon Sep 17 00:00:00 2001 From: nzywucka Date: Wed, 27 Mar 2024 12:09:54 +0100 Subject: [PATCH 28/38] Unnecessary comments and whitespaces have been removed. --- ctapipe_io_magic/tests/test_lidar.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/ctapipe_io_magic/tests/test_lidar.py b/ctapipe_io_magic/tests/test_lidar.py index c16e0d9..fe7268b 100644 --- a/ctapipe_io_magic/tests/test_lidar.py +++ b/ctapipe_io_magic/tests/test_lidar.py @@ -1,4 +1,3 @@ -import copy import os import numpy as np from pathlib import Path @@ -21,10 +20,8 @@ data_dict["20210314_M1_05095172.002_Y_CrabNebula-W0.40+035.root"]["n_events_mc_mono"] = 0 def test_sample_report_laser(sample_report_laser): - # Write your test code here assert sample_report_laser is not None -# Define a fixture to create an instance of ReportLaserContainer for testing @pytest.fixture def sample_report_laser(): from ctapipe_io_magic import MAGICEventSource, ReportLaserContainer @@ -36,8 +33,6 @@ def sample_report_laser(): @pytest.mark.parametrize("dataset", test_calibrated_real) def test_lidar_parameters(dataset): from ctapipe_io_magic import MAGICEventSource - - # Ensure dataset is properly passed from parametrize with MAGICEventSource(input_url=dataset) as source: laser = source.laser for key, report_list in laser.items(): From 338a7a9df653d78615adf5d10894749ff198df2b Mon Sep 17 00:00:00 2001 From: nzywucka Date: Wed, 27 Mar 2024 12:27:58 +0100 Subject: [PATCH 29/38] Implementation of some corrections of lidar parameters. --- ctapipe_io_magic/laser_params_corr_test.py | 215 +++++++++++++++++++++ 1 file changed, 215 insertions(+) create mode 100644 ctapipe_io_magic/laser_params_corr_test.py diff --git a/ctapipe_io_magic/laser_params_corr_test.py b/ctapipe_io_magic/laser_params_corr_test.py new file mode 100644 index 0000000..9d5d252 --- /dev/null +++ b/ctapipe_io_magic/laser_params_corr_test.py @@ -0,0 +1,215 @@ +import logging +import datetime +from pathlib import Path +import os +import numpy as np +import astropy.units as u +from ctapipe_io_magic import MAGICEventSource +import zd_az_corr +import params +from switch_times import switch_times + +logger = logging.getLogger(__name__) +logger.addHandler(logging.StreamHandler()) +logger.setLevel(logging.INFO) + +class ReportLaser: + def magic_reports_reading(self, input_file, process_run=False): + logger.info(f"\nInput file: {input_file}") + event_source = MAGICEventSource(input_file, process_run=process_run) + is_simulation = event_source.is_simulation + logger.info(f"\nIs simulation: {is_simulation}") + + obs_id = event_source.obs_ids[0] + tel_id = event_source.telescope + + logger.info(f"\nObservation ID: {obs_id}") + logger.info(f"Telescope ID: {tel_id}") + + laser = event_source.laser + + for key, report_list in laser.items(): + for report in report_list: + logger.info(f"Accessing parameters for each ReportLaserContainer object") + + return report + + def get_c0(self, time, zenith, range_max): + C_0 = params.gkC_0Period1 # absolute calibration at beginning of 2013 + if params.stime21 > 0: + fFullOverlap = params.gkFullOverlap5 # overlap has moved, after installation of new DAQ, seems that saturation occurs earlier now + if params.stime22 > 0: + fFullOverlap = 400. # overlap has moved, after installation of new DAQ, seems that saturation occurs earlier now + if params.stime232 > 0: + fFullOverlap = params.gkFullOverlap6 # overlap has moved, after installation of new DAQ, seems that saturation occurs earlier now + if params.stime233 > 0: + fFullOverlap = params.gkFullOverlap5 # overlap has moved, after installation of new DAQ, seems that saturation occurs earlier now + if params.stime234 > 0: + fFullOverlap = params.gkFullOverlap7 # overlap has moved, after installation of new DAQ, seems that saturation occurs earlier now + if params.stime26 > 0: + fFullOverlap = params.gkFullOverlap5 # overlap was corrected again, valid for only one night + if params.stime261 > 0: + fFullOverlap = params.gkFullOverlap7 # overlap got worse, probably screw loose + if params.stime27 > 0: + fFullOverlap = 1500. # overlap got worse, probably screw loose + if params.stime285 > 0: + fFullOverlap = 1000. # overlap got worse, probably screw loose + range_max_clouds = 23000. + + zenith_radians = np.deg2rad(zenith) + + if (range_max - 100.) * np.cos(zenith_radians) < range_max_clouds: + print("Reducing maximum range to search for clouds from", range_max_clouds, "to", + (range_max - 100.) * np.cos(zenith_radians), "Zd=", zenith) + range_max_clouds = (range_max - 100.) * np.cos(zenith_radians) + return C_0 + + def apply_time_dependent_corrections(self, datime, c0, coszd, case_index=None): + stimes = [(datime - switch_time).total_seconds() for switch_time in switch_times] + Alphafit_corr = 0.0 + + def apply_zenith_correction(self, c0, coszd): + gkC_0ZenithCorr = 0.038 + zd_corr = np.log(1 - gkC_0ZenithCorr * (1 - coszd)) + c0 += 2 * zd_corr + return True + + def apply_azimuth_correction(self, c0, zenith, azimuth, zd_az_corr): + zenith_threshold = zd_az_corr.correction_rules[0].get('zenith_threshold', None) + azimuth_threshold = zd_az_corr.correction_rules[0].get('azimuth_threshold', None) + zenith_threshold_min = zd_az_corr.correction_rules[0].get('zenith_threshold_min', None) + zenith_threshold_max = zd_az_corr.correction_rules[0].get('zenith_threshold_max', None) + azimuth_threshold_min = zd_az_corr.correction_rules[0].get('azimuth_threshold_min', None) + azimuth_threshold_max = zd_az_corr.correction_rules[0].get('azimuth_threshold_max', None) + + if zenith_threshold is not None and azimuth_threshold is not None: + if (zenith < zenith_threshold * u.deg and + azimuth_threshold[0] * u.deg < azimuth < azimuth_threshold[1] * u.deg): + c0 -= 2 * zd_az_corr.correction_rules[0]['correction_factor'] + return True + elif zenith_threshold_min is not None and zenith_threshold_max is not None: + if (zenith_threshold_min * u.deg < zenith < zenith_threshold_max * u.deg and + azimuth_threshold_min * u.deg < azimuth < azimuth_threshold_max * u.deg): + c0 -= 2 * zd_az_corr.correction_rules[0]['correction_factor'] + return True + return False + + def get_bin_width(self, ifadc, bgsamples, ncollapse): + #fBGSamples = 0 + if bgsamples <= params.gkBGSamplesV1to4: + return 1.0 + if ifadc < 128: + return 0.5 + if ifadc < 256: + return 1.0 + if ifadc < 384: + return 2.0 + if ifadc < 448: + return 4.0 + return 4.0 * ncollapse + + def get_offset_r(self, ifadc, bgsamples, ncollapse): + #fBGSamples = 0 + if bgsamples <= params.gkBGSamplesV1to4: + return 0.0 + if ifadc < 128: + return - bgsamples * params.gkIntegrationWindow * self.get_bin_width(ifadc, bgsamples-1, ncollapse) + r = (128 - bgsamples) * params.gkIntegrationWindow * self.get_bin_width(ifadc, bgsamples-1, ncollapse) + if ifadc < 256: + return r + r += 128 * params.gkIntegrationWindow + if ifadc < 384: + return r + r += 128 * params.gkIntegrationWindow * 2 + if ifadc < 448: + return r + r += 64 * params.gkIntegrationWindow * 4 + return r + + def get_fadcr(self, i, bgsamples): + #fBGSamples = 0 + if bgsamples <= params.gkBGSamplesV1to4: + return i + ifadc = i + bgsamples + if ifadc < 448: + return ifadc % 128 + return ifadc % 64 + + def get_collapsed_signal(self, transmission6km, background1, phecounts, bg, bgvar): + signalsamples = 0 + sig = np.zeros(signalsamples) + if transmission6km < -998.: + for i in range(466, signalsamples): + background1 += np.power((phecounts[i] - phecounts[i - 1]) * 1E6 / np.power(params.gkIntegrationWindow * i, 2), 2) + background1 /= 30.0 + bg = background1 + for i in range(signalsamples): + sig[i] = pheounts[i] * 1E6 / np.power(params.gkIntegrationWindow * i, 2) + background1 + return sig + + def get_range(self, i, bgsamples, ncollapse, t0shift): + ifadc = i + bgsamples + return (self.get_offset_r(ifadc, bgsamples, ncollapse) + params.gkIntegrationWindow * (self.get_fadcr(i, bgsamples) + 0.5) * self.get_bin_width(ifadc, bgsamples, ncollapse) + + t0shift * params.gkIntegrationWindow * self.get_bin_width(ifadc, bgsamples, ncollapse)) + + def get_signal(self, i, bgsamples, phecounts, ncollapse, bg, hwcorr): + ifadc = i + bgsamples + return ((phecounts - self.get_bin_width(ifadc, bgsamples, ncollapse) * bg) * hwcorr + / self.get_bin_width(ifadc, bgsamples, ncollapse)) + +def main(): + test_data = Path(os.getenv("MAGIC_TEST_DATA", "test_data")).absolute() + test_calibrated_real_dir = test_data / "real/calibrated" + input_file = (test_calibrated_real_dir/"20210314_M1_05095172.002_Y_CrabNebula-W0.40+035.root") + report = ReportLaser() + laser_params = report.magic_reports_reading(input_file) + + time = datetime.datetime(2020, 7, 26, 3, 10, 0) + zenith = laser_params['Zenith'] + azimuth = laser_params['Azimuth'] + range_max = laser_params['RangeMax'] + phecounts = laser_params['PheCounts'] + bgsamples = laser_params['BGSamples'] + ncollapse = laser_params['NCollapse'] + transmission6km = laser_params['Transmission6km'] + background1 = laser_params['Background1'] + t0shift = laser_params['T0Shift'] + coszd = np.cos(np.deg2rad(zenith)) + ifadc = 200 + bg = 0.0 + bgvar = 0.0 + + C_0 = report.get_c0(time, zenith, range_max) + print("C_0 value:", C_0) + + case_index = 15 + report.apply_time_dependent_corrections(time, C_0, zenith, case_index) + + c0 = report.apply_zenith_correction(C_0, coszd) + print("Zenith corrected c0:", c0) + + c0 = report.apply_azimuth_correction(C_0, zenith, azimuth, zd_az_corr) + + print("Zenith and Azimuth corrected c0:", c0) + + bin_width = report.get_bin_width(ifadc, bgsamples, ncollapse) + print("Bin Width:", bin_width) + + offset_r = report.get_offset_r(ifadc, bgsamples, ncollapse) + print("Offset R:", offset_r) + + fadc_r = report.get_fadcr(150, bgsamples) + print("FADCR:", fadc_r) + + collapsed_signal = report.get_collapsed_signal(transmission6km, background1, phecounts, bg, bgvar) + + print("Collapsed Signal:", collapsed_signal) + print("Background:", bg) + print("Background Variance:", bgvar) + + range_value = report.get_range(0, bgsamples, ncollapse, t0shift) + print("Range:", range_value) + +if __name__ == "__main__": + main() + From 2ef3ecd3cc229660569060f6e90005060c896ec3 Mon Sep 17 00:00:00 2001 From: nzywucka Date: Thu, 28 Mar 2024 11:18:45 +0100 Subject: [PATCH 30/38] Improvements to read external files have been added --- ctapipe_io_magic/laser_params_corr_test.py | 300 +++++++++++++++++---- 1 file changed, 248 insertions(+), 52 deletions(-) diff --git a/ctapipe_io_magic/laser_params_corr_test.py b/ctapipe_io_magic/laser_params_corr_test.py index 9d5d252..9c574d8 100644 --- a/ctapipe_io_magic/laser_params_corr_test.py +++ b/ctapipe_io_magic/laser_params_corr_test.py @@ -5,9 +5,6 @@ import numpy as np import astropy.units as u from ctapipe_io_magic import MAGICEventSource -import zd_az_corr -import params -from switch_times import switch_times logger = logging.getLogger(__name__) logger.addHandler(logging.StreamHandler()) @@ -35,24 +32,70 @@ def magic_reports_reading(self, input_file, process_run=False): return report def get_c0(self, time, zenith, range_max): - C_0 = params.gkC_0Period1 # absolute calibration at beginning of 2013 - if params.stime21 > 0: - fFullOverlap = params.gkFullOverlap5 # overlap has moved, after installation of new DAQ, seems that saturation occurs earlier now - if params.stime22 > 0: + filename = '/home/zywuckan/ctapipe_io_magic/ctapipe_io_magic/params.txt' + + gkC_0Period1 = None + gkFullOverlap5 = None + gkFullOverlap6 = None + gkFullOverlap7 = None + gkIntegrationWindow = None + gkHWSwitchV1to4 = None + + with open(filename, 'r') as file: + for line in file: + if line.strip(): + variable, value = line.split("=") + variable = variable.strip() + value = value.split("#")[0].strip() + if variable == "stime21": + stime21 = int(value) + elif variable == "stime22": + stime22 = int(value) + elif variable == "stime26": + stime26 = int(value) + elif variable == "stime261": + stime261 = int(value) + elif variable == "stime232": + stime232 = int(value) + elif variable == "stime233": + stime233 = int(value) + elif variable == "stime234": + stime234 = int(value) + elif variable == "stime27": + stime27 = int(value) + elif variable == "stime285": + stime285 = int(value) + elif variable == "gkC_0Period1": + gkC_0Period1 = float(value) + elif variable == "gkFullOverlap5": + gkFullOverlap5 = float(value) + elif variable == "gkFullOverlap6": + gkFullOverlap6 = float(value) + elif variable == "gkFullOverlap7": + gkFullOverlap7 = float(value) + elif variable == "gkIntegrationWindow": + gkIntegrationWindow = float(value) + elif variable == "gkHWSwitchV1to4": + gkHWSwitchV1to4 = int(value) + + C_0 = gkC_0Period1 # absolute calibration at beginning of 2013 + if stime21 > 0: + fFullOverlap = gkFullOverlap5 # overlap has moved, after installation of new DAQ, seems that saturation occurs earlier now + if stime22 > 0: fFullOverlap = 400. # overlap has moved, after installation of new DAQ, seems that saturation occurs earlier now - if params.stime232 > 0: - fFullOverlap = params.gkFullOverlap6 # overlap has moved, after installation of new DAQ, seems that saturation occurs earlier now - if params.stime233 > 0: - fFullOverlap = params.gkFullOverlap5 # overlap has moved, after installation of new DAQ, seems that saturation occurs earlier now - if params.stime234 > 0: - fFullOverlap = params.gkFullOverlap7 # overlap has moved, after installation of new DAQ, seems that saturation occurs earlier now - if params.stime26 > 0: - fFullOverlap = params.gkFullOverlap5 # overlap was corrected again, valid for only one night - if params.stime261 > 0: - fFullOverlap = params.gkFullOverlap7 # overlap got worse, probably screw loose - if params.stime27 > 0: + if stime232 > 0: + fFullOverlap = gkFullOverlap6 # overlap has moved, after installation of new DAQ, seems that saturation occurs earlier now + if stime233 > 0: + fFullOverlap = gkFullOverlap5 # overlap has moved, after installation of new DAQ, seems that saturation occurs earlier now + if stime234 > 0: + fFullOverlap = gkFullOverlap7 # overlap has moved, after installation of new DAQ, seems that saturation occurs earlier now + if stime26 > 0: + fFullOverlap = gkFullOverlap5 # overlap was corrected again, valid for only one night + if stime261 > 0: + fFullOverlap = gkFullOverlap7 # overlap got worse, probably screw loose + if stime27 > 0: fFullOverlap = 1500. # overlap got worse, probably screw loose - if params.stime285 > 0: + if stime285 > 0: fFullOverlap = 1000. # overlap got worse, probably screw loose range_max_clouds = 23000. @@ -65,6 +108,18 @@ def get_c0(self, time, zenith, range_max): return C_0 def apply_time_dependent_corrections(self, datime, c0, coszd, case_index=None): + switch_times = {} + filename = '/home/zywuckan/ctapipe_io_magic/ctapipe_io_magic/switch_times.txt' + with open(filename, 'r') as file: + switch_time = None + for line in file: + if line.startswith("Switch Time:"): + switch_time = datetime.datetime.strptime(line.split(": ")[1].strip(), "%Y-%m-%d %H:%M:%S") + elif line.startswith("Correction Parameters:"): + correction_params = eval(line.split(": ")[1].strip()) + if switch_time: + switch_times[switch_time] = correction_params + stimes = [(datime - switch_time).total_seconds() for switch_time in switch_times] Alphafit_corr = 0.0 @@ -74,29 +129,76 @@ def apply_zenith_correction(self, c0, coszd): c0 += 2 * zd_corr return True - def apply_azimuth_correction(self, c0, zenith, azimuth, zd_az_corr): - zenith_threshold = zd_az_corr.correction_rules[0].get('zenith_threshold', None) - azimuth_threshold = zd_az_corr.correction_rules[0].get('azimuth_threshold', None) - zenith_threshold_min = zd_az_corr.correction_rules[0].get('zenith_threshold_min', None) - zenith_threshold_max = zd_az_corr.correction_rules[0].get('zenith_threshold_max', None) - azimuth_threshold_min = zd_az_corr.correction_rules[0].get('azimuth_threshold_min', None) - azimuth_threshold_max = zd_az_corr.correction_rules[0].get('azimuth_threshold_max', None) - - if zenith_threshold is not None and azimuth_threshold is not None: - if (zenith < zenith_threshold * u.deg and - azimuth_threshold[0] * u.deg < azimuth < azimuth_threshold[1] * u.deg): - c0 -= 2 * zd_az_corr.correction_rules[0]['correction_factor'] - return True - elif zenith_threshold_min is not None and zenith_threshold_max is not None: - if (zenith_threshold_min * u.deg < zenith < zenith_threshold_max * u.deg and - azimuth_threshold_min * u.deg < azimuth < azimuth_threshold_max * u.deg): - c0 -= 2 * zd_az_corr.correction_rules[0]['correction_factor'] + def apply_azimuth_correction(self, c0, zenith, azimuth): + filename = '/home/zywuckan/ctapipe_io_magic/ctapipe_io_magic/zd_az_corr.txt' + rules = [] + with open(filename, 'r') as file: + rule = {} + for line in file: + if line.startswith("Rule"): + if rule: + rules.append(rule) + rule = {} + elif line.startswith("Zenith Threshold:"): + rule['zenith_threshold'] = float(line.split(": ")[1].strip()) + elif line.startswith("Azimuth Threshold:"): + rule['azimuth_threshold'] = eval(line.split(": ")[1].strip()) + elif line.startswith("Zenith Threshold Min:"): + rule['zenith_threshold_min'] = float(line.split(": ")[1].strip()) + elif line.startswith("Zenith Threshold Max:"): + rule['zenith_threshold_max'] = float(line.split(": ")[1].strip()) + elif line.startswith("Azimuth Threshold Min:"): + rule['azimuth_threshold_min'] = float(line.split(": ")[1].strip()) + elif line.startswith("Azimuth Threshold Max:"): + rule['azimuth_threshold_max'] = float(line.split(": ")[1].strip()) + elif line.startswith("Correction Factor:"): + rule['correction_factor'] = float(line.split(": ")[1].strip()) + if rule: + rules.append(rule) + + for rule in rules: + zenith_threshold = rule.get('zenith_threshold', None) + azimuth_threshold = rule.get('azimuth_threshold', None) + zenith_threshold_min = rule.get('zenith_threshold_min', None) + zenith_threshold_max = rule.get('zenith_threshold_max', None) + azimuth_threshold_min = rule.get('azimuth_threshold_min', None) + azimuth_threshold_max = rule.get('azimuth_threshold_max', None) + + if zenith_threshold is not None and azimuth_threshold is not None: + if (zenith < zenith_threshold * u.deg and + azimuth_threshold[0] * u.deg < azimuth < azimuth_threshold[1] * u.deg): + c0 -= 2 * rule['correction_factor'] return True + elif zenith_threshold_min is not None and zenith_threshold_max is not None: + if (zenith_threshold_min * u.deg < zenith < zenith_threshold_max * u.deg and + azimuth_threshold_min * u.deg < azimuth < azimuth_threshold_max * u.deg): + c0 -= 2 * rule['correction_factor'] + return True return False def get_bin_width(self, ifadc, bgsamples, ncollapse): - #fBGSamples = 0 - if bgsamples <= params.gkBGSamplesV1to4: + filename = '/home/zywuckan/ctapipe_io_magic/ctapipe_io_magic/params.txt' + gkBGSamplesV1to4 = None + with open(filename, 'r') as file: + for line in file: + line = line.strip() + if line and not line.startswith("#"): # Skip empty lines and comments + variable, value = line.split("=") + variable = variable.strip() + value = value.split("#")[0].strip() # Remove comments + if variable == "gkBGSamplesV1to4": + try: + gkBGSamplesV1to4 = int(value) + except ValueError: + raise ValueError("Invalid value for gkBGSamplesV1to4 in params.txt") + + if gkBGSamplesV1to4 is None: + raise ValueError("Variable gkBGSamplesV1to4 not found in params.txt") + + if not isinstance(bgsamples, int): + raise ValueError("Invalid type for bgsamples") + + if bgsamples <= gkBGSamplesV1to4: return 1.0 if ifadc < 128: return 0.5 @@ -109,26 +211,75 @@ def get_bin_width(self, ifadc, bgsamples, ncollapse): return 4.0 * ncollapse def get_offset_r(self, ifadc, bgsamples, ncollapse): - #fBGSamples = 0 - if bgsamples <= params.gkBGSamplesV1to4: + filename = '/home/zywuckan/ctapipe_io_magic/ctapipe_io_magic/params.txt' + gkBGSamplesV1to4 = None + gkIntegrationWindow = None + with open(filename, 'r') as file: + for line in file: + exec(line.strip()) + if line.strip(): # Check if line is not empty + variable, value = line.split("=") + variable = variable.strip() + value = value.split("#")[0].strip() # Remove comments + if variable == "gkBGSamplesV1to4": + try: + gkBGSamplesV1to4 = int(value) + except ValueError: + raise ValueError("Invalid value for gkBGSamplesV1to4 in params.txt") + elif variable == "gkIntegrationWindow": # Fix indentation here + try: + gkIntegrationWindow = float(value) + except ValueError: + raise ValueError("Invalid value for gkIntegrationWindow in params.txt") + + if gkBGSamplesV1to4 is None: + raise ValueError("Variable gkBGSamplesV1to4 not found in params.txt") + + if gkIntegrationWindow is None: + raise ValueError("Variable gkIntegrationWindow not found in params.txt") + + if not isinstance(bgsamples, int): + raise ValueError("Invalid type for bgsamples") + + if bgsamples <= gkBGSamplesV1to4: return 0.0 if ifadc < 128: - return - bgsamples * params.gkIntegrationWindow * self.get_bin_width(ifadc, bgsamples-1, ncollapse) - r = (128 - bgsamples) * params.gkIntegrationWindow * self.get_bin_width(ifadc, bgsamples-1, ncollapse) + return - bgsamples * gkIntegrationWindow * self.get_bin_width(ifadc, bgsamples-1, ncollapse) + r = (128 - bgsamples) * gkIntegrationWindow * self.get_bin_width(ifadc, bgsamples-1, ncollapse) if ifadc < 256: return r - r += 128 * params.gkIntegrationWindow + r += 128 * gkIntegrationWindow if ifadc < 384: return r - r += 128 * params.gkIntegrationWindow * 2 + r += 128 * gkIntegrationWindow * 2 if ifadc < 448: return r - r += 64 * params.gkIntegrationWindow * 4 + r += 64 * gkIntegrationWindow * 4 return r def get_fadcr(self, i, bgsamples): - #fBGSamples = 0 - if bgsamples <= params.gkBGSamplesV1to4: + filename = '/home/zywuckan/ctapipe_io_magic/ctapipe_io_magic/params.txt' + gkBGSamplesV1to4 = None + with open(filename, 'r') as file: + for line in file: + exec(line.strip()) + if line.strip(): # Check if line is not empty + variable, value = line.split("=") + variable = variable.strip() + value = value.split("#")[0].strip() # Remove comments + if variable == "gkBGSamplesV1to4": + try: + gkBGSamplesV1to4 = int(value) + except ValueError: + raise ValueError("Invalid value for gkBGSamplesV1to4 in params.txt") + + if gkBGSamplesV1to4 is None: + raise ValueError("Variable gkBGSamplesV1to4 not found in params.txt") + + if not isinstance(bgsamples, int): + raise ValueError("Invalid type for bgsamples") + + if bgsamples <= gkBGSamplesV1to4: return i ifadc = i + bgsamples if ifadc < 448: @@ -136,21 +287,66 @@ def get_fadcr(self, i, bgsamples): return ifadc % 64 def get_collapsed_signal(self, transmission6km, background1, phecounts, bg, bgvar): + filename = '/home/zywuckan/ctapipe_io_magic/ctapipe_io_magic/params.txt' + gkIntegrationWindow = None + with open(filename, 'r') as file: + for line in file: + exec(line.strip()) + if line.strip(): # Check if line is not empty + variable, value = line.split("=") + variable = variable.strip() + value = value.split("#")[0].strip() # Remove comments + if variable == "gkIntegrationWindow": # Fix indentation here + try: + gkIntegrationWindow = float(value) + except ValueError: + raise ValueError("Invalid value for gkIntegrationWindow in params.txt") + + if gkIntegrationWindow is None: + raise ValueError("Variable gkIntegrationWindow not found in params.txt") + signalsamples = 0 sig = np.zeros(signalsamples) if transmission6km < -998.: for i in range(466, signalsamples): - background1 += np.power((phecounts[i] - phecounts[i - 1]) * 1E6 / np.power(params.gkIntegrationWindow * i, 2), 2) + background1 += np.power((phecounts[i] - phecounts[i - 1]) * 1E6 / np.power(gkIntegrationWindow * i, 2), 2) background1 /= 30.0 bg = background1 for i in range(signalsamples): - sig[i] = pheounts[i] * 1E6 / np.power(params.gkIntegrationWindow * i, 2) + background1 + sig[i] = pheounts[i] * 1E6 / np.power(gkIntegrationWindow * i, 2) + background1 return sig def get_range(self, i, bgsamples, ncollapse, t0shift): + filename = '/home/zywuckan/ctapipe_io_magic/ctapipe_io_magic/params.txt' + gkBGSamplesV1to4 = None + gkIntegrationWindow = None + with open(filename, 'r') as file: + for line in file: + exec(line.strip()) + if line.strip(): # Check if line is not empty + variable, value = line.split("=") + variable = variable.strip() + value = value.split("#")[0].strip() + if variable == "gkBGSamplesV1to4": + try: + gkBGSamplesV1to4 = int(value) + except ValueError: + raise ValueError("Invalid value for gkBGSamplesV1to4 in params.txt") + elif variable == "gkIntegrationWindow": # Fix indentation here + try: + gkIntegrationWindow = float(value) + except ValueError: + raise ValueError("Invalid value for gkIntegrationWindow in params.txt") + + if gkBGSamplesV1to4 is None: + raise ValueError("Variable gkBGSamplesV1to4 not found in params.txt") + + if gkIntegrationWindow is None: + raise ValueError("Variable gkIntegrationWindow not found in params.txt") + ifadc = i + bgsamples - return (self.get_offset_r(ifadc, bgsamples, ncollapse) + params.gkIntegrationWindow * (self.get_fadcr(i, bgsamples) + 0.5) * self.get_bin_width(ifadc, bgsamples, ncollapse) - + t0shift * params.gkIntegrationWindow * self.get_bin_width(ifadc, bgsamples, ncollapse)) + return (self.get_offset_r(ifadc, bgsamples, ncollapse) + gkIntegrationWindow * (self.get_fadcr(i, bgsamples) + 0.5) * self.get_bin_width(ifadc, bgsamples, ncollapse) + + t0shift * gkIntegrationWindow * self.get_bin_width(ifadc, bgsamples, ncollapse)) def get_signal(self, i, bgsamples, phecounts, ncollapse, bg, hwcorr): ifadc = i + bgsamples @@ -188,7 +384,7 @@ def main(): c0 = report.apply_zenith_correction(C_0, coszd) print("Zenith corrected c0:", c0) - c0 = report.apply_azimuth_correction(C_0, zenith, azimuth, zd_az_corr) + c0 = report.apply_azimuth_correction(C_0, zenith, azimuth) #, zd_az_corr) print("Zenith and Azimuth corrected c0:", c0) From d699a68c480b2ca3ac4a3b13ee3d060db2fd176a Mon Sep 17 00:00:00 2001 From: nzywucka Date: Thu, 28 Mar 2024 11:20:32 +0100 Subject: [PATCH 31/38] External files have been added --- ctapipe_io_magic/params.txt | 21 +++++++++ ctapipe_io_magic/switch_times.txt | 78 +++++++++++++++++++++++++++++++ ctapipe_io_magic/zd_az_corr.txt | 32 +++++++++++++ 3 files changed, 131 insertions(+) create mode 100644 ctapipe_io_magic/params.txt create mode 100644 ctapipe_io_magic/switch_times.txt create mode 100644 ctapipe_io_magic/zd_az_corr.txt diff --git a/ctapipe_io_magic/params.txt b/ctapipe_io_magic/params.txt new file mode 100644 index 0000000..a80df54 --- /dev/null +++ b/ctapipe_io_magic/params.txt @@ -0,0 +1,21 @@ +stime21 = -1 +stime22 = -1 +stime26 = -1 +stime261 = -1 +stime232 = -1 +stime233 = -1 +stime234 = -1 +stime27 = -1 +stime285 = -1 + +gkC_0Period1 = 16.43 # log (6.825E12*0.04*48*1.545E-6) - log(ptMAGIC) + +gkFullOverlap5 = 780.0 # start region of full overlap and end of ADC saturation (m) (after installation of new DAQ) +gkFullOverlap6 = 2000.0 # start region of full overlap and end of ADC saturation (m) (after installation of new DAQ) +gkFullOverlap7 = 1500.0 # start region of full overlap and end of ADC saturation (m) (after installation of new DAQ) + +gkIntegrationWindow = 47.954 # integration window of one FADC slice (in m): 0.299792458 (m/ns) / 1.00027 / 2 * 5 (ns) * 64 (output bins) +gkHWSwitchV1to4 = 96 # position of switch betw. amplitude sensing and p.e. counting (FADC slice) + +gkBGSamplesV1to4 = 16 # Number of BG samples before LIDAR signal (starting from 0) + diff --git a/ctapipe_io_magic/switch_times.txt b/ctapipe_io_magic/switch_times.txt new file mode 100644 index 0000000..91ec51b --- /dev/null +++ b/ctapipe_io_magic/switch_times.txt @@ -0,0 +1,78 @@ +Switch Time: 2018-01-12 12:00:00 +Correction Parameters: (0.311, 0.0) + +Switch Time: 2018-02-11 12:00:00 +Correction Parameters: (0.302, 0.016) + +Switch Time: 2018-06-27 12:00:00 +Correction Parameters: (0.323, 0.389) + +Switch Time: 2018-08-26 12:00:00 +Correction Parameters: (0.375, 0.0) + +Switch Time: 2018-09-13 12:00:00 +Correction Parameters: (0.394, 0.127) + +Switch Time: 2019-03-14 12:00:00 +Correction Parameters: (0.462, 0.11) + +Switch Time: 2019-05-08 12:00:00 +Correction Parameters: (0.489, 0.145) + +Switch Time: 2020-02-18 02:00:00 +Correction Parameters: (0.661, 0.0) + +Switch Time: 2020-02-20 12:00:00 +Correction Parameters: (0.573, 0.03) + +Switch Time: 2020-03-13 12:00:00 +Correction Parameters: (0.58, 0.75) + +Switch Time: 2020-06-10 00:00:00 +Correction Parameters: (0.637, 0.555) + +Switch Time: 2020-07-21 12:00:00 +Correction Parameters: (0.358, 0.5) + +Switch Time: 2020-07-26 01:22:00 +Correction Parameters: (0.53, 0.0) + +Switch Time: 2020-07-26 02:46:00 +Correction Parameters: (0.385, 0.0) + +Switch Time: 2020-07-26 03:10:00 +Correction Parameters: (0.447, 0.0) + +Switch Time: 2020-07-29 12:00:00 +Correction Parameters: (0.36, 0.1) + +Switch Time: 2020-09-09 12:00:00 +Correction Parameters: (0.99, 0.0) + +Switch Time: 2020-09-19 02:30:00 +Correction Parameters: (0.38, 0.001) + +Switch Time: 2020-09-19 23:10:00 +Correction Parameters: (0.451, 0.0) + +Switch Time: 2020-10-06 12:00:00 +Correction Parameters: (0.423, 0.209) + +Switch Time: 2021-02-10 12:00:00 +Correction Parameters: (0.502, 0.112) + +Switch Time: 2021-09-19 14:00:00 +Correction Parameters: (0.43, 0.03) + +Switch Time: 2022-01-01 14:00:00 +Correction Parameters: (0.0, 0.03) + +Switch Time: 2022-05-05 12:00:00 +Correction Parameters: (0.0, 0.03) + +Switch Time: 2023-05-17 12:00:00 +Correction Parameters: (0.0, 0.03) + +Switch Time: 2023-07-03 12:00:00 +Correction Parameters: (0.0, 0.03) + diff --git a/ctapipe_io_magic/zd_az_corr.txt b/ctapipe_io_magic/zd_az_corr.txt new file mode 100644 index 0000000..a5f6c5c --- /dev/null +++ b/ctapipe_io_magic/zd_az_corr.txt @@ -0,0 +1,32 @@ +Rule 1: +Zenith Threshold: 50 +Azimuth Threshold: [157.5, 202.5] +Correction Factor: 0.17 + +Rule 2: +Zenith Threshold: 40 +Azimuth Threshold: [135, 225] +Correction Factor: 0.24 + +Rule 3: +Zenith Threshold: 30 +Azimuth Threshold: [112.5, 247.5] +Correction Factor: 0.24 + +Rule 4: +Zenith Threshold: 20 +Azimuth Threshold: [90, 270] +Correction Factor: 0.24 + +Rule 5: +Zenith Threshold: 10 +Azimuth Threshold: [337.5, 360] +Correction Factor: 0.2 + +Rule 6: +Zenith Threshold Min: 10 +Zenith Threshold Max: 50 +Azimuth Threshold Min: 337.5 +Azimuth Threshold Max: 360 +Correction Factor: 0.11 + From fadf6d4292469179e8a6cbaf8dc5eaf6b5ee3959 Mon Sep 17 00:00:00 2001 From: nzywucka Date: Thu, 28 Mar 2024 11:50:51 +0100 Subject: [PATCH 32/38] Reading BadReport and State parameters added. --- ctapipe_io_magic/__init__.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/ctapipe_io_magic/__init__.py b/ctapipe_io_magic/__init__.py index 4d55c02..a8cb5e3 100644 --- a/ctapipe_io_magic/__init__.py +++ b/ctapipe_io_magic/__init__.py @@ -99,8 +99,8 @@ class ReportLaserContainer(Container): UniqueID = Field(List[Any], 'No.') Bits = Field(List[Any], 'ID') MJD = Field(List[Any], 'Modified Julian Date') -# BadReport = Field(List[Any], 'Bad Report') # to be improved -# State = Field(List[Any], 'State') # to be improved + BadReport = Field(List[Any], 'Bad Report') # to be improved + State = Field(List[Any], 'State') # to be improved IsOffsetCorrection = Field(List[Any], 'Is Offset Correction') IsOffsetFitted = Field(List[Any], 'Is Offset Fitted') IsBGCorrection = Field(List[Any], 'Is BG Correction') @@ -1108,8 +1108,8 @@ def parse_laser_info(self): 'MReportLaser.MReport.fBits', 'MTimeLaser.fMjd', 'MTimeLaser.fTime.fMilliSec', - #'MReportLaser.MReport.fBadReport', # to be improved - #'MReportLaser.MReport.fState', # to be improved + 'MReportLaser.MReport.fBadReport', # to be improved + 'MReportLaser.MReport.fState', # to be improved 'MReportLaser.fIsOffsetCorrection', 'MReportLaser.fIsOffsetFitted', 'MReportLaser.fIsBGCorrection', @@ -1181,8 +1181,8 @@ def parse_laser_info(self): laser = ReportLaserContainer() laser.UniqueID = laser_info_runh['MReportLaser.MReport.fUniqueID'] laser.Bits = laser_info_runh['MReportLaser.MReport.fBits'] - #laser.BadReport = laser_info_runh['MReportLaser.MReport.fBadReport'] # to be improved - #laser.State = laser_info_runh['MReportLaser.MReport.fState'] # to be improved + laser.BadReport = laser_info_runh['MReportLaser.MReport.fBadReport'] # to be improved + laser.State = laser_info_runh['MReportLaser.MReport.fState'] # to be improved laser.IsOffsetCorrection = laser_info_runh['MReportLaser.fIsOffsetCorrection'] laser.IsOffsetFitted = laser_info_runh['MReportLaser.fIsOffsetFitted'] laser.IsBGCorrection = laser_info_runh['MReportLaser.fIsBGCorrection'] From fce73742717fbfd25bbab419fdba8b33128d4087 Mon Sep 17 00:00:00 2001 From: nzywucka Date: Tue, 2 Apr 2024 18:56:44 +0200 Subject: [PATCH 33/38] Reading unique reports changed. --- ctapipe_io_magic/__init__.py | 264 +++++++++++++-------------- ctapipe_io_magic/tests/test_lidar.py | 51 ++---- 2 files changed, 148 insertions(+), 167 deletions(-) diff --git a/ctapipe_io_magic/__init__.py b/ctapipe_io_magic/__init__.py index a8cb5e3..6140d3c 100644 --- a/ctapipe_io_magic/__init__.py +++ b/ctapipe_io_magic/__init__.py @@ -96,66 +96,66 @@ class ReportLaserContainer(Container): """ Container for Magic laser parameters """ - UniqueID = Field(List[Any], 'No.') - Bits = Field(List[Any], 'ID') - MJD = Field(List[Any], 'Modified Julian Date') - BadReport = Field(List[Any], 'Bad Report') # to be improved - State = Field(List[Any], 'State') # to be improved - IsOffsetCorrection = Field(List[Any], 'Is Offset Correction') - IsOffsetFitted = Field(List[Any], 'Is Offset Fitted') - IsBGCorrection = Field(List[Any], 'Is BG Correction') - IsT0ShiftFitted = Field(List[Any], 'Is T0 Shift Fitted') - IsUseGDAS = Field(List[Any], 'Is Use GDAS') - IsUpwardMoving = Field(List[Any], 'Is Upward Moving') - OverShoot = Field(List[Any], 'Over Shoot') - UnderShoot = Field(List[Any], 'Under Shoot') - BGSamples = Field(List[np.float32], 'BG Samples') - Transmission3km = Field(List[np.float32], 'Transmission at 3 km') - Transmission6km = Field(List[np.float32], 'Transmission at 6 km') - Transmission9km = Field(List[np.float32], 'Transmission at 9 km') - Transmission12km = Field(List[np.float32], 'Transmission at 12 km') - Zenith = Field(List[Angle], 'Zenith angle', unit=u.deg) - Azimuth = Field(List[Angle], 'Azimuth angle', unit=u.deg) - FullOverlap = Field(List[np.float32], 'Full Overlap') - EndGroundLayer = Field(List[np.float32], 'End Ground Layer') - GroundLayerTrans = Field(List[np.float32], 'Ground Layer Trans') - Calimaness = Field(List[np.float32], 'Calimaness') - CloudLayerAlt = Field(List[np.float32], 'Altitude of cloud layer') - CloudLayerDens = Field(List[np.float32], 'Density of cloud layer') - Klett_k = Field(List[np.float32], 'Klett k') - PheCounts = Field(List[List[np.float32]], 'Phe Counts') - Offset = Field(List[np.float32], 'Offset') - Offset_Calculated = Field(List[np.float32], 'Offset calculated') - Offset_Fitted = Field(List[np.float32], 'Offset fitted') - Offset2 = Field(List[np.float32], 'Offset 2') - Offset3 = Field(List[np.float32], 'Offset 3') - Background1 = Field(List[np.float32], 'Background 1') - Background2 = Field(List[np.float32], 'Background 2') - BackgroundErr1 = Field(List[np.float32], 'Background error 1') - BackgroundErr2 = Field(List[np.float32], 'Background error 2') - RangeMax = Field(List[np.float32], 'Range max') - RangeMax_Clouds = Field(List[np.float32], 'Range max clouds') - ErrorCode = Field(List[Any], 'Error code') - ScaleHeight_fit = Field(List[np.float32], 'Scale Height fit') - Alpha_fit = Field(List[np.float32], 'Alpha fit') - Chi2Alpha_fit = Field(List[np.float32], 'Chi2 Alpha fit') - Alpha_firstY = Field(List[np.float32], 'Alpha first Y') - Alpha_Junge = Field(List[np.float32], 'Alpha Junge') - PBLHeight = Field(List[np.float32], 'PBL Height') - Chi2Full_fit = Field(List[np.float32], 'Chi2 Full fit') - SignalSamples = Field(List[np.float32], 'Signal Samples') - HWSwitch = Field(List[np.float32], 'HW Switch') - HWSwitchMaxOffset = Field(List[np.float32], 'HW Switch Max Offset') - NCollapse = Field(List[np.float32], 'N Collapse') - Shots = Field(List[np.float32], 'Shots') - T0Shift = Field(List[np.float32], 'T0 Shift') - Interval_0 = Field(List[np.float32], 'Interval 0') - RCS_min_perfect = Field(List[np.float32], 'RCS min perfect') - RCS_min_clouds = Field(List[np.float32], 'RCS min cloud') - RCS_min_mol = Field(List[np.float32], 'RCS min mol') - LIDAR_ratio = Field(List[np.float32], 'LIDAR ratio') - LIDAR_ratio_Cloud = Field(List[np.float32], 'LIDAR ratio cloud') - LIDAR_ratio_Junge = Field(List[np.float32], 'LIDAR ratio Junge') + 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""" @@ -1108,8 +1108,8 @@ def parse_laser_info(self): 'MReportLaser.MReport.fBits', 'MTimeLaser.fMjd', 'MTimeLaser.fTime.fMilliSec', - 'MReportLaser.MReport.fBadReport', # to be improved - 'MReportLaser.MReport.fState', # to be improved + 'MReportLaser.MReport.fBadReport', + 'MReportLaser.MReport.fState', 'MReportLaser.fIsOffsetCorrection', 'MReportLaser.fIsOffsetFitted', 'MReportLaser.fIsBGCorrection', @@ -1165,82 +1165,82 @@ def parse_laser_info(self): 'MReportLaser.fLIDAR_ratio_Junge', ] - unique_reports = {} + old_mjd = None + unique_reports = [] for rootf in self.files_: try: laser_info_runh = rootf['Laser'].arrays( - laser_info_array_list_runh, library="np" - ) + laser_info_array_list_runh, library="np") mjd_value = laser_info_runh['MTimeLaser.fMjd'] millisec_value = laser_info_runh['MTimeLaser.fTime.fMilliSec'] + 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] - for mjd_values, millisec_values in zip(mjd_value, millisec_value): - unique_key = (mjd_values, millisec_values) - if unique_key not in unique_reports: - unique_reports[unique_key] = [] - laser = ReportLaserContainer() - laser.UniqueID = laser_info_runh['MReportLaser.MReport.fUniqueID'] - laser.Bits = laser_info_runh['MReportLaser.MReport.fBits'] - laser.BadReport = laser_info_runh['MReportLaser.MReport.fBadReport'] # to be improved - laser.State = laser_info_runh['MReportLaser.MReport.fState'] # to be improved - laser.IsOffsetCorrection = laser_info_runh['MReportLaser.fIsOffsetCorrection'] - laser.IsOffsetFitted = laser_info_runh['MReportLaser.fIsOffsetFitted'] - laser.IsBGCorrection = laser_info_runh['MReportLaser.fIsBGCorrection'] - laser.IsT0ShiftFitted = laser_info_runh['MReportLaser.fIsT0ShiftFitted'] - laser.IsUseGDAS = laser_info_runh['MReportLaser.fIsUseGDAS'] - laser.IsUpwardMoving = laser_info_runh['MReportLaser.fIsUpwardMoving'] - laser.OverShoot = int(laser_info_runh['MReportLaser.fOverShoot']) - laser.UnderShoot = int(laser_info_runh['MReportLaser.fUnderShoot']) - laser.BGSamples = int(laser_info_runh['MReportLaser.fBGSamples']) - laser.Transmission3km = laser_info_runh['MReportLaser.fTransmission3km'] - laser.Transmission6km = laser_info_runh['MReportLaser.fTransmission6km'] - laser.Transmission9km = laser_info_runh['MReportLaser.fTransmission9km'] - laser.Transmission12km = laser_info_runh['MReportLaser.fTransmission12km'] - laser.Zenith = laser_info_runh['MReportLaser.fZenith']* u.deg - laser.Azimuth = laser_info_runh['MReportLaser.fAzimuth']* u.deg - laser.FullOverlap = laser_info_runh['MReportLaser.fFullOverlap'] - laser.EndGroundLayer = laser_info_runh['MReportLaser.fEndGroundLayer'] - laser.GroundLayerTrans = laser_info_runh['MReportLaser.fGroundLayerTrans'] - laser.Klett_k = laser_info_runh['MReportLaser.fKlett_k'] - laser.PheCounts = laser_info_runh['MReportLaser.fPheCounts'] - laser.Calimaness = laser_info_runh['MReportLaser.fCalimaness'] - laser.CloudLayerAlt = laser_info_runh['MReportLaser.fCloudLayerAlt'] - laser.CloudLayerDens = laser_info_runh['MReportLaser.fCloudLayerDens'] - laser.Offset = laser_info_runh['MReportLaser.fOffset'] - laser.Offset_Calculated = laser_info_runh['MReportLaser.fOffset_Calculated'] - laser.Offset_Fitted = laser_info_runh['MReportLaser.fOffset_Fitted'] - laser.Offset2 = laser_info_runh['MReportLaser.fOffset2'] - laser.Offset3 = laser_info_runh['MReportLaser.fOffset3'] - laser.Background1 = laser_info_runh['MReportLaser.fBackground1'] - laser.Background2 = laser_info_runh['MReportLaser.fBackground2'] - laser.BackgroundErr1 = laser_info_runh['MReportLaser.fBackgroundErr1'] - laser.BackgroundErr2 = laser_info_runh['MReportLaser.fBackgroundErr2'] - laser.RangeMax = laser_info_runh['MReportLaser.fRangeMax'] - laser.RangeMax_Clouds = laser_info_runh['MReportLaser.fRangeMax_Clouds'] - laser.ErrorCode = laser_info_runh['MReportLaser.fErrorCode'] - laser.ScaleHeight_fit = laser_info_runh['MReportLaser.fScaleHeight_fit'] - laser.Alpha_fit = laser_info_runh['MReportLaser.fChi2Alpha_fit'] - laser.Chi2Alpha_fit = laser_info_runh['MReportLaser.fChi2Alpha_fit'] - laser.Alpha_firstY = laser_info_runh['MReportLaser.fAlpha_firstY'] - laser.Alpha_Junge = laser_info_runh['MReportLaser.fAlpha_Junge'] - laser.PBLHeight = laser_info_runh['MReportLaser.fPBLHeight'] - laser.Chi2Full_fit = laser_info_runh['MReportLaser.fChi2Full_fit'] - laser.HWSwitchMaxOffset = laser_info_runh['MReportLaser.fHWSwitchMaxOffset'] - laser.NCollapse = laser_info_runh['MReportLaser.fNCollapse'] - laser.Shots = laser_info_runh['MReportLaser.fShots'] - laser.T0Shift = laser_info_runh['MReportLaser.fT0Shift'] - laser.Interval_0 = laser_info_runh['MReportLaser.fInterval_0'] - laser.RCS_min_perfect = laser_info_runh['MReportLaser.fRCS_min_perfect'] - laser.RCS_min_clouds = laser_info_runh['MReportLaser.fRCS_min_clouds'] - laser.RCS_min_mol = laser_info_runh['MReportLaser.fRCS_min_mol'] - laser.LIDAR_ratio = laser_info_runh['MReportLaser.fLIDAR_ratio'] - laser.LIDAR_ratio_Cloud = laser_info_runh['MReportLaser.fLIDAR_ratio_Cloud'] - laser.LIDAR_ratio_Junge = laser_info_runh['MReportLaser.fLIDAR_ratio_Junge'] - - millisec_seconds = millisec_values * 1e-3 - combined_mjd_value = mjd_values + millisec_seconds / 86400 - laser.MJD = combined_mjd_value - unique_reports[unique_key].append(laser) except KeyError as e: print(f"Required key not found in the file {rootf}: {e}") continue diff --git a/ctapipe_io_magic/tests/test_lidar.py b/ctapipe_io_magic/tests/test_lidar.py index fe7268b..f36a685 100644 --- a/ctapipe_io_magic/tests/test_lidar.py +++ b/ctapipe_io_magic/tests/test_lidar.py @@ -6,43 +6,24 @@ 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", -] - -data_dict = dict() - -data_dict["20210314_M1_05095172.002_Y_CrabNebula-W0.40+035.root"] = dict() - -data_dict["20210314_M1_05095172.002_Y_CrabNebula-W0.40+035.root"]["n_events_tot"] = 500 -data_dict["20210314_M1_05095172.002_Y_CrabNebula-W0.40+035.root"]["n_events_stereo"] = 452 -data_dict["20210314_M1_05095172.002_Y_CrabNebula-W0.40+035.root"]["n_events_pedestal"] = 48 -data_dict["20210314_M1_05095172.002_Y_CrabNebula-W0.40+035.root"]["n_events_mc_mono"] = 0 - -def test_sample_report_laser(sample_report_laser): - assert sample_report_laser is not None +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, ReportLaserContainer - dataset = (test_calibrated_real_dir/"20210314_M1_05095172.002_Y_CrabNebula-W0.40+035.root") - with MAGICEventSource(input_url=dataset) as source: - print(ReportLaserContainer()) - return ReportLaserContainer() - -@pytest.mark.parametrize("dataset", test_calibrated_real) -def test_lidar_parameters(dataset): from ctapipe_io_magic import MAGICEventSource - with MAGICEventSource(input_url=dataset) as source: + with MAGICEventSource(input_url=test_calibrated_real,process_run=True) as source: laser = source.laser - for key, report_list in laser.items(): - for report in report_list: - assert report.Transmission12km == pytest.approx(0.89, abs=8.9e-07) - assert report.IsUseGDAS == pytest.approx([False]) # bool - azimuth_degrees = report.Azimuth.value # Convert Angle to degrees - assert azimuth_degrees == pytest.approx(276.63) # Angle - flattened_array = np.concatenate(report.PheCounts) - phe_counts = [int(str(value).replace('\n', '').strip()) for value in flattened_array] - assert 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.Bits == (50331656) # Any - + 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: + 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) == [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.Bits == 50331656 + assert report.MJD == pytest.approx(59286.91349633, rel = 1e-10, abs = 1e-6) From 25136f0c81d287a150fc2c24a84c00bcd6a55844 Mon Sep 17 00:00:00 2001 From: nzywucka Date: Tue, 2 Apr 2024 19:11:34 +0200 Subject: [PATCH 34/38] Remove assert from tests in codacy `# Please enter the commit message for your changes. Lines starting --- .codacy.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.codacy.yml b/.codacy.yml index 260dd03..2c99101 100644 --- a/.codacy.yml +++ b/.codacy.yml @@ -6,3 +6,5 @@ engines: pyflakes: disable: - F999 +assert_used: skips: ['*_test.py', '*test_*.py'] + From 9f14bb55fd69956bce72fac65e36035294e611d6 Mon Sep 17 00:00:00 2001 From: nzywucka Date: Tue, 2 Apr 2024 19:17:55 +0200 Subject: [PATCH 35/38] Skipping assert checks in the test files. Trailing whitespaces and unused libraries have been removed. --- .codacy.yml | 3 ++- ctapipe_io_magic/__init__.py | 5 ++--- ctapipe_io_magic/tests/test_lidar.py | 1 - 3 files changed, 4 insertions(+), 5 deletions(-) diff --git a/.codacy.yml b/.codacy.yml index 2c99101..ed42fe1 100644 --- a/.codacy.yml +++ b/.codacy.yml @@ -3,8 +3,9 @@ engines: pylint: enabled: true python_version: 3 + assert_used: skips: ['*_test.py', '*test_*.py'] pyflakes: disable: - F999 -assert_used: skips: ['*_test.py', '*test_*.py'] + assert_used: skips: ['*_test.py', '*test_*.py'] diff --git a/ctapipe_io_magic/__init__.py b/ctapipe_io_magic/__init__.py index 6140d3c..362ae95 100644 --- a/ctapipe_io_magic/__init__.py +++ b/ctapipe_io_magic/__init__.py @@ -15,7 +15,6 @@ from decimal import Decimal from astropy import units as u from astropy.time import Time -from astropy.coordinates import Angle from pkg_resources import resource_filename from ctapipe.io import EventSource, DataLevel @@ -1175,8 +1174,8 @@ def parse_laser_info(self): millisec_value = laser_info_runh['MTimeLaser.fTime.fMilliSec'] 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: + 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] diff --git a/ctapipe_io_magic/tests/test_lidar.py b/ctapipe_io_magic/tests/test_lidar.py index f36a685..bb920fd 100644 --- a/ctapipe_io_magic/tests/test_lidar.py +++ b/ctapipe_io_magic/tests/test_lidar.py @@ -1,5 +1,4 @@ import os -import numpy as np from pathlib import Path import pytest From 20091c743283b4108894738057da6f093f36c380 Mon Sep 17 00:00:00 2001 From: nzywucka Date: Tue, 2 Apr 2024 19:23:20 +0200 Subject: [PATCH 36/38] Skipping assert checks in the test files. --- .codacy.yml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/.codacy.yml b/.codacy.yml index ed42fe1..55a070d 100644 --- a/.codacy.yml +++ b/.codacy.yml @@ -3,9 +3,8 @@ engines: pylint: enabled: true python_version: 3 - assert_used: skips: ['*_test.py', '*test_*.py'] pyflakes: disable: - F999 - assert_used: skips: ['*_test.py', '*test_*.py'] + assert_used: skips: ['*/*/test_*.py', '*/test_*.py', 'test_*.py'] From 4f8597c91476dcbe06f4a2278e86196d125c630b Mon Sep 17 00:00:00 2001 From: nzywucka Date: Wed, 3 Apr 2024 09:56:42 +0200 Subject: [PATCH 37/38] Tests for reading lidar parameters have been removed. --- .../tests/test_magic_event_source.py | 44 +------------------ 1 file changed, 1 insertion(+), 43 deletions(-) diff --git a/ctapipe_io_magic/tests/test_magic_event_source.py b/ctapipe_io_magic/tests/test_magic_event_source.py index f57988f..fbac975 100644 --- a/ctapipe_io_magic/tests/test_magic_event_source.py +++ b/ctapipe_io_magic/tests/test_magic_event_source.py @@ -286,7 +286,7 @@ def test_allowed_tels(): dataset = ( test_calibrated_real_dir - / "20210314_M1_05095172.002_Y_CrabNebula-W0.40+035.root" #"20210314_M1_05095172.001_Y_CrabNebula-W0.40+035.root" + / "20210314_M1_05095172.001_Y_CrabNebula-W0.40+035.root" ) allowed_tels = {1} with MAGICEventSource( @@ -559,45 +559,3 @@ def test_broken_subruns_missing_arrays(): input_url=input_file, process_run=True, ) -# def test_eventseeker(): -# dataset = get_dataset_path("20131004_M1_05029747.003_Y_MagicCrab-W0.40+035.root") -# -# with MAGICEventSource(input_url=dataset) as source: -# seeker = EventSeeker(source) -# event = seeker.get_event_index(0) -# assert event.count == 0 -# assert event.index.event_id == 29795 -# -# event = seeker.get_event_index(2) -# assert event.count == 2 -# assert event.index.event_id == 29798 -# -# def test_eventcontent(): -# dataset = get_dataset_path("20131004_M1_05029747.003_Y_MagicCrab-W0.40+035.root") -# -# with MAGICEventSource(input_url=dataset) as source: -# seeker = EventSeeker(source) -# event = seeker.get_event_index(0) -# assert event.dl1.tel[1].image[0] == -0.53125 -# assert event.dl1.tel[1].peak_time[0] == 49.125 - - -#@pytest.mark.parametrize("dataset", test_calibrated_all) -#def test_lidar_parameters(dataset): -# try: -# with MAGICEventSource(input_url=dataset) as source: -# params = ReportLaserContainer() -# for item, event in params: -# #params.parse_laser_info() -# # assert params.laser.MJD == 59286.91349633102 -# assert params.Transmission12km == 0.89 -# # assert params.laser.Azimuth == u.Quantity(277, u.deg) -# #assert params.laser.BadReport == False -# except KeyInFileError: -# pytest.skip("Skipping test for file without 'Laser' key.") - -#dataset = ( -# test_calibrated_real_dir -# / "20210314_M1_05095172.002_Y_CrabNebula-W0.40+035.root" -#) - From 5221ce41cdf056a5f4bf5dadf2220ffe1c60195b Mon Sep 17 00:00:00 2001 From: nzywucka Date: Wed, 30 Oct 2024 10:37:02 +0100 Subject: [PATCH 38/38] laser_test_new.py: improved implementation of corrections of LIDAR parameters azimuth_conditions.txt: parameters to perform azimuth corrections switch_times.txt and switchtimes1.txt: parameters to perform time corrections --- ctapipe_io_magic/azimuth_conditions.txt | 15 ++ ctapipe_io_magic/laser_test_new.py | 306 ++++++++++++++++++++++++ ctapipe_io_magic/switch_times.txt | 113 +++------ ctapipe_io_magic/switchtimes1.txt | 16 ++ 4 files changed, 372 insertions(+), 78 deletions(-) create mode 100644 ctapipe_io_magic/azimuth_conditions.txt create mode 100644 ctapipe_io_magic/laser_test_new.py create mode 100644 ctapipe_io_magic/switchtimes1.txt diff --git a/ctapipe_io_magic/azimuth_conditions.txt b/ctapipe_io_magic/azimuth_conditions.txt new file mode 100644 index 0000000..0139515 --- /dev/null +++ b/ctapipe_io_magic/azimuth_conditions.txt @@ -0,0 +1,15 @@ +# zenith conditions (1,2) +# azimuth conditions (3,4) +# correction (5) +0. 50. 157.5 202.5 0.17 +0. 40. 135. 225. 0.24 +0. 30. 112.5 247.5 0.24 +0. 20. 90. 270. 0.24 +10. 90. 337.5 360. 0.2 +10. 50. 292.5 360. 0.11 +10. 40. 0. 4.5 0.11 +10. 30. 0. 67.5 0.11 +20. 40. 247.5 270. 0.11 +15. 50. 225. 247.5 0.05 +30. 60. 202.5 225. 0.05 + diff --git a/ctapipe_io_magic/laser_test_new.py b/ctapipe_io_magic/laser_test_new.py new file mode 100644 index 0000000..e45deef --- /dev/null +++ b/ctapipe_io_magic/laser_test_new.py @@ -0,0 +1,306 @@ +import datetime +import numpy as np +import astropy.units as u +from scipy.stats import t +from scipy.stats import iqr +from scipy.optimize import curve_fit +from ctapipe_io_magic import MAGICEventSource +import matplotlib.pyplot as plt +import matplotlib.dates as mdates + +import logging +from pathlib import Path +import os + +logger = logging.getLogger(__name__) +logger.addHandler(logging.StreamHandler()) +logger.setLevel(logging.INFO) + +class MReportLaser: + def __init__(self): + self.gkRangeMax2 = 31000. + self.gkShotsNewLaser = 25000 + self.gkShotsOldLaser = 50000 + self.fRangeMax = 0.0 + self.fShots = 0 + self.fRangeMax_Clouds = 0.0 + self.fZenith = 50.0 * u.deg + self.fAzimuth = 200.0 * u.deg + self.fLog = None + self.Alphafit_corr = 0. + + self.gkIntegrationWindow = 47.954 + self.gkHWSwitchV1to4 = 96 + self.gkSecsPerYear = 3600.*24*365 + self.gkC_0ZenithCorr = 0.038 + self.gkBins = 512 + self.gkMaxCloudLayers = 3 + self.gkUpdate1Vers = 201206050 + self.gkUpdate2Vers = 201302221 + self.gkUpdate3Vers = 201304250 + self.gkUpdate4Vers = 201401091 + self.gkUpdate5Vers = 201705110 + self.gkC_0Elterman = 0.081 + self.coszd = np.cos(self.fZenith) + + """ + Code completion, need to decide which part of it should be kept + """ + + self.fCloudHM = [0.] * self.gkMaxCloudLayers + self.fCloudHStd = [0.] * self.gkMaxCloudLayers + self.fCloudLR = [0.] * self.gkMaxCloudLayers + self.fCloudFWHM = [0.] * self.gkMaxCloudLayers + self.fCloudBase = [0.] * self.gkMaxCloudLayers + self.fCloudTop = [0.] * self.gkMaxCloudLayers + self.fCloudTopT = [0.] * self.gkMaxCloudLayers + self.fCloudTrans = [1.] * self.gkMaxCloudLayers + self.fPheCounts = [0.] * self.gkBins + self.fCalimaness = 0. + self.fCloudLayerAlt = 0. + self.fCloudLayerDens = 0. + self.fTransmission3km = 0. + self.fTransmission6km = 0. + self.fTransmission9km = 0. + self.fTransmission12km = 0. + self.fErrorCode = 0. + self.gkBGSamplesV1to4 = 16 + self.gkBGSamplesV5 = 31 + self.gkSignalSamplesV5 = 480 + self.gkHWSwitchV5 = 256 + self.fBGSamples = self.gkBGSamplesV5 + self.fSignalSamples = self.gkSignalSamplesV5 + self.fHWSwitch = self.gkHWSwitchV5 + self.gkC_0TempCorr3 = 0.0002 + self.temp_mean = 8.89 + self.gkC_0HumCorr = 0.00023 + + self.ifadc = 1 + self.collapse = 1 + + self.signalsamples = self.fSignalSamples + self.fIsBGCorrection = 0 + self.fBackground1 = 0.0 + self.fBackground2 = 0.0 + self.fBackgroundErr1 = 0.0 + self.fBackgroundErr2 = 0.0 + self.fNCollapse = 4 + + def magic_reports_reading(self, input_file, process_run=False): + """Read laser parameters from MAGIC files""" + logger.info(f"\nInput file: {input_file}") + event_source = MAGICEventSource(input_file, process_run=process_run) + laser = event_source.laser + weather = event_source.weather + return laser, weather #{'laser': laser, 'weather': weather} + + def reset_clouds(self): + for arr in (self.fCloudHM, self.fCloudHStd, self.fCloudLR, self.fCloudFWHM, + self.fCloudBase, self.fCloudTop, self.fCloudTopT, self.fCloudTrans): + arr[:] = [0.] * self.gkMaxCloudLayers + + def interprete_body(self, str, ver): + interpreters = { + self.gkUpdate1Vers: self.interprete_clouds, + self.gkUpdate4Vers: self.interprete_cloud_base, + self.gkUpdate3Vers: self.interprete_error_code, + self.gkUpdate1Vers: self.interprete_pointing, + } + for i in range(self.gkBins): + n = sscanf(str, "%08x" if ver >= self.gkUpdate1Vers else "%06x") + if n != 1: + self.fLog << _warn_ << GetDescriptor() << f": could not read photon counts for bin {i}, instead got: {str}" << endl + return kCONTINUE + str = str[6 if ver < self.gkUpdate1Vers else 8:] + + for v, interpreter in interpreters.items(): + if ver >= v: + if not interpreter(str): + return kCONTINUE + + if str != "OVER": + self.fLog << _warn_ << f"WARNING - 'OVER' tag not found (instead: {str})" << endl + return kTRUE + + def interprete_clouds(self, str): + self.fCalimaness, self.fCloudLayerAlt, self.fCloudLayerDens, str = map(float, str.split(' ', 3)) + return kTRUE + + def interprete_cloud_base(self, str): + self.fCloudLayerAlt, str = map(float, str.split(' ', 1)) + return kTRUE + + def interprete_error_code(self, str): + self.fErrorCode, str = map(int, str.split(' ', 1)) + return kTRUE + + def interprete_pointing(self, str): + self.fZenith, self.fAzimuth, str = map(float, str.split(' ', 2)) + return kTRUE + + def print(self, o): + self.fLog << GetDescriptor() << ":" << endl + for i, phe in enumerate(self.fPheCounts): + self.fLog << f"i={i} binw: {GetBinWidth(i):.0f} range: {(GetOffsetR(i) + gkIntegrationWindow * (GetFADCR(i - self.fBGSamples) + 0.5) * GetBinWidth(i)) * 0.001:.2f} km counts: {phe / GetBinWidth(i)}" << endl + + self.fLog << endl + self.fLog << f"Calimaness: {self.fCalimaness}, Altiude of first cloud layer: {self.fCloudLayerAlt}, Optical thickness cloud layer: {self.fCloudLayerDens}" << endl + self.fLog << f"Online analysis: Transmission3km: {self.fTransmission3km}, Transmission6km: {self.fTransmission6km}, Transmission9km: {self.fTransmission9km}, Transmission12km: {self.fTransmission12km}" << endl + self.fLog << endl + + def eval_gdas(self, x, par): + hlo, costheta, logPTMAGIC, h = par[1:] + return par[0] - logPTMAGIC + self.fGDAS.EvalF(hlo, h, costheta, fBeta_mol_0) + + def read_switch_times(self, filename): + overlaps = {} + C_0s = {} + with open(filename, 'r') as file: + for line in file: + if line[0] != "#": + parts = line.strip().split() + if len(parts) >= 4: + dt_str = ' '.join(parts[:2]) + overlap = float(parts[2].rstrip('.')) + C_0 = float(parts[3]) + self.switch_time = datetime.datetime.strptime(dt_str, '%Y-%m-%d %H:%M:%S') + overlaps[self.switch_time] = overlap + C_0s[self.switch_time] = C_0 + return overlaps, C_0s + + def read_zd_az_conditions(self, filename): + with open(filename, 'r') as file: + azimuth_conditions = [list(map(float, line.strip().split())) for line in file if not line.startswith('#')] + return azimuth_conditions + + """ + C_0 parameter calculation based on time + """ + def GetC0(self, time, overlaps, C_0s): + if time: + for switchtime, overlap in overlaps.items(): + stime = int((time - switchtime).total_seconds()) + if stime > 0: + FullOverlap = overlap + C_0 = C_0s[switchtime] + + self.fRangeMax = 40000. if self.fZenith > 70. * u.deg else self.gkRangeMax2 + self.fShots = self.gkShotsNewLaser + self.fRangeMax_Clouds = 23000. + + self.fRangeMax_Clouds = np.min([(self.fRangeMax - 100.) * np.cos(self.fZenith), self.fRangeMax_Clouds]) + + return C_0, FullOverlap + + def ApplyZenithCorrection(self, c0): + print("C_0_prev:", c0) + zd_corr = np.log(1 - self.gkC_0ZenithCorr * (1 - self.coszd)) + c0 += 2 * zd_corr + return True + + def ApplyAzimuthCorrection(self, c0, zd_az_conditions): + for condition in zd_az_conditions: + zenith_min, zenith_max, azimuth_min, azimuth_max, correction = condition + if (zenith_min <= self.fZenith.value <= zenith_max and + azimuth_min <= self.fAzimuth.value <= azimuth_max): + c0 -= 2 * correction + return True + return False + + def ApplyTimeCorrection(self, c0, time, c1, c2, zd_az_conditions): + if time: + for switchtime, corr1 in c1.items(): + stime = int((time - switchtime).total_seconds()) + if stime < 0: + correction1 = corr1 + correction2 = c2[switchtime] + if (time - switchtime28).total_seconds() < 0 < (time - switchtime275).total_seconds(): + self.ApplyAzimuthCorrection(c0, zd_az_conditions) + c0 -= 2 * (correction1 + correction2 * prev_stime / self.gkSecsPerYear) + 2.0 / self.coszd * self.Alphafit_corr + return c0 + prev_stime = stime + return c0 + + def apply_temperature_correction(self, c0, temperature): + """valid since 17.11.2016""" + c0 -= self.gkC_0TempCorr3 * (temperature - self.temp_mean) + + return True, c0 + + def apply_humidity_correction(self, c0, humidity): + if humidity < 0: + return False + + c0 -= self.gkC_0HumCorr * (humidity - 30.0) + + return True, c0 + + def test_c0_drift(self, c1, c2, zd_az_conditions, filename): + time0 = datetime.datetime(2019, 1, 1, 0, 0, 0) + times = [] + c0_values = [] + with open(filename, 'w') as file: + for i in range(int(5.5 * 365)): + c0 = 0. + time1 = time0 + datetime.timedelta(days=i) + c0 = self.ApplyTimeCorrection(c0, time1, c1, c2, zd_az_conditions) + times.append(time1) + c0_values.append(c0) + file.write(f"{time1.year}-{time1.month}-{time1.day} {time1.hour}:{time1.minute} {c0}\n") + print(f"{time1.year}-{time1.month}-{time1.day} {time1.hour}:{time1.minute} {c0:.5f}") + plt.rcParams['font.size'] = 20 + fig, ax = plt.subplots() + ax.plot(times, c0_values, marker='o', linestyle='-', label='C0 Drift') + + ax.xaxis.set_major_locator(mdates.YearLocator()) + ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d')) + ax.set_xlabel('Time') + ax.set_ylabel('C_0 Correction') + #ax.legend() + plt.xticks() #(rotation=45) + plt.tight_layout() + plt.savefig("c0_python.png") + plt.show() + +# Sample usage: +def main(): + input_file = '/home/zywuckan/ctapipe_io_magic/ctapipe_io_magic/tests/test_data/real/calibrated/20210314_M1_05095172.002_Y_CrabNebula-W0.40+035.root' + report = MReportLaser() + laser_list, weather_list = report.magic_reports_reading(input_file) # Unpack the tuple + first_laser = laser_list[0] + first_weather = weather_list[0] + time = datetime.datetime(2020, 7, 26, 3, 10, 0) + zenith = first_laser.Zenith #['Zenith'] + azimuth = first_laser['Azimuth'] + temperature = first_weather['Temperature'] #* u.deg_C + print(temperature) + humidity = first_weather['Humidity'] + print(humidity) + i = 0 + overlaps, C_0s = report.read_switch_times('switchtimes1.txt') + c0, _ = report.GetC0(time, overlaps, C_0s) + c0, overlap = report.GetC0(time, overlaps, C_0s) + print("C0:", c0) + print("overlap:", overlap) + + zenith_corr = report.ApplyZenithCorrection(c0) + print(zenith_corr) + + zd_az_conditions = report.read_zd_az_conditions('azimuth_conditions.txt') + azimuth_corr = report.ApplyAzimuthCorrection(c0, zd_az_conditions) + print("azimuth correction:", azimuth_corr) + + c1, c2 = report.read_switch_times('switch_times.txt') + time_corr2 = report.ApplyTimeCorrection(c0, time, c1, c2, zd_az_conditions) + print("time correction:", time_corr2) + + temp = report.apply_temperature_correction(c0, temperature) + print(temp) + + humid = report.apply_humidity_correction(c0, humidity) + print(humid) + + report.test_c0_drift(c1, c2, zd_az_conditions, 'c0_drift_data.txt') + +if __name__ == "__main__": + main() diff --git a/ctapipe_io_magic/switch_times.txt b/ctapipe_io_magic/switch_times.txt index 91ec51b..4830fd6 100644 --- a/ctapipe_io_magic/switch_times.txt +++ b/ctapipe_io_magic/switch_times.txt @@ -1,78 +1,35 @@ -Switch Time: 2018-01-12 12:00:00 -Correction Parameters: (0.311, 0.0) - -Switch Time: 2018-02-11 12:00:00 -Correction Parameters: (0.302, 0.016) - -Switch Time: 2018-06-27 12:00:00 -Correction Parameters: (0.323, 0.389) - -Switch Time: 2018-08-26 12:00:00 -Correction Parameters: (0.375, 0.0) - -Switch Time: 2018-09-13 12:00:00 -Correction Parameters: (0.394, 0.127) - -Switch Time: 2019-03-14 12:00:00 -Correction Parameters: (0.462, 0.11) - -Switch Time: 2019-05-08 12:00:00 -Correction Parameters: (0.489, 0.145) - -Switch Time: 2020-02-18 02:00:00 -Correction Parameters: (0.661, 0.0) - -Switch Time: 2020-02-20 12:00:00 -Correction Parameters: (0.573, 0.03) - -Switch Time: 2020-03-13 12:00:00 -Correction Parameters: (0.58, 0.75) - -Switch Time: 2020-06-10 00:00:00 -Correction Parameters: (0.637, 0.555) - -Switch Time: 2020-07-21 12:00:00 -Correction Parameters: (0.358, 0.5) - -Switch Time: 2020-07-26 01:22:00 -Correction Parameters: (0.53, 0.0) - -Switch Time: 2020-07-26 02:46:00 -Correction Parameters: (0.385, 0.0) - -Switch Time: 2020-07-26 03:10:00 -Correction Parameters: (0.447, 0.0) - -Switch Time: 2020-07-29 12:00:00 -Correction Parameters: (0.36, 0.1) - -Switch Time: 2020-09-09 12:00:00 -Correction Parameters: (0.99, 0.0) - -Switch Time: 2020-09-19 02:30:00 -Correction Parameters: (0.38, 0.001) - -Switch Time: 2020-09-19 23:10:00 -Correction Parameters: (0.451, 0.0) - -Switch Time: 2020-10-06 12:00:00 -Correction Parameters: (0.423, 0.209) - -Switch Time: 2021-02-10 12:00:00 -Correction Parameters: (0.502, 0.112) - -Switch Time: 2021-09-19 14:00:00 -Correction Parameters: (0.43, 0.03) - -Switch Time: 2022-01-01 14:00:00 -Correction Parameters: (0.0, 0.03) - -Switch Time: 2022-05-05 12:00:00 -Correction Parameters: (0.0, 0.03) - -Switch Time: 2023-05-17 12:00:00 -Correction Parameters: (0.0, 0.03) - -Switch Time: 2023-07-03 12:00:00 -Correction Parameters: (0.0, 0.03) - +# date up to correction is valid +# time +# correction parameter 1 +# correction Parameter 2 +2017-12-10 12:00:00 0.146 0.389 #stime16 +2018-01-12 12:00:00 0.304 0.553 #stime17 +2018-02-11 12:00:00 0.311 0.0 #stime18 +2018-06-27 12:00:00 0.302 0.16 #stime19 +2018-08-26 12:00:00 0.323 0.389 #stime20 +2018-09-13 12:00:00 0.375 0.0 #stime201 +2019-03-14 12:00:00 0.394 0.127 #stime21 +2019-05-08 12:00:00 0.462 0.11 #stime22 +2020-02-18 02:00:00 0.489 0.145 #stime229 +2020-02-20 12:00:00 0.661 0.0 #stime23 +2020-03-13 12:00:00 0.573 0.03 #stime2300 +2020-06-10 00:00:00 0.58 0.75 #stime230 +2020-07-21 12:00:00 0.637 0.555 #stime231 +2020-07-26 01:22:00 0.358 0.5 #stime232 +2020-07-26 02:46:00 0.53 0.0 #stime233 +2020-07-26 03:10:00 0.385 0.0 #stime234 +2020-07-29 12:00:00 0.447 0.0 #stime24 +2020-09-09 12:00:00 0.36 0.1 #stime25 +2020-09-19 02:30:00 0.99 0.0 #stime26 +2020-09-19 23:10:00 0.38 0.001 #stime261 +2020-10-06 12:00:00 0.451 0.0 #stime262 +2021-02-10 12:00:00 0.423 0.209 #stime263 +2021-09-19 14:00:00 0.502 0.112 #stime27 +#2022-01-01 14:00:00 0.43 0.03 #stime275 +2022-05-05 12:00:00 0.439 0.03 #stime28 +#2022-05-05 12:00:00 0.439 0.03 #stime28 calculated as 0.43 + 0.03 * (switchtime275 - switchtime27) +#2023-05-17 12:00:00 0.449 0.03 #stime285 calculated as 0.439 + 0.03 * (switchtime28 - switchtime275) +#2023-05-17 12:00:00 0.449 0.03 #stime285 +#2023-07-03 12:00:00 0.0 0.03 #stime29 calculated as 0.449 + 0.03 * (switchtime285 - switchtime28) +2023-07-03 12:00:00 0.0 0.03 #stime29 +2025-01-01 12:00:00 0.20 0.03 #final calculations diff --git a/ctapipe_io_magic/switchtimes1.txt b/ctapipe_io_magic/switchtimes1.txt new file mode 100644 index 0000000..eb0bb7e --- /dev/null +++ b/ctapipe_io_magic/switchtimes1.txt @@ -0,0 +1,16 @@ +# date and time +# fullOverlap +# C_0 +2017-02-13 12:00:00 340. 18.87 +2017-04-02 12:00:00 420. 18.87 +2019-03-22 12:00:00 780. 18.87 +2019-11-02 12:00:00 400. 18.87 +2020-07-26 01:22:00 2000. 18.87 +2020-07-26 02:46:00 780. 18.87 +2020-07-26 03:10:00 1500. 18.87 +2020-09-19 02:00:00 780. 18.87 +2020-09-19 23:10:00 1500. 18.87 +2020-10-26 12:00:00 1500. 18.87 +2023-05-17 12:00:00 1000. 18.9 + +