diff --git a/environment.yml b/environment.yml index d97b4c72..8852206c 100644 --- a/environment.yml +++ b/environment.yml @@ -6,9 +6,9 @@ dependencies: - astropy=5.2 - python=3.11 # nail the python version, so conda does not try upgrading / dowgrading - ctapipe=0.19 + - protozfits=2.2 - eventio - corsikaio - - protozfits=2 - zeromq - ipython - numba diff --git a/setup.cfg b/setup.cfg index 2d8f6d01..9a633efa 100644 --- a/setup.cfg +++ b/setup.cfg @@ -32,7 +32,7 @@ zip_safe = False install_requires= astropy~=5.2 ctapipe~=0.19.0 - protozfits~=2.0 + protozfits~=2.2 numpy>=1.20 [options.package_data] diff --git a/src/ctapipe_io_lst/__init__.py b/src/ctapipe_io_lst/__init__.py index 11bc7f37..1f814c90 100644 --- a/src/ctapipe_io_lst/__init__.py +++ b/src/ctapipe_io_lst/__init__.py @@ -3,6 +3,7 @@ EventSource for LSTCam protobuf-fits.fz-files. """ from ctapipe.instrument.subarray import EarthLocation +import logging import numpy as np from astropy import units as u from pkg_resources import resource_filename @@ -31,10 +32,10 @@ from ctapipe_io_lst.ground_frame import ground_frame_from_earth_location from .multifiles import MultiFiles -from .containers import LSTArrayEventContainer, LSTServiceContainer, LSTEventContainer +from .containers import LSTArrayEventContainer, LSTServiceContainer, LSTEventContainer, LSTCameraContainer from .version import __version__ from .calibration import LSTR0Corrections -from .event_time import EventTimeCalculator +from .event_time import EventTimeCalculator, cta_high_res_to_time from .pointing import PointingSource from .anyarray_dtypes import ( CDTS_AFTER_37201_DTYPE, @@ -50,6 +51,9 @@ ) +log = logging.getLogger(__name__) + + __all__ = ['LSTEventSource', '__version__'] @@ -278,9 +282,6 @@ def __init__(self, input_url=None, **kwargs): parent=self, ) self.camera_config = self.multi_file.camera_config - self.run_id = self.camera_config.configuration_id - self.tel_id = self.camera_config.telescope_id - self.run_start = Time(self.camera_config.date, format='unix') self.dvr_applied = self.multi_file.dvr_applied reference_location = EarthLocation( @@ -288,18 +289,38 @@ def __init__(self, input_url=None, **kwargs): lat=self.reference_position_lat * u.deg, height=self.reference_position_height * u.m, ) + + self.cta_r1 = self.multi_file.cta_r1 + if self.cta_r1: + self.tel_id = self.camera_config.tel_id + self.local_run_id = self.camera_config.local_run_id + self.module_id_map = self.camera_config.module_id_map + self.pixel_id_map = self.camera_config.pixel_id_map + self.run_start = Time(self.camera_config.config_time_s, format="unix") + self.n_pixels = self.camera_config.num_pixels + self.n_samples = self.camera_config.num_samples_nominal + self.lst_service = self.fill_lst_service_container_ctar1(self.camera_config) + else: + self.tel_id = self.camera_config.telescope_id + self.local_run_id = self.camera_config.configuration_id + self.module_id_map = self.camera_config.lstcam.expected_modules_id + self.pixel_id_map = self.camera_config.expected_pixels_id + self.run_start = Time(self.camera_config.date, format="unix") + self.n_pixels = self.camera_config.num_pixels + self.n_samples = self.camera_config.num_samples + self.lst_service = self.fill_lst_service_container(self.tel_id, self.camera_config) + self._subarray = self.create_subarray(self.tel_id, reference_location) self.r0_r1_calibrator = LSTR0Corrections( subarray=self._subarray, parent=self ) self.time_calculator = EventTimeCalculator( subarray=self.subarray, - run_id=self.run_id, - expected_modules_id=self.camera_config.lstcam.expected_modules_id, + run_id=self.local_run_id, + expected_modules_id=self.module_id_map, parent=self, ) self.pointing_source = PointingSource(subarray=self.subarray, parent=self) - self.lst_service = self.fill_lst_service_container(self.tel_id, self.camera_config) target_info = {} pointing_mode = PointingMode.UNKNOWN @@ -312,17 +333,17 @@ def __init__(self, input_url=None, **kwargs): pointing_mode = PointingMode.TRACK self._scheduling_blocks = { - self.run_id: SchedulingBlockContainer( - sb_id=np.uint64(self.run_id), + self.local_run_id: SchedulingBlockContainer( + sb_id=np.uint64(self.local_run_id), producer_id=f"LST-{self.tel_id}", pointing_mode=pointing_mode, ) } self._observation_blocks = { - self.run_id: ObservationBlockContainer( - obs_id=np.uint64(self.run_id), - sb_id=np.uint64(self.run_id), + self.local_run_id: ObservationBlockContainer( + obs_id=np.uint64(self.local_run_id), + sb_id=np.uint64(self.local_run_id), producer_id=f"LST-{self.tel_id}", actual_start_time=self.run_start, **target_info @@ -414,6 +435,99 @@ def create_subarray(tel_id=1, reference_location=None): return subarray + def fill_from_cta_r1(self, array_event, zfits_event): + tel_id = self.tel_id + + # TODO: decide from applied preprocessing options, not gains + if zfits_event.num_channels == 2: + shape = (zfits_event.num_channels, zfits_event.num_pixels, zfits_event.num_samples) + array_event.r0.tel[self.tel_id] = R0CameraContainer( + waveform=zfits_event.waveform.reshape(shape), + ) + array_event.r1.tel[self.tel_id] = R1CameraContainer() + else: + has_high_gain = (zfits_event.pixel_status & PixelStatus.HIGH_GAIN_STORED).astype(bool) + selected_gain_channel = np.where(has_high_gain, 0, 1) + + shape = (zfits_event.num_pixels, zfits_event.num_samples) + array_event.r1.tel[self.tel_id] = R1CameraContainer( + waveform=zfits_event.waveform.reshape(shape), + selected_gain_channel=selected_gain_channel, + ) + array_event.r0.tel[self.tel_id] = R0CameraContainer() + + array_event.lst.tel[self.tel_id] = self.fill_lst_from_ctar1(zfits_event) + + trigger = array_event.trigger + trigger.time = cta_high_res_to_time(zfits_event.event_time_s, zfits_event.event_time_qns) + trigger.tels_with_trigger = [tel_id] + trigger.tel[tel_id].time = trigger.time + trigger.event_type = self._event_type_from_trigger_bits(zfits_event.event_type) + + def fill_lst_from_ctar1(self, zfits_event): + evt = LSTEventContainer( + pixel_status=zfits_event.pixel_status, + first_capacitor_id=zfits_event.first_cell_id, + calibration_monitoring_id=zfits_event.calibration_monitoring_id, + local_clock_counter=zfits_event.module_hires_local_clock_counter, + ) + + if zfits_event.debug is not None: + debug = zfits_event.debug + evt.module_status = debug.module_status + evt.extdevices_presence = debug.extdevices_presence + evt.chips_flags = debug.chips_flags + evt.charges_hg = debug.charges_gain1 + evt.charges_lg = debug.charges_gain2 + + # unpack Dragon counters + counters = debug.counters.view(DRAGON_COUNTERS_DTYPE) + evt.pps_counter = counters['pps_counter'] + evt.tenMHz_counter = counters['tenMHz_counter'] + evt.event_counter = counters['event_counter'] + evt.trigger_counter = counters['trigger_counter'] + evt.local_clock_counter = counters['local_clock_counter'] + + # if TIB data are there + if evt.extdevices_presence & 1: + tib = debug.tib_data.view(TIB_DTYPE)[0] + evt.tib_event_counter = tib['event_counter'] + evt.tib_pps_counter = tib['pps_counter'] + evt.tib_tenMHz_counter = parse_tib_10MHz_counter(tib['tenMHz_counter']) + evt.tib_stereo_pattern = tib['stereo_pattern'] + evt.tib_masked_trigger = tib['masked_trigger'] + + # if UCTS data are there + if evt.extdevices_presence & 2: + cdts = debug.cdts_data.view(CDTS_AFTER_37201_DTYPE)[0] + evt.ucts_timestamp = cdts["timestamp"] + evt.ucts_address = cdts["address"] + evt.ucts_event_counter = cdts["event_counter"] + evt.ucts_busy_counter = cdts["busy_counter"] + evt.ucts_pps_counter = cdts["pps_counter"] + evt.ucts_clock_counter = cdts["clock_counter"] + evt.ucts_trigger_type = cdts["trigger_type"] + evt.ucts_white_rabbit_status = cdts["white_rabbit_status"] + evt.ucts_stereo_pattern = cdts["stereo_pattern"] + evt.ucts_num_in_bunch = cdts["num_in_bunch"] + evt.ucts_cdts_version = cdts["cdts_version"] + + # if SWAT data are there + if evt.extdevices_presence & 4: + # unpack SWAT data + swat = debug.swat_data.view(SWAT_DTYPE)[0] + evt.swat_assigned_event_id = swat["assigned_event_id"] + evt.swat_trigger_id = swat["trigger_id"] + evt.swat_trigger_type = swat["trigger_type"] + evt.swat_trigger_time_s = swat["trigger_time_s"] + evt.swat_trigger_time_qns = swat["trigger_time_qns"] + evt.swat_readout_requested = swat["readout_requested"] + evt.swat_data_available = swat["data_available"] + evt.swat_hardware_stereo_trigger_mask = swat["hardware_stereo_trigger_mask"] + evt.swat_negative_flag = swat["negative_flag"] + + return LSTCameraContainer(evt=evt, svc=self.lst_service) + def _generator(self): # container for LST data @@ -423,7 +537,6 @@ def _generator(self): array_event.meta['origin'] = 'LSTCAM' # also add service container to the event section - array_event.lst.tel[self.tel_id].svc = self.lst_service # initialize general monitoring container self.initialize_mon_container(array_event) @@ -432,17 +545,23 @@ def _generator(self): for count, zfits_event in enumerate(self.multi_file): array_event.count = count array_event.index.event_id = zfits_event.event_id - array_event.index.obs_id = self.obs_ids[0] + array_event.index.obs_id = self.local_run_id # Skip "empty" events that occur at the end of some runs if zfits_event.event_id == 0: self.log.warning('Event with event_id=0 found, skipping') continue - self.fill_r0r1_container(array_event, zfits_event) - self.fill_lst_event_container(array_event, zfits_event) - if self.trigger_information: - self.fill_trigger_info(array_event) + array_event.lst.tel[self.tel_id].svc = self.lst_service + + if self.cta_r1: + self.fill_from_cta_r1(array_event, zfits_event) + else: + self.fill_r0r1_container(array_event, zfits_event) + self.fill_lst_event_container(array_event, zfits_event) + + if self.trigger_information: + self.fill_trigger_info(array_event) self.fill_mon_container(array_event, zfits_event) @@ -479,6 +598,7 @@ def is_compatible(file_path): try: with fits.open(file_path) as hdul: if "Events" not in hdul: + log.debug("FITS file does not contain an EVENTS HDU, returning False") return False header = hdul["Events"].header @@ -486,19 +606,39 @@ def is_compatible(file_path): value for key, value in header.items() if 'TTYPE' in key } - except OSError: + except Exception as e: + log.debug(f"Error trying to open input file as fits: {e}") return False + if header["XTENSION"] != "BINTABLE": + log.debug(f"EVENTS HDU is not a bintable") + return False - is_protobuf_zfits_file = ( - (header['XTENSION'] == 'BINTABLE') - and (header['ZTABLE'] is True) - and (header['ORIGIN'] == 'CTA') - and (header['PBFHEAD'] in ('R1.CameraEvent', 'ProtoR1.CameraEvent')) - ) + if not header.get("ZTABLE", False): + log.debug(f"ZTABLE is not in header or False") + return False + + if header.get("ORIGIN", "") != "CTA": + log.debug("ORIGIN != CTA") + return False - is_lst_file = 'lstcam_counters' in ttypes - return is_protobuf_zfits_file & is_lst_file + proto_class = header.get("PBFHEAD") + if proto_class is None: + log.debug("Missing PBFHEAD key") + return False + + supported_protos = { + "R1.CameraEvent", + "ProtoR1.CameraEvent", + "CTAR1.Event", + "R1v1.Event", + } + if proto_class not in supported_protos: + log.debug(f"Unsupported PBDHEAD: {proto_class}") + return False + + # TODO: how to differentiate nectarcam from lstcam? + return True @staticmethod def fill_lst_service_container(tel_id, camera_config): @@ -509,6 +649,7 @@ def fill_lst_service_container(tel_id, camera_config): """ return LSTServiceContainer( telescope_id=tel_id, + local_run_id=camera_config.configuration_id, cs_serial=camera_config.cs_serial, configuration_id=camera_config.configuration_id, date=camera_config.date, @@ -524,6 +665,38 @@ def fill_lst_service_container(tel_id, camera_config): pre_proc_algorithms=camera_config.lstcam.pre_proc_algorithms, ) + @staticmethod + def fill_lst_service_container_ctar1(camera_config): + """ + Fill LSTServiceContainer with specific LST service data data + (from the CameraConfig table of zfit file) + + """ + return LSTServiceContainer( + telescope_id=camera_config.tel_id, + local_run_id=camera_config.local_run_id, + date=camera_config.config_time_s, + configuration_id=camera_config.camera_config_id, + pixel_ids=camera_config.pixel_id_map, + module_ids=camera_config.module_id_map, + num_modules=camera_config.num_modules, + num_pixels=camera_config.num_pixels, + num_channels=camera_config.num_channels, + num_samples=camera_config.num_samples_nominal, + + data_model_version=camera_config.data_model_version, + calibration_service_id=camera_config.calibration_service_id, + calibration_algorithm_id=camera_config.calibration_algorithm_id, + + # in debug in CTA R1 debug + cs_serial=camera_config.debug.cs_serial, + idaq_version=camera_config.debug.evb_version, + cdhs_version=camera_config.debug.cdhs_version, + tdp_type=camera_config.debug.tdp_type, + tdp_action=camera_config.debug.tdp_action, + ttype_pattern=camera_config.debug.ttype_pattern, + ) + def fill_lst_event_container(self, array_event, zfits_event): """ Fill LSTEventContainer with specific LST service data @@ -589,19 +762,6 @@ def fill_lst_event_container(self, array_event, zfits_event): lst_evt.ucts_trigger_type = cdts[5] lst_evt.ucts_white_rabbit_status = cdts[6] - # if SWAT data are there - if lst_evt.extdevices_presence & 4: - # unpack SWAT data - unpacked_swat = zfits_event.lstcam.swat_data.view(SWAT_DTYPE)[0] - lst_evt.swat_timestamp = unpacked_swat[0] - lst_evt.swat_counter1 = unpacked_swat[1] - lst_evt.swat_counter2 = unpacked_swat[2] - lst_evt.swat_event_type = unpacked_swat[3] - lst_evt.swat_camera_flag = unpacked_swat[4] - lst_evt.swat_camera_event_num = unpacked_swat[5] - lst_evt.swat_array_flag = unpacked_swat[6] - lst_evt.swat_array_event_num = unpacked_swat[7] - # unpack Dragon counters counters = zfits_event.lstcam.counters.view(DRAGON_COUNTERS_DTYPE) lst_evt.pps_counter = counters['pps_counter'] @@ -767,9 +927,9 @@ def fill_r0r1_camera_container(self, zfits_event): Missing or broken pixels are filled using maxval of the waveform dtype. """ - n_pixels = self.camera_config.num_pixels - n_samples = self.camera_config.num_samples - expected_pixels = self.camera_config.expected_pixels_id + n_pixels = self.n_pixels + n_samples = self.n_samples + pixel_id_map = self.pixel_id_map has_low_gain = (zfits_event.pixel_status & PixelStatus.LOW_GAIN_STORED).astype(bool) has_high_gain = (zfits_event.pixel_status & PixelStatus.HIGH_GAIN_STORED).astype(bool) @@ -807,11 +967,11 @@ def fill_r0r1_camera_container(self, zfits_event): " the run for which this happened." ) - reordered_waveform = np.full((N_PIXELS, N_SAMPLES), fill, dtype=dtype) - reordered_waveform[expected_pixels[stored_pixels]] = waveform + reordered_waveform = np.full((N_PIXELS, n_samples), fill, dtype=dtype) + reordered_waveform[pixel_id_map[stored_pixels]] = waveform reordered_selected_gain = np.full(N_PIXELS, -1, dtype=np.int8) - reordered_selected_gain[expected_pixels] = selected_gain + reordered_selected_gain[pixel_id_map] = selected_gain r0 = R0CameraContainer() r1 = R1CameraContainer( @@ -822,7 +982,7 @@ def fill_r0r1_camera_container(self, zfits_event): reshaped_waveform = zfits_event.waveform.reshape(N_GAINS, -1, n_samples) # re-order the waveform following the expected_pixels_id values reordered_waveform = np.full((N_GAINS, N_PIXELS, N_SAMPLES), fill, dtype=dtype) - reordered_waveform[:, expected_pixels[stored_pixels], :] = reshaped_waveform + reordered_waveform[:, pixel_id_map[stored_pixels], :] = reshaped_waveform r0 = R0CameraContainer(waveform=reordered_waveform) r1 = R1CameraContainer() @@ -864,10 +1024,10 @@ def fill_mon_container(self, array_event, zfits_event): status_container = array_event.mon.tel[self.tel_id].pixel_status # reorder the array - expected_pixels_id = self.camera_config.expected_pixels_id + pixel_id_map = self.pixel_id_map reordered_pixel_status = np.zeros(N_PIXELS, dtype=zfits_event.pixel_status.dtype) - reordered_pixel_status[expected_pixels_id] = zfits_event.pixel_status + reordered_pixel_status[pixel_id_map] = zfits_event.pixel_status channel_info = get_channel_info(reordered_pixel_status) status_container.hardware_failing_pixels[:] = channel_info == 0 diff --git a/src/ctapipe_io_lst/anyarray_dtypes.py b/src/ctapipe_io_lst/anyarray_dtypes.py index 7ec96d3f..d7762c4e 100644 --- a/src/ctapipe_io_lst/anyarray_dtypes.py +++ b/src/ctapipe_io_lst/anyarray_dtypes.py @@ -50,16 +50,18 @@ ('unknown', np.uint8), # called arbitraryInformation in C-Struct ]).newbyteorder('<') + SWAT_DTYPE = np.dtype([ - ('timestamp', np.uint64), - ('counter1', np.uint32), - ('counter2', np.uint32), - ('event_type', np.uint8), - ('camera_flag', np.uint8), - ('camera_event_num', np.uint32), - ('array_flag', np.uint8), - ('event_num', np.uint32), -]).newbyteorder('<') + ("assigned_event_id", np.uint32), + ("trigger_id", np.uint64), + ("trigger_type", np.uint8), + ("trigger_time_s", np.uint32), + ("trigger_time_qns", np.uint32), + ("readout_requested", np.bool_), + ("data_available", np.bool_), + ("hardware_stereo_trigger_mask", np.uint16), + ("negative_flag", np.uint8), +], align=True).newbyteorder('<') def parse_tib_10MHz_counter(counter): diff --git a/src/ctapipe_io_lst/containers.py b/src/ctapipe_io_lst/containers.py index 87e4c16f..e863e085 100644 --- a/src/ctapipe_io_lst/containers.py +++ b/src/ctapipe_io_lst/containers.py @@ -1,6 +1,7 @@ """ Container structures for data that should be read or written to disk """ +import numpy as np from ctapipe.core import Container, Field, Map from ctapipe.containers import ArrayEventContainer from functools import partial @@ -21,41 +22,58 @@ class LSTServiceContainer(Container): """ # Data from the CameraConfig table + # present in CTA R1 telescope_id = Field(-1, "telescope id") - cs_serial = Field(None, "serial number of the camera server") - configuration_id = Field(None, "id of the CameraConfiguration") - date = Field(None, "NTP start of run date") + local_run_id = Field(-1, "local run id") + date = Field(None, "config_time_s in data model") + configuration_id = Field(None, "camera_config_id in the date model") + pixel_ids = Field([], "pixel_id_map in the data model") + module_ids = Field([], "module_id_map in the data model") + + num_modules = Field(-1, "number of modules") num_pixels = Field(-1, "number of pixels") + num_channels = Field(-1, "number of channels") num_samples = Field(-1, "num samples") - pixel_ids = Field([], "id of the pixels in the waveform array") + data_model_version = Field(None, "data model version") + calibration_service_id = Field(-1, "calibration service id") + calibration_algorithm_id = Field(-1, "calibration service id") - idaq_version = Field(0, "idaq version") + # in debug in CTA R1 debug + cs_serial = Field(None, "serial number of the camera server") + idaq_version = Field(0, "idaq/evb version") cdhs_version = Field(0, "cdhs version") + tdp_type = Field(None, "tdp type") + tdp_action = Field(None, "tdp action") + ttype_pattern = Field(None, "ttype_pattern") + + # only in old R1 algorithms = Field(None, "algorithms") pre_proc_algorithms = Field(None, "pre processing algorithms") - module_ids = Field([], "module ids") - num_modules = Field(-1, "number of modules") class LSTEventContainer(Container): """ Container for Fields that are specific to each LST event """ + # present in CTA R1 but not in ctapipe R1CameraEvent + pixel_status = Field(None, "status of the pixels (n_pixels)", dtype=np.uint8) + first_capacitor_id = Field(None, "first capacitor id") + calibration_monitoring_id = Field(None, "calibration id of applied pre-calibration") + local_clock_counter = Field(None, "Dragon local 133 MHz counter (n_modules)") - # Data from the CameraEvent table - configuration_id = Field(None, "id of the CameraConfiguration") - event_id = Field(None, "local id of the event") - tel_event_id = Field(None, "global id of the event") - pixel_status = Field([], "status of the pixels (n_pixels)") - ped_id = Field(None, "tel_event_id of the event used for pedestal substraction") - module_status = Field([], "status of the modules (n_modules)") + # in debug event + module_status = Field(None, "status of the modules (n_modules)") extdevices_presence = Field(None, "presence of data for external devices") - - tib_event_counter = Field(-1, "TIB event counter") - tib_pps_counter = Field(-1, "TIB pps counter") - tib_tenMHz_counter = Field(-1, "TIB 10 MHz counter") - tib_stereo_pattern = Field(0, "TIB stereo pattern") + chips_flags = Field(None, "chips flags") + charges_hg = Field(None, "charges of high gain channel") + charges_lg = Field(None, "charges of low gain channel") + tdp_action = Field(None, "tdp action") + + tib_event_counter = Field(np.uint32(0), "TIB event counter", dtype=np.uint32) + tib_pps_counter = Field(np.uint16(0), "TIB pps counter", dtype=np.uint16) + tib_tenMHz_counter = Field(np.uint32(0), "TIB 10 MHz counter", dtype=np.uint32) + tib_stereo_pattern = Field(np.uint16(0), "TIB stereo pattern", dtype=np.uint16) tib_masked_trigger = Field(0, "TIB trigger mask") ucts_event_counter = Field(-1, "UCTS event counter") @@ -71,26 +89,31 @@ class LSTEventContainer(Container): ucts_num_in_bunch = Field(-1, "UCTS num in bunch (for debugging)") ucts_cdts_version = Field(-1, "UCTS CDTS version") - swat_timestamp = Field(-1, "SWAT timestamp") - swat_counter1 = Field(-1, "SWAT event counter 1") - swat_counter2 = Field(-1, "SWAT event counter 2") - swat_event_type = Field(0, "SWAT event type") - swat_camera_flag = Field(-1, "SWAT camera flag ") - swat_camera_event_num = Field(-1, "SWAT camera event number") - swat_array_flag = Field(-1, "SWAT array negative flag") - swat_array_event_num = Field(-1, "SWAT array event number") + swat_assigned_event_id = Field(np.uint32(0), "SWAT assigned event id") + swat_trigger_id = Field(np.uint64(0), "SWAT trigger id") + swat_trigger_type = Field(np.uint8(0), "SWAT trigger type") + swat_trigger_time_s = Field(np.uint32(0), "SWAT trigger_time_s") + swat_trigger_time_qns = Field(np.uint32(0), "SWAT trigger_time_qns") + swat_readout_requested = Field(np.bool_(False), "SWAT readout requested") + swat_data_available = Field(np.bool_(False), "SWAT data available") + swat_hardware_stereo_trigger_mask = Field(np.uint16(0), "SWAT hardware stereo trigger mask") + swat_negative_flag = Field(np.uint8(0), "SWAT negative flag") pps_counter= Field(None, "Dragon pulse per second counter (n_modules)") tenMHz_counter = Field(None, "Dragon 10 MHz counter (n_modules)") event_counter = Field(None, "Dragon event counter (n_modules)") trigger_counter = Field(None, "Dragon trigger counter (n_modules)") - local_clock_counter = Field(None, "Dragon local 133 MHz counter (n_modules)") - chips_flags = Field(None, "chips flags") - first_capacitor_id = Field(None, "first capacitor id") + # Only in old R1 + configuration_id = Field(None, "id of the CameraConfiguration") + event_id = Field(None, "global id of the event") + tel_event_id = Field(None, "local id of the event") + ped_id = Field(None, "tel_event_id of the event used for pedestal substraction") + drs_tag_status = Field(None, "DRS tag status") drs_tag = Field(None, "DRS tag") + # custom here ucts_jump = Field(False, "A ucts jump happened in the current event") diff --git a/src/ctapipe_io_lst/event_time.py b/src/ctapipe_io_lst/event_time.py index e4570783..df0c1492 100644 --- a/src/ctapipe_io_lst/event_time.py +++ b/src/ctapipe_io_lst/event_time.py @@ -19,7 +19,57 @@ TimeUnixTai.epoch_scale = 'tai' +EPOCH = Time(0, format='unix_tai', scale='tai') +DAY_TO_S = 86400 CENTRAL_MODULE = 132 +S_TO_QNS = 4e9 +QNS_TO_S = 0.25e-9 + + +def cta_high_res_to_time(seconds, quarter_nanoseconds): + '''Convert cta high resolution timestamp to astropy Time''' + # unix_tai accepts two floats for maximum precision + # we can just pass integral and fractional part + fractional_seconds = quarter_nanoseconds * QNS_TO_S + return Time( + seconds, + fractional_seconds, + format='unix_tai', + # this is only for displaying iso timestamp, not any actual precision + precision=9, + ) + + +def to_seconds(days): + '''Returns whole and fractional seconds from a number in days''' + seconds = days * DAY_TO_S + return np.divmod(seconds, 1) + + +def time_to_cta_high(time): + '''Convert astropy Time to cta high precision timestamp''' + # make sure we are in TAI + time = time.tai + + # internally, astropy always uses jd values + # jd1 is integral and jd2 is in [-0.5, 0.5] + # we get the integral and fractional seconds from both jd values + # relative to the epoch + seconds_jd1, fractional_seconds_jd1 = to_seconds(time.jd1 - EPOCH.jd1) + seconds_jd2, fractional_seconds_jd2 = to_seconds(time.jd2 - EPOCH.jd2) + + # add up the integral number of seconds + seconds = seconds_jd1 + seconds_jd2 + + # convert fractional seconds to quarter nanoseconds + fractional_seconds = fractional_seconds_jd1 + fractional_seconds_jd2 + quarter_nanoseconds = np.round(fractional_seconds * S_TO_QNS) + + # convert to uint32 + seconds = np.asanyarray(seconds, dtype=np.uint32) + quarter_nanoseconds = np.asanyarray(quarter_nanoseconds, dtype=np.uint32) + return seconds, quarter_nanoseconds + # fix for https://github.com/ipython/traitlets/issues/637 diff --git a/src/ctapipe_io_lst/multifiles.py b/src/ctapipe_io_lst/multifiles.py index 69003be3..91bb11b0 100644 --- a/src/ctapipe_io_lst/multifiles.py +++ b/src/ctapipe_io_lst/multifiles.py @@ -1,5 +1,4 @@ import re -import warnings from pathlib import Path from collections import namedtuple, defaultdict from itertools import count @@ -160,9 +159,10 @@ def _load_next_subrun(self, stream): self._files.pop(stream).close() Provenance().add_input_file(str(path), "R0") - self._files[stream] = File(str(path)) + file_ = File(str(path)) + self._files[stream] = file_ self.log.info("Opened file %s", path) - self._events_tables[stream] = self._files[stream].Events + self._events_tables[stream] = file_.Events self._headers[stream] = self._events_tables[stream].header dvr_applied = self._headers[stream].get("LSTDVR", False) if self.dvr_applied is None: @@ -175,7 +175,13 @@ def _load_next_subrun(self, stream): self._events.put_nowait(NextEvent(event.event_id, event, stream)) # make sure we have a camera config - config = next(self._files[stream].CameraConfig) + if hasattr(file_, "CameraConfig"): + config = next(file_.CameraConfig) + self.cta_r1 = False + else: + # new files use CameraConfiguration + self.cta_r1 = True + config = next(file_.CameraConfiguration) if self.camera_config is None: self.camera_config = config diff --git a/src/ctapipe_io_lst/tests/test_cta_r1.py b/src/ctapipe_io_lst/tests/test_cta_r1.py new file mode 100644 index 00000000..7ebc4d29 --- /dev/null +++ b/src/ctapipe_io_lst/tests/test_cta_r1.py @@ -0,0 +1,285 @@ +import os +from pathlib import Path +from contextlib import ExitStack +import pytest +import numpy as np +from ctapipe.io import EventSource +from astropy.time import Time +import astropy.units as u + +import protozfits +from protozfits.CTA_R1_pb2 import CameraConfiguration, Event +from protozfits.Debug_R1_pb2 import DebugEvent, DebugCameraConfiguration +from protozfits.CoreMessages_pb2 import AnyArray +from traitlets.config import Config +from ctapipe_io_lst import LSTEventSource +from ctapipe_io_lst.anyarray_dtypes import CDTS_AFTER_37201_DTYPE, TIB_DTYPE +from ctapipe_io_lst.constants import CLOCK_FREQUENCY_KHZ +from ctapipe.image.toymodel import WaveformModel, Gaussian +import socket + +from ctapipe_io_lst.event_time import time_to_cta_high + + +test_data = Path(os.getenv('LSTCHAIN_TEST_DATA', 'test_data')) +test_drs4_pedestal_path = test_data / 'real/monitoring/PixelCalibration/LevelA/drs4_baseline/20200218/v0.8.2.post2.dev48+gb1343281/drs4_pedestal.Run02005.0000.h5' + + +ANY_ARRAY_TYPE_TO_NUMPY_TYPE = { + 1: np.int8, + 2: np.uint8, + 3: np.int16, + 4: np.uint16, + 5: np.int32, + 6: np.uint32, + 7: np.int64, + 8: np.uint64, + 9: np.float32, + 10: np.float64, +} + +DTYPE_TO_ANYARRAY_TYPE = {v: k for k, v in ANY_ARRAY_TYPE_TO_NUMPY_TYPE.items()} + + +subarray = LSTEventSource.create_subarray(tel_id=1) +GEOMETRY = subarray.tel[1].camera.geometry +waveform_model = WaveformModel( + reference_pulse=subarray.tel[1].camera.readout.reference_pulse_shape[0], + reference_pulse_sample_width=subarray.tel[1].camera.readout.reference_pulse_sample_width, + sample_width=1 * u.ns, +) + + +def weighted_average(a, b, weight_a, weight_b): + return (a * weight_a + b * weight_b) / (weight_a + weight_b) + + +def create_shower(rng): + x = rng.uniform(-0.5, 0.5) * u.m + y = rng.uniform(-0.5, 0.5) * u.m + width = rng.uniform(0.01, 0.05) * u.m + length = rng.uniform(2 * width.to_value(u.m), 5 * width.to_value(u.m)) * u.m + psi = rng.uniform(0, 360) * u.deg + + + area = np.pi * (width.to_value(u.m) * length.to_value(u.m)) + intensity = 5e4 * area + 2e6 * area**2 + + model = Gaussian(x, y, length, width, psi) + image, signal, noise = model.generate_image(GEOMETRY, intensity=intensity, nsb_level_pe=3) + + long = (GEOMETRY.pix_x - x) * np.cos(psi) + (GEOMETRY.pix_y - y) * np.sin(psi) + + peak_time = rng.uniform(0, 40, size=len(GEOMETRY)) + signal_peak_time = np.clip(20 + 30 * (long / (3 * u.m)).to_value(u.one), 0, 40) + + mask = signal > 0 + peak_time[mask] = weighted_average( + peak_time[mask], signal_peak_time[mask], noise[mask], signal[mask] + ) + + return image, peak_time + + +def create_waveform(image, peak_time, num_samples=40, gains=(86, 5), offset=400): + r1 = waveform_model.get_waveform(image, peak_time, num_samples) + return np.array([r1 * gain + offset for gain in gains]).astype(np.uint16) + + +def create_flat_field(rng): + image = rng.uniform(65, 75, len(GEOMETRY)) + peak_time = rng.uniform(18, 22, len(GEOMETRY)) + return image, peak_time + +def create_pedestal(rng): + image = rng.uniform(-2, 2, len(GEOMETRY)) + peak_time = rng.uniform(0, 40, len(GEOMETRY)) + return image, peak_time + + +def to_anyarray(array): + type_ = DTYPE_TO_ANYARRAY_TYPE[array.dtype.type] + return AnyArray(type=type_, data=array.tobytes()) + + +@pytest.fixture(scope="session") +def dummy_cta_r1_dir(tmp_path_factory): + return tmp_path_factory.mktemp("dummy_cta_r1") + + +@pytest.fixture(scope="session") +def dummy_cta_r1(dummy_cta_r1_dir): + # with protozfits.File("test_data/real/R0/20200218/LST-1.1.Run02006.0004.fits.fz") as f: + # old_camera_config = f.CameraConfig[0] + + stream_paths = [ + dummy_cta_r1_dir / "LST-1.1.Run10000.0000.fits.fz", + dummy_cta_r1_dir / "LST-1.2.Run10000.0000.fits.fz", + dummy_cta_r1_dir / "LST-1.3.Run10000.0000.fits.fz", + dummy_cta_r1_dir / "LST-1.4.Run10000.0000.fits.fz", + ] + + num_samples = 40 + num_pixels = 1855 + num_modules = 265 + ucts_address = np.frombuffer(socket.inet_pton(socket.AF_INET, "10.0.100.123"), dtype=np.uint32)[0] + + run_start = Time("2023-05-16T16:06:31.123") + camera_config = CameraConfiguration( + tel_id=1, + local_run_id=10000, + config_time_s=run_start.unix, + camera_config_id=1, + num_modules=num_modules, + num_pixels=num_pixels, + num_channels=2, + data_model_version="1.0", + num_samples_nominal=num_samples, + pixel_id_map=to_anyarray(np.arange(num_pixels).astype(np.uint16)), + module_id_map=to_anyarray(np.arange(num_modules).astype(np.uint16)), + debug=DebugCameraConfiguration( + cs_serial="???", + evb_version="evb-dummy", + cdhs_version="evb-dummy", + ) + ) + + rng = np.random.default_rng() + + streams = [] + with ExitStack() as stack: + for stream_path in stream_paths: + stream = stack.enter_context(protozfits.ProtobufZOFits(n_tiles=5, rows_per_tile=20, compression_block_size_kb=64 * 1024)) + stream.open(str(stream_path)) + stream.move_to_new_table("CameraConfiguration") + stream.write_message(camera_config) + stream.move_to_new_table("Events") + streams.append(stream) + + event_time = run_start + event_rate = 8000 / u.s + + module_hires_local_clock_counter = np.zeros(num_modules, dtype=np.uint64) + for event_count in range(100): + if event_count % 20 == 18: + # flatfield + event_type = 0 + trigger_type = 0b0000_0001 + image, peak_time = create_flat_field(rng) + elif event_count % 20 == 19: + # pedestal + event_type = 2 + trigger_type = 0b0100_0000 + image, peak_time = create_pedestal(rng) + else: + # air shower + event_type = 32 + trigger_type = 0b0000_0100 + image, peak_time = create_shower(rng) + + delta = rng.exponential(1 / event_rate.to_value(1 / u.s)) * u.s + event_time = event_time + delta + event_time_s, event_time_qns = time_to_cta_high(event_time) + + waveform = create_waveform(image, peak_time, num_samples) + + if event_type == 32: + num_channels = 1 + pixel_status = rng.choice([0b1000, 0b0100], p=[0.001, 0.999], size=num_pixels).astype(np.uint8) + + low_gain = pixel_status == 0b1000 + gain_selected_waveform = waveform[0].copy() + gain_selected_waveform[low_gain] = waveform[1, low_gain] + waveform = gain_selected_waveform + else: + num_channels = 2 + pixel_status = np.full(num_pixels, 0b1100, dtype=np.uint8) + + first_cell_id = rng.choice(4096, size=2120).astype(np.uint8) + + jitter = rng.choice([-2, 1, 0, 1, -2], size=num_modules) + module_hires_local_clock_counter[:] += np.uint64(CLOCK_FREQUENCY_KHZ * delta.to_value(u.ms) + jitter) + + # uint64 unix_tai, whole ns + timestamp = np.uint64(event_time_s) * np.uint64(1e9) + np.uint64(np.round(event_time_qns / 4.0)) + cdts_data = ( + timestamp, + ucts_address, + event_count + 1, # event_counter + 0, # busy_counter + 0, # pps_counter + 0, # clock_counter + trigger_type, # trigger_type + 0, # white_rabbit_status + 0, # stereo_pattern + event_count % 10, # num_in_bunch + 1234, # cdts_version + ) + cdts_data = np.array([cdts_data], dtype=CDTS_AFTER_37201_DTYPE) + + tib_data = ( + event_count + 1, + 0, # pps_counter + 0, # tenMHz_counter + 0, # stereo pattern + trigger_type, + ) + tib_data = np.array([tib_data], dtype=TIB_DTYPE) + + event = Event( + event_id=event_count + 1, + tel_id=1, + local_run_id=10000, + event_type=event_type, + event_time_s=int(event_time_s), + event_time_qns=int(event_time_qns), + num_channels=num_channels, + num_pixels=num_pixels, + num_samples=num_samples, + pixel_status=to_anyarray(pixel_status), + waveform=to_anyarray(waveform), + first_cell_id=to_anyarray(first_cell_id), + module_hires_local_clock_counter=to_anyarray(module_hires_local_clock_counter), + debug=DebugEvent( + module_status=to_anyarray(np.ones(num_modules, dtype=np.uint8)), + extdevices_presence=0b011, + cdts_data=to_anyarray(cdts_data.view(np.uint8)), + tib_data=to_anyarray(tib_data.view(np.uint8)), + ) + ) + + + streams[event_count % len(streams)].write_message(event) + + return stream_paths[0] + + +def test_is_compatible(dummy_cta_r1): + from ctapipe_io_lst import LSTEventSource + assert LSTEventSource.is_compatible(dummy_cta_r1) + + +def test_no_calibration(dummy_cta_r1): + with EventSource(dummy_cta_r1, apply_drs4_corrections=False, pointing_information=False) as source: + n_events = 0 + for e in source: + n_events += 1 + assert n_events == 100 + + +def test_drs4_calibration(dummy_cta_r1): + config = Config({ + 'LSTEventSource': { + 'pointing_information': False, + 'LSTR0Corrections': { + 'drs4_pedestal_path': test_drs4_pedestal_path, + }, + }, + }) + + with EventSource(dummy_cta_r1, config=config) as source: + + n_events = 0 + for e in source: + n_events += 1 + assert n_events == 100