diff --git a/CHECLabPy/core/io/simtel.py b/CHECLabPy/core/io/simtel.py index fdc695b..404b69d 100644 --- a/CHECLabPy/core/io/simtel.py +++ b/CHECLabPy/core/io/simtel.py @@ -48,7 +48,7 @@ def __init__(self, path, max_events=None): ) self.seeker = EventSeeker(reader) - first_event = self.seeker[0] + first_event = self.seeker.get_event_index(0) tels = list(first_event.r0.tels_with_data) self.tel = tels[0] shape = first_event.r0.tel[self.tel].waveform.shape @@ -62,7 +62,7 @@ def __init__(self, path, max_events=None): self.mapping = get_clp_mapping_from_tc_mapping(tc_mapping) n_rows = self.mapping.metadata['n_rows'] n_columns = self.mapping.metadata['n_columns'] - camera_geom = first_event.inst.subarray.tel[tels[0]].camera + camera_geom = reader.subarray.tel[tels[0]].camera.geometry engineering_frame = EngineeringCameraFrame(n_mirrors=2) engineering_geom = camera_geom.transform_to(engineering_frame) pix_x = engineering_geom.pix_x.value @@ -84,7 +84,7 @@ def __init__(self, path, max_events=None): def _build_waveform(self, event): self._fill_event_containers(event) samples = event.r1.tel[self.tel].waveform[self.pixel_order] - mc_true = event.mc.tel[self.tel].photo_electron_image[self.pixel_order] + mc_true = event.mc.tel[self.tel].true_image[self.pixel_order] waveform = SimtelWaveform( samples, iev=self._iev, @@ -95,7 +95,7 @@ def _build_waveform(self, event): return waveform def _get_event(self, iev): - event = self.seeker[iev] + event = self.seeker.get_event_index(iev) return self._build_waveform(event) def __iter__(self): @@ -119,8 +119,8 @@ def is_compatible(path): def _fill_event_containers(self, event): self._iev = event.count - self._t_cpu = pd.to_datetime(event.trig.gps_time.value, unit='s') - self.run_id = event.r0.obs_id + self._t_cpu = pd.to_datetime(event.trigger.tel[self.tel].time.value, unit='s') + self.run_id = event.index.obs_id self.mc = dict( iev=self._iev, diff --git a/CHECLabPy/utils/mapping.py b/CHECLabPy/utils/mapping.py index a0dfd14..b013ff2 100644 --- a/CHECLabPy/utils/mapping.py +++ b/CHECLabPy/utils/mapping.py @@ -1,5 +1,8 @@ import warnings import numpy as np +from ctapipe.instrument import SubarrayDescription, TelescopeDescription, \ + CameraDescription, CameraReadout, OpticsDescription +from astropy import units as u def get_clp_mapping_from_tc_mapping(tc_mapping): @@ -169,11 +172,25 @@ def get_ctapipe_camera_geometry(mapping): pix_area=None, pix_type='rectangular', frame=EngineeringCameraFrame(n_mirrors=2), - sampling_rate=u.Quantity(1, u.GHz), ) return camera +def get_ctapipe_subarray(mapping): + tel_id = 0 + geom = get_ctapipe_camera_geometry(mapping) + camera = CameraDescription("CHEC", geom, CameraReadout.from_name("CHEC")) + optics = OpticsDescription.from_name("SST-ASTRI") + subarray = SubarrayDescription( + "test array", + tel_positions={tel_id: np.zeros(3) * u.m}, + tel_descriptions={ + tel_id: TelescopeDescription("SST-ASTRI_CHEC", "SST", optics, camera) + }, + ) + return subarray, tel_id + + def get_row_column(pix_x, pix_y): """ Get the row and column of rectangular pixels of a camera that are diff --git a/CHECLabPy/waveform_reducers/ctapipe_integrators.py b/CHECLabPy/waveform_reducers/ctapipe_integrators.py index 1d14832..5734309 100644 --- a/CHECLabPy/waveform_reducers/ctapipe_integrators.py +++ b/CHECLabPy/waveform_reducers/ctapipe_integrators.py @@ -1,5 +1,6 @@ from CHECLabPy.core.reducer import WaveformReducer, column -from CHECLabPy.utils.mapping import get_ctapipe_camera_geometry +from CHECLabPy.utils.mapping import get_ctapipe_subarray +import numpy as np class CtapipeLocalPeakIntegrator(WaveformReducer): @@ -69,20 +70,21 @@ def __init__(self, n_pixels, n_samples, mapping=None, **kwargs): .format(self.__class__.__name__, self.columns)) raise ImportError(msg) - camera = get_ctapipe_camera_geometry(mapping) - + subarray, self._tel_id = get_ctapipe_subarray(mapping) + n_pixels = subarray.tel[self._tel_id].camera.geometry.n_pixels + self.channel = np.zeros(n_pixels, dtype=int) self.window_size = self.kwargs.get("window_size", 6) self.window_shift = self.kwargs.get("window_shift", 3) self.integrator = NeighborPeakWindowSum( window_shift=self.window_shift, window_width=self.window_size, lwt=0, + subarray=subarray, ) - self.integrator.neighbors = camera.neighbor_matrix_where def _prepare(self, waveforms): super()._prepare(waveforms) - charge, pulse_time = self.integrator(waveforms) + charge, pulse_time = self.integrator(waveforms, self._tel_id, self.channel) self.t = pulse_time self.charge = charge