diff --git a/.codacy.yml b/.codacy.yml index 260dd03..55a070d 100644 --- a/.codacy.yml +++ b/.codacy.yml @@ -6,3 +6,5 @@ engines: pyflakes: disable: - F999 + assert_used: skips: ['*/*/test_*.py', '*/test_*.py', 'test_*.py'] + diff --git a/ctapipe_io_magic/__init__.py b/ctapipe_io_magic/__init__.py index 4b73773..362ae95 100644 --- a/ctapipe_io_magic/__init__.py +++ b/ctapipe_io_magic/__init__.py @@ -10,6 +10,7 @@ import uproot import logging import numpy as np +from typing import List, Any from pathlib import Path from decimal import Decimal from astropy import units as u @@ -17,7 +18,7 @@ from pkg_resources import resource_filename from ctapipe.io import EventSource, DataLevel -from ctapipe.core import Provenance +from ctapipe.core import Provenance, Container, Field from ctapipe.core.traits import Bool, UseEnum from ctapipe.coordinates import CameraFrame @@ -92,6 +93,68 @@ DATA_MAGIC_LST_TRIGGER: EventType.SUBARRAY, } +class ReportLaserContainer(Container): + """ Container for Magic laser parameters """ + UniqueID = Field(Any, 'No.') + Bits = Field(Any, 'ID') + MJD = Field(np.float64, 'Modified Julian Date') + BadReport = Field(Any, 'Bad Report') + State = Field(Any, 'State') + IsOffsetCorrection = Field(Any, 'Is Offset Correction') + IsOffsetFitted = Field(Any, 'Is Offset Fitted') + IsBGCorrection = Field(Any, 'Is BG Correction') + IsT0ShiftFitted = Field(Any, 'Is T0 Shift Fitted') + IsUseGDAS = Field(Any, 'Is Use GDAS') + IsUpwardMoving = Field(Any, 'Is Upward Moving') + OverShoot = Field(Any, 'Over Shoot') + UnderShoot = Field(Any, 'Under Shoot') + BGSamples = Field(Any, 'BG Samples') + Transmission3km = Field(Any, 'Transmission at 3 km') + Transmission6km = Field(Any, 'Transmission at 6 km') + Transmission9km = Field(Any, 'Transmission at 9 km') + Transmission12km = Field(Any, 'Transmission at 12 km') + Zenith = Field(Any, 'Zenith angle', unit=u.deg) + Azimuth = Field(Any, 'Azimuth angle', unit=u.deg) + FullOverlap = Field(Any, 'Full Overlap') + EndGroundLayer = Field(Any, 'End Ground Layer') + GroundLayerTrans = Field(Any, 'Ground Layer Trans') + Calimaness = Field(Any, 'Calimaness') + CloudLayerAlt = Field(Any, 'Altitude of cloud layer') + CloudLayerDens = Field(Any, 'Density of cloud layer') + Klett_k = Field(Any, 'Klett k') + PheCounts = Field(List[int], 'Phe Counts') + Offset = Field(Any, 'Offset') + Offset_Calculated = Field(Any, 'Offset calculated') + Offset_Fitted = Field(Any, 'Offset fitted') + Offset2 = Field(Any, 'Offset 2') + Offset3 = Field(Any, 'Offset 3') + Background1 = Field(Any, 'Background 1') + Background2 = Field(Any, 'Background 2') + BackgroundErr1 = Field(Any, 'Background error 1') + BackgroundErr2 = Field(Any, 'Background error 2') + RangeMax = Field(Any, 'Range max') + RangeMax_Clouds = Field(Any, 'Range max clouds') + ErrorCode = Field(Any, 'Error code') + ScaleHeight_fit = Field(Any, 'Scale Height fit') + Alpha_fit = Field(Any, 'Alpha fit') + Chi2Alpha_fit = Field(Any, 'Chi2 Alpha fit') + Alpha_firstY = Field(Any, 'Alpha first Y') + Alpha_Junge = Field(Any, 'Alpha Junge') + PBLHeight = Field(Any, 'PBL Height') + Chi2Full_fit = Field(Any, 'Chi2 Full fit') + SignalSamples = Field(Any, 'Signal Samples') + HWSwitch = Field(Any, 'HW Switch') + HWSwitchMaxOffset = Field(Any, 'HW Switch Max Offset') + NCollapse = Field(Any, 'N Collapse') + Shots = Field(Any, 'Shots') + T0Shift = Field(Any, 'T0 Shift') + Interval_0 = Field(Any, 'Interval 0') + RCS_min_perfect = Field(Any, 'RCS min perfect') + RCS_min_clouds = Field(Any, 'RCS min cloud') + RCS_min_mol = Field(Any, 'RCS min mol') + LIDAR_ratio = Field(Any, 'LIDAR ratio') + LIDAR_ratio_Cloud = Field(Any, 'LIDAR ratio cloud') + LIDAR_ratio_Junge = Field(Any, 'LIDAR ratio Junge') def load_camera_geometry(): """Load camera geometry from bundled resources of this repo""" @@ -223,6 +286,9 @@ def __init__(self, input_url=None, config=None, parent=None, **kwargs): self.metadata = self.parse_metadata_info() + if not self.is_simulation: + self.laser = self.parse_laser_info() + # Retrieving the data level (so far HARDCODED Sorcerer) self.datalevel = DataLevel.DL0 @@ -381,13 +447,13 @@ def get_run_info_from_name(file_name): elif re.match(mask_data_superstar, file_name) is not None: parsed_info = re.match(mask_data_superstar, file_name) telescope = None - run_number = int(parsed_info.grou(1)) + run_number = int(parsed_info.group(1)) datalevel = MARSDataLevel.SUPERSTAR is_mc = False elif re.match(mask_data_melibea, file_name) is not None: parsed_info = re.match(mask_data_melibea, file_name) telescope = None - run_number = int(parsed_info.grou(1)) + run_number = int(parsed_info.group(1)) datalevel = MARSDataLevel.MELIBEA is_mc = False elif re.match(mask_mc_calibrated, file_name) is not None: @@ -1035,6 +1101,150 @@ def parse_metadata_info(self): return metadata + def parse_laser_info(self): + laser_info_array_list_runh = [ + 'MReportLaser.MReport.fUniqueID', + 'MReportLaser.MReport.fBits', + 'MTimeLaser.fMjd', + 'MTimeLaser.fTime.fMilliSec', + 'MReportLaser.MReport.fBadReport', + 'MReportLaser.MReport.fState', + 'MReportLaser.fIsOffsetCorrection', + 'MReportLaser.fIsOffsetFitted', + 'MReportLaser.fIsBGCorrection', + 'MReportLaser.fIsT0ShiftFitted', + 'MReportLaser.fIsUseGDAS', + 'MReportLaser.fIsUpwardMoving', + 'MReportLaser.fOverShoot', + 'MReportLaser.fUnderShoot', + 'MReportLaser.fBGSamples', + 'MReportLaser.fTransmission3km', + 'MReportLaser.fTransmission6km', + 'MReportLaser.fTransmission9km', + 'MReportLaser.fTransmission12km', + 'MReportLaser.fZenith', + 'MReportLaser.fAzimuth', + 'MReportLaser.fFullOverlap', + 'MReportLaser.fEndGroundLayer', + 'MReportLaser.fGroundLayerTrans', + 'MReportLaser.fKlett_k', + 'MReportLaser.fPheCounts', + 'MReportLaser.fCalimaness', + 'MReportLaser.fCloudLayerAlt', + 'MReportLaser.fCloudLayerDens', + 'MReportLaser.fOffset', + 'MReportLaser.fOffset_Calculated', + 'MReportLaser.fOffset_Fitted', + 'MReportLaser.fOffset2', + 'MReportLaser.fOffset3', + 'MReportLaser.fBackground1', + 'MReportLaser.fBackground2', + 'MReportLaser.fBackgroundErr1', + 'MReportLaser.fBackgroundErr2', + 'MReportLaser.fRangeMax', + 'MReportLaser.fRangeMax_Clouds', + 'MReportLaser.fErrorCode', + 'MReportLaser.fScaleHeight_fit', + 'MReportLaser.fChi2Alpha_fit', + 'MReportLaser.fChi2Alpha_fit', + 'MReportLaser.fAlpha_firstY', + 'MReportLaser.fAlpha_Junge', + 'MReportLaser.fPBLHeight', + 'MReportLaser.fChi2Full_fit', + 'MReportLaser.fHWSwitchMaxOffset', + 'MReportLaser.fNCollapse', + 'MReportLaser.fShots', + 'MReportLaser.fT0Shift', + 'MReportLaser.fInterval_0', + 'MReportLaser.fRCS_min_perfect', + 'MReportLaser.fRCS_min_clouds', + 'MReportLaser.fRCS_min_mol', + 'MReportLaser.fLIDAR_ratio', + 'MReportLaser.fLIDAR_ratio_Cloud', + 'MReportLaser.fLIDAR_ratio_Junge', + ] + + old_mjd = None + unique_reports = [] + for rootf in self.files_: + try: + laser_info_runh = rootf['Laser'].arrays( + laser_info_array_list_runh, library="np") + mjd_value = laser_info_runh['MTimeLaser.fMjd'] + millisec_value = laser_info_runh['MTimeLaser.fTime.fMilliSec'] + 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] + + except KeyError as e: + print(f"Required key not found in the file {rootf}: {e}") + continue + return unique_reports + def parse_simulation_header(self): """ Parse the simulation information from the RunHeaders tree. 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_params_corr_test.py b/ctapipe_io_magic/laser_params_corr_test.py new file mode 100644 index 0000000..9c574d8 --- /dev/null +++ b/ctapipe_io_magic/laser_params_corr_test.py @@ -0,0 +1,411 @@ +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 + +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): + 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 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 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): + 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 + + 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): + 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): + 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 + 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): + 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 * 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 * gkIntegrationWindow + if ifadc < 384: + return r + r += 128 * gkIntegrationWindow * 2 + if ifadc < 448: + return r + r += 64 * gkIntegrationWindow * 4 + return r + + def get_fadcr(self, i, bgsamples): + 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: + return ifadc % 128 + 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(gkIntegrationWindow * i, 2), 2) + background1 /= 30.0 + bg = background1 + for i in range(signalsamples): + 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) + 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 + 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() + 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/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..4830fd6 --- /dev/null +++ b/ctapipe_io_magic/switch_times.txt @@ -0,0 +1,35 @@ +# 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 + + diff --git a/ctapipe_io_magic/tests/test_lidar.py b/ctapipe_io_magic/tests/test_lidar.py new file mode 100644 index 0000000..bb920fd --- /dev/null +++ b/ctapipe_io_magic/tests/test_lidar.py @@ -0,0 +1,28 @@ +import os +from pathlib import Path + +import pytest + +test_data = Path(os.getenv("MAGIC_TEST_DATA", "test_data")).absolute() +test_calibrated_real_dir = test_data / "real/calibrated" +test_calibrated_real = test_calibrated_real_dir / "20210314_M1_05095172.002_Y_CrabNebula-W0.40+035.root" + +@pytest.fixture +def sample_report_laser(): + from ctapipe_io_magic import MAGICEventSource + with MAGICEventSource(input_url=test_calibrated_real,process_run=True) as source: + laser = source.laser + return laser + +def test_lidar_parameters(sample_report_laser): + '''Comparing the read Lidar report parameters with hardcoded values''' + laser = sample_report_laser + assert len(laser) == 1 + for report in laser: + 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) 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 +