From f9fcc0f145a19801135fa6e7249305840ebbb8e8 Mon Sep 17 00:00:00 2001 From: Marcel Strzys Date: Tue, 13 Jun 2023 19:58:00 +0900 Subject: [PATCH 01/18] feat: added DL3 data reader for LST dl3 files - added an event class for LST dl3 files - tested that class can read LST dl3 files - entire chain not tested - DL3 format for other instruments (e.g. MAGIC) can be slightly different and mya require a different reader ToDo: - MAGIC DL3 reader - test entire chain with LST DL3 files --- src/pybkgmodel/data.py | 214 ++++++++++++++++++++++++++++++++++++++--- 1 file changed, 200 insertions(+), 14 deletions(-) diff --git a/src/pybkgmodel/data.py b/src/pybkgmodel/data.py index ce95324..67685c3 100644 --- a/src/pybkgmodel/data.py +++ b/src/pybkgmodel/data.py @@ -3,6 +3,7 @@ import numpy import pandas import uproot +from astropy.io import fits import astropy.time import astropy.units as u @@ -28,12 +29,14 @@ def find_run_neighbours(target_run, run_list, time_delta, pointing_delta): """ neihbours = filter( - lambda run_: (abs(run_.mjd_start - target_run.mjd_stop)*u.d < time_delta) or (abs(run_.mjd_stop - target_run.mjd_start)*u.d < time_delta), + lambda run_: (abs(run_.mjd_start - target_run.mjd_stop)*u.d < time_delta) or + (abs(run_.mjd_stop - target_run.mjd_start)*u.d < time_delta), run_list ) neihbours = filter( - lambda run_: target_run.tel_pointing_start.icrs.separation(run_.tel_pointing_start.icrs) < pointing_delta, + lambda run_: target_run.tel_pointing_start.icrs.separation(run_.tel_pointing_start.icrs) + < pointing_delta, neihbours ) @@ -41,11 +44,13 @@ def find_run_neighbours(target_run, run_list, time_delta, pointing_delta): class EventSample: + """_summary_ + """ def __init__( self, event_ra, event_dec, event_energy, pointing_ra, pointing_dec, pointing_az, pointing_zd, - mjd, delta_t + mjd, delta_t, eff_obs_time ): self.__event_ra = event_ra self.__event_dec = event_dec @@ -56,7 +61,10 @@ def __init__( self.__pointing_zd = pointing_zd self.__mjd = mjd self.__delta_t = delta_t - self.__eff_obs_time = self.calc_eff_obs_time() + if eff_obs_time is None: + self.__eff_obs_time = self.calc_eff_obs_time() + else: + self.__eff_obs_time = eff_obs_time @property def delta_t(self): @@ -103,6 +111,13 @@ def mjd(self): return self.__mjd def calc_eff_obs_time(self): + """_summary_ + + Returns + ------- + _type_ + _description_ + """ mjd_sorted = numpy.sort(self.__mjd) time_diff = numpy.diff(mjd_sorted) @@ -132,6 +147,8 @@ def calc_eff_obs_time(self): class EventFile: + """_summary_ + """ file_name = '' obs_id = None @@ -197,7 +214,14 @@ def mjd(self): return self.events.mjd -class MagicEventFile(EventFile): +class MagicRootEventFile(EventFile): + """_summary_ + + Parameters + ---------- + EventFile : _type_ + _description_ + """ def __init__(self, file_name, cuts=None): super().__init__(file_name, cuts) @@ -354,13 +378,21 @@ def load_events(cls, file_name, cuts): event_data['pointing_az'], event_data['pointing_zd'], event_data['mjd'], - event_data['delta_t'] + event_data['delta_t'], + None ) return event_sample -class LstEventFile(EventFile): +class LstDl2EventFile(EventFile): + """_summary_ + + Parameters + ---------- + EventFile : _type_ + _description_ + """ def __init__(self, file_name, cuts=None): super().__init__(file_name, cuts) @@ -399,7 +431,8 @@ def load_events(cls, file_name, cuts): Returns ------- dict: - A dictionary with the even properties: charge / arrival time data, trigger, direction etc. + A dictionary with the even properties: charge / arrival time data, + trigger, direction etc. """ data_units = { @@ -475,7 +508,7 @@ def load_events(cls, file_name, cuts): event_data['event_ra'] = coords.ra.to(data_units['event_ra']).value event_data['event_dec'] = coords.dec.to(data_units['event_dec']).value - except: + except KeyError: # The file is likely corrupted, so return empty arrays print("The file is corrupted or is missing the event tree. Empty arrays will be returned.") for key in data_names_mapping: @@ -501,23 +534,176 @@ def load_events(cls, file_name, cuts): event_data['pointing_az'], event_data['pointing_zd'], event_data['mjd'], - event_data['delta_t'] + event_data['delta_t'], + None ) return event_sample +class LstDl3EventFile(EventFile): + """_summary_ + + Parameters + ---------- + EventFile : _type_ + _description_ + """ + def __init__(self, file_name, cuts=None): + super().__init__(file_name, cuts) + + self.file_name = file_name + self.obs_id = self.get_obs_id(file_name) + self.events = self.load_events(file_name, cuts) + + @classmethod + def is_compatible(cls, file_name): + _, ext = os.path.splitext(file_name) + compatible = ext.lower() == ".fits" + return compatible + + @classmethod + def get_obs_id(cls, file_name): + parsed = re.findall('.*dl3_LST-1.Run(\d+).h5', file_name) + if parsed: + obs_id = int(parsed[0]) + else: + raise RuntimeError(f'Can not find observations ID in {file_name}') + + return obs_id + + @classmethod + def load_events(cls, file_name, cuts): + """_summary_ + + Parameters + ---------- + file_name : _type_ + _description_ + cuts : _type_ + _description_ + + Returns + ------- + _type_ + _description_ + """ + + data_units = { + 'event_ra': u.hourangle, + 'event_dec': u.deg, + 'event_energy': u.GeV, + 'pointing_ra': u.deg, + 'pointing_dec': u.deg, + 'pointing_az': u.deg, + 'pointing_zd': u.deg, + 'mjd': u.d, + 'delta_t': u.s, + 'gammaness':u.one + } + + data_names_mapping = { + 'EVENT_ID': 'daq_event_number', + 'RA': 'event_ra', + 'DEC': 'event_dec', + 'GAMMANESS': 'gammaness', + 'ENERGY': 'event_energy', + 'TIME':'mjd', + 'AZ_PNT': 'pointing_az', + 'ZD_PNT': 'pointing_zd', + 'RA_PNT':'pointing_ra', + 'DEC_PNT':'pointing_dec' + } + + with fits.open(file_name, memmap=False) as input_file: + try: + evt_head = input_file["EVENTS"].header + evt_data = pandas.DataFrame(input_file["EVENTS"].data) + + event_data = {} + + for key in data_names_mapping: + name = data_names_mapping[key] + + if key in evt_data.keys(): + event_data[name] = evt_data[key].to_numpy() + + # Event times need to be converted from LST Epoch + LST_EPOCH = astropy.time.Time('2018-10-01T00:00:00', scale='utc') + event_data['mjd'] = astropy.time.Time(event_data['mjd'], format='unix') + event_data['mjd'] = astropy.time.Time((event_data['mjd'].unix + LST_EPOCH.unix), + scale='utc', + format='unix' + ).mjd * data_units['mjd'] + + # Compute the telescope pointing positions for each event + lst_time = astropy.time.Time(event_data['mjd'], format='mjd') + lst_loc = EarthLocation(lat=28.761758*u.deg, + lon=-17.890659*u.deg, + height=2200*u.m) + alt_az_frame = AltAz(obstime=lst_time, location=lst_loc) + coords = SkyCoord(evt_head['RA_PNT'] *u.deg, + evt_head['DEC_PNT'] *u.deg, + frame='icrs') + + altaz_pointing = coords.transform_to(alt_az_frame) + + event_data['pointing_zd'] = 90 - altaz_pointing.alt.to( + data_units['pointing_zd'] + ).value + event_data['pointing_az'] = altaz_pointing.az.to(data_units['pointing_az']).value + event_data['pointing_ra'] = [evt_head['RA_PNT']] * len(event_data['pointing_zd']) + event_data['pointing_ra'] = numpy.array(event_data['pointing_ra']) + event_data['pointing_dec'] = [evt_head['DEC_PNT']] * len(event_data['pointing_zd']) + event_data['pointing_dec'] = numpy.array(event_data['pointing_dec']) + + except KeyError: + print(f"File {file_name} corrupted or missing the Events hdu." + + "Empty arrays will be returned.") + + finite = [numpy.isfinite(event_data[key]) for key in event_data if event_data[key] is not None] + all_finite = numpy.prod(finite, axis=0, dtype=bool) + + for key, item in event_data.items(): + if item is not None: + event_data[key] = item[all_finite] + + if key in data_units: + event_data[key] = item * data_units[key] + + + event_sample = EventSample( + event_data['event_ra'], + event_data['event_dec'], + event_data['event_energy'], + event_data['pointing_ra'], + event_data['pointing_dec'], + event_data['pointing_az'], + event_data['pointing_zd'], + event_data['mjd'], + None, + numpy.array(evt_head['LIVETIME']) * u.s, + ) + + return event_sample class RunSummary: + """_summary_ + + Raises + ------ + RuntimeError + _description_ + """ __obs_id = None __file_name = None __tel_pointing_start = None __tel_pointing_stop = None def __init__(self, file_name): - if MagicEventFile.is_compatible(file_name): - events = MagicEventFile(file_name) - elif LstEventFile.is_compatible(file_name): - events = LstEventFile(file_name) + if MagicRootEventFile.is_compatible(file_name): + events = MagicRootEventFile(file_name) + elif LstDl2EventFile.is_compatible(file_name): + events = LstDl2EventFile(file_name) else: raise RuntimeError(f"Unsupported file format for '{file_name}'.") From 475cdba33630bef5f4c294b56978a64ef7697f0e Mon Sep 17 00:00:00 2001 From: Marcel Strzys Date: Tue, 13 Jun 2023 20:22:41 +0900 Subject: [PATCH 02/18] - just made the LstDl3event class usable by the higher classes --- src/pybkgmodel/model.py | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/src/pybkgmodel/model.py b/src/pybkgmodel/model.py index b6db989..78949f0 100644 --- a/src/pybkgmodel/model.py +++ b/src/pybkgmodel/model.py @@ -3,7 +3,10 @@ from astropy.coordinates import SkyCoord -from pybkgmodel.data import MagicEventFile, LstEventFile +from pybkgmodel.data import (MagicRootEventFile, + LstDl2EventFile, + LstDl3EventFile + ) from pybkgmodel.data import find_run_neighbours from pybkgmodel.camera import RectangularCameraImage @@ -41,17 +44,23 @@ def read_runs(self, target_run, neighbours, cuts): ------ RuntimeError Raise if a run in an unsupported format is provided. - Currenttly supported formats are DL2 for LST and ROOT for MAGIC. + Currenttly supported formats are DL2 and DL3 for LST and ROOT for MAGIC. """ - if MagicEventFile.is_compatible(target_run.file_name): + if MagicRootEventFile.is_compatible(target_run.file_name): evtfiles = [ - MagicEventFile(run.file_name, cuts=cuts) + MagicRootEventFile(run.file_name, cuts=cuts) for run in (target_run,) + neighbours ] return evtfiles - elif LstEventFile.is_compatible(target_run.file_name): + elif LstDl2EventFile.is_compatible(target_run.file_name): evtfiles = [ - LstEventFile(run.file_name, cuts=cuts) + LstDl2EventFile(run.file_name, cuts=cuts) + for run in (target_run,) + neighbours + ] + return evtfiles + elif LstDl3EventFile.is_compatible(target_run.file_name): + evtfiles = [ + LstDl3EventFile(run.file_name) for run in (target_run,) + neighbours ] return evtfiles From 61925ca50d6ca53be2df86ee57e3335331475551 Mon Sep 17 00:00:00 2001 From: Marcel Strzys Date: Sun, 18 Jun 2023 01:16:46 +0900 Subject: [PATCH 03/18] bug fix: tested the DL3 reader with LST-1 dl3 files Found and corrected the following issues, some concern also previous commits to the main branch: - issue of the unti of the event time variable - wrong unit of the ra event variable - dl2 file reader had an error if the ra variable of the pointing position was not precomputed - wrong overwriting of the setter function for the bkgmap maker in the bkg classes --- src/pybkgmodel/camera.py | 12 ++++++------ src/pybkgmodel/data.py | 14 ++++++++------ src/pybkgmodel/processing.py | 8 ++++---- 3 files changed, 18 insertions(+), 16 deletions(-) diff --git a/src/pybkgmodel/camera.py b/src/pybkgmodel/camera.py index 62401df..3a4166b 100644 --- a/src/pybkgmodel/camera.py +++ b/src/pybkgmodel/camera.py @@ -169,7 +169,7 @@ def __init__(self, counts, xedges, yedges, energy_edges, center=None, mask=None, ny = yedges.size - 1 if mask is None: - mask = numpy.ones((nx, ny), dtype=numpy.bool) + mask = numpy.ones((nx, ny), dtype=bool) if exposure is None: exposure = numpy.ones((nx, ny), dtype=numpy.float) * u.s @@ -244,7 +244,7 @@ def differential_rate(self, index=None): Parameters ---------- index: float - Power law spectral index to assume. + Power law spectral index to assume. If none, will be dynamically determined assuming a "node function" for the spectral shape. @@ -322,14 +322,14 @@ def mask_region(self, region): dummy_wcs.wcs.crval = [0, -90] dummy_wcs.wcs.ctype = ["RA---AIR", "DEC--AIR"] dummy_wcs.wcs.set_pv([(2, 1, 45.0)]) - + in_region = region.contains(self.pixel_coords, dummy_wcs) self.mask[in_region] = False - + def mask_reset(self): """_summary_ """ - self.mask = numpy.ones((self.xedges.size - 1, self.yedges.size - 1), dtype=numpy.bool) + self.mask = numpy.ones((self.xedges.size - 1, self.yedges.size - 1), dtype=bool) def plot(self, energy_bin_id=0, ax_unit='deg', val_unit='1/s', **kwargs): @@ -394,7 +394,7 @@ def get_pixel_areas(self): for i in range(nx): for j in range(ny): area[i, j] = solid_angle_lat_lon_rectangle(self.xedges[i], self.xedges[i+1], self.yedges[j], self.yedges[j+1]) - + return area def to_hdu(self, name='BACKGROUND'): diff --git a/src/pybkgmodel/data.py b/src/pybkgmodel/data.py index 67685c3..476f839 100644 --- a/src/pybkgmodel/data.py +++ b/src/pybkgmodel/data.py @@ -496,13 +496,13 @@ def load_events(cls, file_name, cuts): lst_loc = EarthLocation(lat=28.761758*u.deg, lon=-17.890659*u.deg, height=2200*u.m) alt_az_frame = AltAz(obstime=lst_time, location=lst_loc) - if 'pointing_ra' not in event_data: + if event_data['pointing_ra'] is None: coords = SkyCoord(alt=data['alt_tel'].to_numpy()*u.rad, az=data['az_tel'].to_numpy()*u.rad, frame=alt_az_frame).icrs event_data['pointing_ra'] = coords.ra.to(data_units['pointing_ra']).value event_data['pointing_dec'] = coords.dec.to(data_units['pointing_dec']).value - if 'event_ra' not in event_data: + if event_data['event_ra'] is None: coords = SkyCoord(alt=data['reco_alt'].to_numpy()*u.rad, az=data['reco_az'].to_numpy()*u.rad, frame=alt_az_frame).icrs event_data['event_ra'] = coords.ra.to(data_units['event_ra']).value @@ -563,7 +563,7 @@ def is_compatible(cls, file_name): @classmethod def get_obs_id(cls, file_name): - parsed = re.findall('.*dl3_LST-1.Run(\d+).h5', file_name) + parsed = re.findall('.*dl3_LST-1.Run(\d+).fits', file_name) if parsed: obs_id = int(parsed[0]) else: @@ -589,9 +589,9 @@ def load_events(cls, file_name, cuts): """ data_units = { - 'event_ra': u.hourangle, + 'event_ra': u.deg, 'event_dec': u.deg, - 'event_energy': u.GeV, + 'event_energy': u.TeV, 'pointing_ra': u.deg, 'pointing_dec': u.deg, 'pointing_az': u.deg, @@ -633,7 +633,7 @@ def load_events(cls, file_name, cuts): event_data['mjd'] = astropy.time.Time((event_data['mjd'].unix + LST_EPOCH.unix), scale='utc', format='unix' - ).mjd * data_units['mjd'] + ).mjd # Compute the telescope pointing positions for each event lst_time = astropy.time.Time(event_data['mjd'], format='mjd') @@ -704,6 +704,8 @@ def __init__(self, file_name): events = MagicRootEventFile(file_name) elif LstDl2EventFile.is_compatible(file_name): events = LstDl2EventFile(file_name) + elif LstDl3EventFile.is_compatible(file_name): + events = LstDl3EventFile(file_name) else: raise RuntimeError(f"Unsupported file format for '{file_name}'.") diff --git a/src/pybkgmodel/processing.py b/src/pybkgmodel/processing.py index 4e9514d..341f365 100644 --- a/src/pybkgmodel/processing.py +++ b/src/pybkgmodel/processing.py @@ -519,7 +519,7 @@ def __init__(self, def bkg_map_maker(self, maker): if not isinstance(maker, WobbleMap): raise TypeError(f"Maker must be of type {WobbleMap}") - BkgMakerBase.bkg_map_maker = maker + BkgMakerBase.bkg_map_maker.fset(self, maker) class StackedWobbleMap(Stacked): """ @@ -645,7 +645,7 @@ def __init__(self, def bkg_map_maker(self, maker): if not isinstance(maker, WobbleMap): raise TypeError(f"Maker must be of type {WobbleMap}") - BkgMakerBase.bkg_map_maker = maker + BkgMakerBase.bkg_map_maker.fset(self, maker) class RunwiseExclusionMap(Runwise): """ @@ -781,7 +781,7 @@ def __init__(self, def bkg_map_maker(self, maker): if not isinstance(maker, ExclusionMap): raise TypeError(f"Maker must be of type {ExclusionMap}") - BkgMakerBase.bkg_map_maker = maker + BkgMakerBase.bkg_map_maker.fset(self, maker) class StackedExclusionMap(Stacked): """ @@ -917,4 +917,4 @@ def __init__(self, def bkg_map_maker(self, maker): if not isinstance(maker, ExclusionMap): raise TypeError(f"Maker must be of type {ExclusionMap}") - BkgMakerBase.bkg_map_maker = maker + BkgMakerBase.bkg_map_maker.fset(self, maker) From 6a9153366cf080d0174e8bead95f180d9a2205a7 Mon Sep 17 00:00:00 2001 From: Marcel Strzys Date: Wed, 21 Jun 2023 00:17:42 +0900 Subject: [PATCH 04/18] bug fix: wrong attribute in the stacked maps class: get_maps() attempted to set attribute bkg_maps instead of the correct _bkg_maps --- src/pybkgmodel/processing.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pybkgmodel/processing.py b/src/pybkgmodel/processing.py index 341f365..e75169b 100644 --- a/src/pybkgmodel/processing.py +++ b/src/pybkgmodel/processing.py @@ -390,7 +390,7 @@ def get_maps(self): self.out_dir, f"{self.out_prefix}stacked_bkg_map.fits" ) - self.bkg_maps = {stacked_name: stacked_map} + self._bkg_maps = {stacked_name: stacked_map} self.write_maps(bkg_maps=self.bkg_maps, overwrite=self.overwrite) return self.bkg_maps From d0fa65a728060c57b199e4a7095f11f1a52273de Mon Sep 17 00:00:00 2001 From: Marcel Strzys Date: Sun, 25 Jun 2023 15:18:38 +0900 Subject: [PATCH 05/18] - just remove unnecessary comma --- src/pybkgmodel/data.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pybkgmodel/data.py b/src/pybkgmodel/data.py index 476f839..e611df4 100644 --- a/src/pybkgmodel/data.py +++ b/src/pybkgmodel/data.py @@ -681,7 +681,7 @@ def load_events(cls, file_name, cuts): event_data['pointing_zd'], event_data['mjd'], None, - numpy.array(evt_head['LIVETIME']) * u.s, + numpy.array(evt_head['LIVETIME']) * u.s ) return event_sample From 4c79f9702ad016f939b4f6ae22e39e2393c3e497 Mon Sep 17 00:00:00 2001 From: Marcel Strzys Date: Tue, 19 Mar 2024 18:42:35 +0900 Subject: [PATCH 06/18] : implement the dl3 reader - Modified the reader so it can read DL3 file regardless the instrument. - Data units are read from hdu header as long as provided. - Tested with LST and MAGIC DL3 files. TODO: - Observatory location not part of the mandatory header keywords. May need to be provided in input file by the user. Will be part of future pull request - Drift method implemented but not tested due to lack of corresponding data. --- src/pybkgmodel/data.py | 148 +++++++++++++++++++++++------------------ 1 file changed, 85 insertions(+), 63 deletions(-) diff --git a/src/pybkgmodel/data.py b/src/pybkgmodel/data.py index e611df4..5d8f630 100644 --- a/src/pybkgmodel/data.py +++ b/src/pybkgmodel/data.py @@ -540,7 +540,7 @@ def load_events(cls, file_name, cuts): return event_sample -class LstDl3EventFile(EventFile): +class DL3EventFile(EventFile): """_summary_ Parameters @@ -548,12 +548,12 @@ class LstDl3EventFile(EventFile): EventFile : _type_ _description_ """ - def __init__(self, file_name, cuts=None): - super().__init__(file_name, cuts) + def __init__(self, file_name): + super().__init__(file_name) self.file_name = file_name self.obs_id = self.get_obs_id(file_name) - self.events = self.load_events(file_name, cuts) + self.events = self.load_events(file_name) @classmethod def is_compatible(cls, file_name): @@ -563,7 +563,7 @@ def is_compatible(cls, file_name): @classmethod def get_obs_id(cls, file_name): - parsed = re.findall('.*dl3_LST-1.Run(\d+).fits', file_name) + parsed = re.findall('.*dl3_LST-1.Run(\d+).*', file_name) if parsed: obs_id = int(parsed[0]) else: @@ -572,105 +572,127 @@ def get_obs_id(cls, file_name): return obs_id @classmethod - def load_events(cls, file_name, cuts): + def load_events(cls, file_name): """_summary_ Parameters ---------- - file_name : _type_ - _description_ - cuts : _type_ - _description_ + file_name: str + Name of the DL3 file to use. Returns ------- - _type_ - _description_ + dict: + A dictionary with the event properties: arrival time, direction, and energy. """ - data_units = { - 'event_ra': u.deg, - 'event_dec': u.deg, - 'event_energy': u.TeV, - 'pointing_ra': u.deg, - 'pointing_dec': u.deg, - 'pointing_az': u.deg, - 'pointing_zd': u.deg, - 'mjd': u.d, - 'delta_t': u.s, - 'gammaness':u.one - } - data_names_mapping = { 'EVENT_ID': 'daq_event_number', 'RA': 'event_ra', 'DEC': 'event_dec', 'GAMMANESS': 'gammaness', 'ENERGY': 'event_energy', - 'TIME':'mjd', - 'AZ_PNT': 'pointing_az', - 'ZD_PNT': 'pointing_zd', - 'RA_PNT':'pointing_ra', - 'DEC_PNT':'pointing_dec' + 'TIME':'mjd' } with fits.open(file_name, memmap=False) as input_file: try: - evt_head = input_file["EVENTS"].header - evt_data = pandas.DataFrame(input_file["EVENTS"].data) + evt_hdu = input_file["EVENTS"] + evt_head = evt_hdu.header + evt_data = pandas.DataFrame(evt_hdu.data) event_data = {} - for key in data_names_mapping: - name = data_names_mapping[key] + for key, name in data_names_mapping.items(): if key in evt_data.keys(): event_data[name] = evt_data[key].to_numpy() + unit_name = f"TUNIT{evt_data.columns.get_loc(key)+1}" + + try: + event_data[name] *= u.Unit(evt_head[unit_name]) + except KeyError: + event_data[name] *= u.one + + # Event times need to be converted from Instrument reference epoch + ref_epoch = astropy.time.Time(evt_head['MJDREFI']+evt_head['MJDREFF'], format='mjd') - # Event times need to be converted from LST Epoch - LST_EPOCH = astropy.time.Time('2018-10-01T00:00:00', scale='utc') event_data['mjd'] = astropy.time.Time(event_data['mjd'], format='unix') - event_data['mjd'] = astropy.time.Time((event_data['mjd'].unix + LST_EPOCH.unix), + event_data['mjd'] = astropy.time.Time((event_data['mjd'].unix + ref_epoch.unix), scale='utc', format='unix' - ).mjd + ).mjd * u.d - # Compute the telescope pointing positions for each event - lst_time = astropy.time.Time(event_data['mjd'], format='mjd') - lst_loc = EarthLocation(lat=28.761758*u.deg, - lon=-17.890659*u.deg, - height=2200*u.m) - alt_az_frame = AltAz(obstime=lst_time, location=lst_loc) - coords = SkyCoord(evt_head['RA_PNT'] *u.deg, - evt_head['DEC_PNT'] *u.deg, - frame='icrs') - altaz_pointing = coords.transform_to(alt_az_frame) + evt_time = astropy.time.Time(event_data['mjd'], format='mjd') + + # TODO: current observatory location only La Palma, no mandatory header keyword + obs_loc = EarthLocation(lat=28.761758*u.deg, + lon=-17.890659*u.deg, + height=2200*u.m) + + if evt_head['OBS_MODE'] == 'POINTING': + + alt_az_frame = AltAz(obstime=evt_time, + location=obs_loc) + + coords = SkyCoord(evt_head['RA_PNT'] *u.deg, + evt_head['DEC_PNT'] *u.deg, + frame='icrs') + + altaz_pointing = coords.transform_to(alt_az_frame) + + event_data['pointing_zd'] = 90 * u.deg - altaz_pointing.alt + event_data['pointing_az'] = altaz_pointing.az.to(u.deg) + - event_data['pointing_zd'] = 90 - altaz_pointing.alt.to( - data_units['pointing_zd'] - ).value - event_data['pointing_az'] = altaz_pointing.az.to(data_units['pointing_az']).value - event_data['pointing_ra'] = [evt_head['RA_PNT']] * len(event_data['pointing_zd']) - event_data['pointing_ra'] = numpy.array(event_data['pointing_ra']) - event_data['pointing_dec'] = [evt_head['DEC_PNT']] * len(event_data['pointing_zd']) - event_data['pointing_dec'] = numpy.array(event_data['pointing_dec']) + event_data['pointing_ra'] = numpy.array( + [evt_head['RA_PNT']] \ + * len(event_data['pointing_zd']) + ) * u.deg + event_data['pointing_dec'] = numpy.array( + [evt_head['DEC_PNT']] \ + * len(event_data['pointing_zd']) + ) * u.deg + + elif evt_head['OBS_MODE'] == 'DRIFT': + # TODO: function not tested yet, since no data at hand to test, hence + # preliminary implementation + + coords = SkyCoord(evt_head['ALT_PNT'] *u.deg, + evt_head['AZ_PNT'] *u.deg, + obstime=astropy.time.Time(event_data['mjd'], format='mjd'), + location=obs_loc, + frame='altaz') + + radec_pointing = coords.transform_to('icrs') + + event_data['pointing_zd'] = 90 * u.deg - coords.alt + event_data['pointing_az'] = coords.az + event_data['pointing_ra'] = radec_pointing.ra + event_data['pointing_dec'] = radec_pointing.dec + + else: + print("Observation mode currently not supported.\ + Function will return zeros for the event location.\ + Supported modes: POINTING, DRIFT") + + event_data['pointing_zd'] = numpy.zeros(len(event_data['mjd'])) * u.deg + event_data['pointing_az'] = numpy.zeros(len(event_data['mjd'])) * u.deg + event_data['pointing_ra'] = numpy.zeros(len(event_data['mjd'])) * u.deg + event_data['pointing_dec'] = numpy.zeros(len(event_data['mjd'])) * u.deg except KeyError: print(f"File {file_name} corrupted or missing the Events hdu." + "Empty arrays will be returned.") - finite = [numpy.isfinite(event_data[key]) for key in event_data if event_data[key] is not None] + finite = [numpy.isfinite(item) for key, item in event_data.items() if item is not None] all_finite = numpy.prod(finite, axis=0, dtype=bool) for key, item in event_data.items(): if item is not None: event_data[key] = item[all_finite] - if key in data_units: - event_data[key] = item * data_units[key] - - event_sample = EventSample( event_data['event_ra'], event_data['event_dec'], @@ -704,8 +726,8 @@ def __init__(self, file_name): events = MagicRootEventFile(file_name) elif LstDl2EventFile.is_compatible(file_name): events = LstDl2EventFile(file_name) - elif LstDl3EventFile.is_compatible(file_name): - events = LstDl3EventFile(file_name) + elif DL3EventFile.is_compatible(file_name): + events = DL3EventFile(file_name) else: raise RuntimeError(f"Unsupported file format for '{file_name}'.") From bdf3989a9113fb8397bc621ca91ae97dfff1645b Mon Sep 17 00:00:00 2001 From: Marcel Strzys Date: Thu, 21 Mar 2024 18:08:59 +0900 Subject: [PATCH 07/18] Minor modifications to the reader for improved compatibility and performance - added some functionalities suggested by Lukas Nickel, but adapted them to the current version - changed the function for retrieving the Obs ID from filename to reading it from the hdu header - changed the is_compatible, so it also identifies compressed files - added the erfaAstromInterpolator for the coordinate transform for speed up - added docstrings to all functions related to the DL3 reader --- src/pybkgmodel/data.py | 63 ++++++++++++++++++++++++++++++++++-------- 1 file changed, 51 insertions(+), 12 deletions(-) diff --git a/src/pybkgmodel/data.py b/src/pybkgmodel/data.py index 5d8f630..b71d28b 100644 --- a/src/pybkgmodel/data.py +++ b/src/pybkgmodel/data.py @@ -1,5 +1,6 @@ import os import re +from pathlib import Path import numpy import pandas import uproot @@ -8,6 +9,7 @@ import astropy.units as u from astropy.coordinates import SkyCoord, EarthLocation, AltAz +from astropy.coordinates.erfa_astrom import erfa_astrom, ErfaAstromInterpolator def find_run_neighbours(target_run, run_list, time_delta, pointing_delta): @@ -541,12 +543,13 @@ def load_events(cls, file_name, cuts): return event_sample class DL3EventFile(EventFile): - """_summary_ + """Reader for DL3 data file compliant with the GADF. For details see + https://gamma-astro-data-formats.readthedocs.io/ Parameters ---------- - EventFile : _type_ - _description_ + file_name: str + Name of the DL3 file to use. """ def __init__(self, file_name): super().__init__(file_name) @@ -557,23 +560,57 @@ def __init__(self, file_name): @classmethod def is_compatible(cls, file_name): - _, ext = os.path.splitext(file_name) - compatible = ext.lower() == ".fits" + """Function checks whether the file is a fits file according to the file extension. + + Parameters + ---------- + file_name: str + Name of the DL3 file to use. + + Returns + ------- + bool + True if fits file according to file extension + """ + + ext = Path(file_name).suffixes + compatible = ".fits" in ext + return compatible @classmethod def get_obs_id(cls, file_name): - parsed = re.findall('.*dl3_LST-1.Run(\d+).*', file_name) - if parsed: - obs_id = int(parsed[0]) - else: - raise RuntimeError(f'Can not find observations ID in {file_name}') + """Function reads the observation id from the fits file. + + Parameters + ---------- + file_name: str + Name of the DL3 file to use. + + Returns + ------- + int + Observation ID number of the fits file + + Raises + ------ + RuntimeError + In case Observation ID is not found or Events HDU layer is missing. In this case the file is not complient to GADF. + """ + + obs_id = None + + with fits.open(file_name, memmap=False) as input_file: + try: + obs_id = int(input_file["EVENTS"].header['OBS_ID']) + except Exception as error: + raise RuntimeError(f'Can not find observations ID in {file_name}') from error return obs_id @classmethod def load_events(cls, file_name): - """_summary_ + """Function to read the events from the DL3 file and to compute the missing variables needed by pybkgmodel Parameters ---------- @@ -595,7 +632,8 @@ def load_events(cls, file_name): 'TIME':'mjd' } - with fits.open(file_name, memmap=False) as input_file: + with fits.open(file_name, memmap=False) as input_file, \ + erfa_astrom.set(ErfaAstromInterpolator(1 * u.s)): try: evt_hdu = input_file["EVENTS"] evt_head = evt_hdu.header @@ -607,6 +645,7 @@ def load_events(cls, file_name): if key in evt_data.keys(): event_data[name] = evt_data[key].to_numpy() + # Fits header keywords start counting from 1 unit_name = f"TUNIT{evt_data.columns.get_loc(key)+1}" try: From 0358b265b0cce581b4a3fc253345536c534dea6f Mon Sep 17 00:00:00 2001 From: Marcel Strzys Date: Thu, 21 Mar 2024 18:20:56 +0900 Subject: [PATCH 08/18] adapted the model.py file to use the new DL3 reader. --- src/pybkgmodel/model.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/pybkgmodel/model.py b/src/pybkgmodel/model.py index 78949f0..e3ed8e0 100644 --- a/src/pybkgmodel/model.py +++ b/src/pybkgmodel/model.py @@ -5,7 +5,7 @@ from pybkgmodel.data import (MagicRootEventFile, LstDl2EventFile, - LstDl3EventFile + DL3EventFile ) from pybkgmodel.data import find_run_neighbours @@ -44,7 +44,7 @@ def read_runs(self, target_run, neighbours, cuts): ------ RuntimeError Raise if a run in an unsupported format is provided. - Currenttly supported formats are DL2 and DL3 for LST and ROOT for MAGIC. + Currenttly supported formats are DL3 according to GADF, DL2 for LST and ROOT for MAGIC. """ if MagicRootEventFile.is_compatible(target_run.file_name): evtfiles = [ @@ -58,9 +58,9 @@ def read_runs(self, target_run, neighbours, cuts): for run in (target_run,) + neighbours ] return evtfiles - elif LstDl3EventFile.is_compatible(target_run.file_name): + elif DL3EventFile.is_compatible(target_run.file_name): evtfiles = [ - LstDl3EventFile(run.file_name) + DL3EventFile(run.file_name) for run in (target_run,) + neighbours ] return evtfiles From 61e2e9adeb5b31b007545b4963bbd60440b818d2 Mon Sep 17 00:00:00 2001 From: Marcel Strzys Date: Thu, 21 Mar 2024 18:51:30 +0900 Subject: [PATCH 09/18] Added the Obs mode WOBBLE in addition to POINTING. Seems the Obs modes are not to well defined e.g. the H.E.S.S. public data use WOBBLE as the obs mode. Hence, added it to support such files until the allowed modes are better defined. --- src/pybkgmodel/data.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pybkgmodel/data.py b/src/pybkgmodel/data.py index b71d28b..819efee 100644 --- a/src/pybkgmodel/data.py +++ b/src/pybkgmodel/data.py @@ -670,7 +670,7 @@ def load_events(cls, file_name): lon=-17.890659*u.deg, height=2200*u.m) - if evt_head['OBS_MODE'] == 'POINTING': + if evt_head['OBS_MODE'] in ('POINTING', 'WOBBLE'): alt_az_frame = AltAz(obstime=evt_time, location=obs_loc) From 4a46b01583ddabacc74ace9832d676b05df1be34 Mon Sep 17 00:00:00 2001 From: Marcel Strzys Date: Fri, 7 Jun 2024 18:36:54 +0900 Subject: [PATCH 10/18] refactor: change variables and data handling to accommodate comments from the corresponding PR (#7). - change name from Dl2 to DL2 to agree with the DL3 naming - changed to way the event times are read from DL3 files to be consistent with the 'mjd' naming of the dictionary - the functionality of the code was not changed --- src/pybkgmodel/data.py | 14 +++++++------- src/pybkgmodel/model.py | 6 +++--- src/pybkgmodel/processing.py | 7 ++++--- 3 files changed, 14 insertions(+), 13 deletions(-) diff --git a/src/pybkgmodel/data.py b/src/pybkgmodel/data.py index 819efee..5566a4b 100644 --- a/src/pybkgmodel/data.py +++ b/src/pybkgmodel/data.py @@ -387,7 +387,8 @@ def load_events(cls, file_name, cuts): return event_sample -class LstDl2EventFile(EventFile): +#Adapt +class LstDL2EventFile(EventFile): """_summary_ Parameters @@ -629,7 +630,6 @@ def load_events(cls, file_name): 'DEC': 'event_dec', 'GAMMANESS': 'gammaness', 'ENERGY': 'event_energy', - 'TIME':'mjd' } with fits.open(file_name, memmap=False) as input_file, \ @@ -656,13 +656,13 @@ def load_events(cls, file_name): # Event times need to be converted from Instrument reference epoch ref_epoch = astropy.time.Time(evt_head['MJDREFI']+evt_head['MJDREFF'], format='mjd') - event_data['mjd'] = astropy.time.Time(event_data['mjd'], format='unix') - event_data['mjd'] = astropy.time.Time((event_data['mjd'].unix + ref_epoch.unix), + event_data['mjd'] = astropy.time.Time((evt_data['TIME'].to_numpy() + + ref_epoch.unix), scale='utc', format='unix' ).mjd * u.d - + # Adapt evt_time = astropy.time.Time(event_data['mjd'], format='mjd') # TODO: current observatory location only La Palma, no mandatory header keyword @@ -763,8 +763,8 @@ class RunSummary: def __init__(self, file_name): if MagicRootEventFile.is_compatible(file_name): events = MagicRootEventFile(file_name) - elif LstDl2EventFile.is_compatible(file_name): - events = LstDl2EventFile(file_name) + elif LstDL2EventFile.is_compatible(file_name): + events = LstDL2EventFile(file_name) elif DL3EventFile.is_compatible(file_name): events = DL3EventFile(file_name) else: diff --git a/src/pybkgmodel/model.py b/src/pybkgmodel/model.py index e3ed8e0..03e7e44 100644 --- a/src/pybkgmodel/model.py +++ b/src/pybkgmodel/model.py @@ -4,7 +4,7 @@ from astropy.coordinates import SkyCoord from pybkgmodel.data import (MagicRootEventFile, - LstDl2EventFile, + LstDL2EventFile, DL3EventFile ) from pybkgmodel.data import find_run_neighbours @@ -52,9 +52,9 @@ def read_runs(self, target_run, neighbours, cuts): for run in (target_run,) + neighbours ] return evtfiles - elif LstDl2EventFile.is_compatible(target_run.file_name): + elif LstDL2EventFile.is_compatible(target_run.file_name): evtfiles = [ - LstDl2EventFile(run.file_name, cuts=cuts) + LstDL2EventFile(run.file_name, cuts=cuts) for run in (target_run,) + neighbours ] return evtfiles diff --git a/src/pybkgmodel/processing.py b/src/pybkgmodel/processing.py index e75169b..7f925f4 100644 --- a/src/pybkgmodel/processing.py +++ b/src/pybkgmodel/processing.py @@ -55,8 +55,9 @@ class BkgMakerBase: """ - A class used to store the settings from the configuation file and to - facilitate the generation of background maps. + Base class for all processing classes, which store the settings from the + configuation file and facilitate the generation of background maps. + Not intended for direct usage. Attributes ---------- @@ -355,7 +356,7 @@ def get_maps(self): Returns ------- - bkg_maps : dict + bkg_maps : dict Dictionary containing the bkg maps and output names for each run. """ From 57913dea4041e58d05ab22609cefa95e7312162f Mon Sep 17 00:00:00 2001 From: Marcel Strzys Date: Wed, 31 Jul 2024 17:55:32 +0900 Subject: [PATCH 11/18] FIX: - changed the time handling following Max' recommendations in the discussion of the PR - changed the way the fits file compatibility is checked following Michele's recommendations in the discussion of the PR - changed the setter for the bkg_map_makers following Ievgen's recommendation TODO: - Test with data --- src/pybkgmodel/data.py | 25 ++++++++++++++----------- src/pybkgmodel/processing.py | 32 ++++++++++++++++++++++++-------- 2 files changed, 38 insertions(+), 19 deletions(-) diff --git a/src/pybkgmodel/data.py b/src/pybkgmodel/data.py index 5566a4b..ed979dd 100644 --- a/src/pybkgmodel/data.py +++ b/src/pybkgmodel/data.py @@ -574,9 +574,17 @@ def is_compatible(cls, file_name): True if fits file according to file extension """ - ext = Path(file_name).suffixes - compatible = ".fits" in ext + ext = Path(file_name) + + print("PATH", ext) + try: + with fits.open(ext) as file: + pass + compatible = True + except OSError: + compatible = False + return compatible @classmethod @@ -654,17 +662,12 @@ def load_events(cls, file_name): event_data[name] *= u.one # Event times need to be converted from Instrument reference epoch - ref_epoch = astropy.time.Time(evt_head['MJDREFI']+evt_head['MJDREFF'], format='mjd') + ref_epoch = astropy.time.Time(evt_head['MJDREFI'], evt_head['MJDREFF'], format='mjd') - event_data['mjd'] = astropy.time.Time((evt_data['TIME'].to_numpy() - + ref_epoch.unix), - scale='utc', - format='unix' + event_data['mjd'] = astropy.time.Time((evt_data['TIME'].quantity + + ref_epoch) ).mjd * u.d - # Adapt - evt_time = astropy.time.Time(event_data['mjd'], format='mjd') - # TODO: current observatory location only La Palma, no mandatory header keyword obs_loc = EarthLocation(lat=28.761758*u.deg, lon=-17.890659*u.deg, @@ -672,7 +675,7 @@ def load_events(cls, file_name): if evt_head['OBS_MODE'] in ('POINTING', 'WOBBLE'): - alt_az_frame = AltAz(obstime=evt_time, + alt_az_frame = AltAz(obstime=event_data['mjd'], location=obs_loc) coords = SkyCoord(evt_head['RA_PNT'] *u.deg, diff --git a/src/pybkgmodel/processing.py b/src/pybkgmodel/processing.py index 7f925f4..66a5f4a 100644 --- a/src/pybkgmodel/processing.py +++ b/src/pybkgmodel/processing.py @@ -516,11 +516,15 @@ def __init__(self, pointing_delta=self.pointing_delta ) - @BkgMakerBase.bkg_map_maker.setter + @property + def bkg_map_maker(self): + return self.__bkg_map_maker + + @bkg_map_maker.setter def bkg_map_maker(self, maker): if not isinstance(maker, WobbleMap): raise TypeError(f"Maker must be of type {WobbleMap}") - BkgMakerBase.bkg_map_maker.fset(self, maker) + self.__bkg_map_maker = maker class StackedWobbleMap(Stacked): """ @@ -642,11 +646,15 @@ def __init__(self, pointing_delta=self.pointing_delta ) - @BkgMakerBase.bkg_map_maker.setter + @property + def bkg_map_maker(self): + return self.__bkg_map_maker + + @bkg_map_maker.setter def bkg_map_maker(self, maker): if not isinstance(maker, WobbleMap): raise TypeError(f"Maker must be of type {WobbleMap}") - BkgMakerBase.bkg_map_maker.fset(self, maker) + self.__bkg_map_maker = maker class RunwiseExclusionMap(Runwise): """ @@ -778,11 +786,15 @@ def __init__(self, pointing_delta=self.pointing_delta ) - @BkgMakerBase.bkg_map_maker.setter + @property + def bkg_map_maker(self): + return self.__bkg_map_maker + + @bkg_map_maker.setter def bkg_map_maker(self, maker): if not isinstance(maker, ExclusionMap): raise TypeError(f"Maker must be of type {ExclusionMap}") - BkgMakerBase.bkg_map_maker.fset(self, maker) + self.__bkg_map_maker = maker class StackedExclusionMap(Stacked): """ @@ -914,8 +926,12 @@ def __init__(self, pointing_delta=self.pointing_delta ) - @BkgMakerBase.bkg_map_maker.setter + @property + def bkg_map_maker(self): + return self.__bkg_map_maker + + @bkg_map_maker.setter def bkg_map_maker(self, maker): if not isinstance(maker, ExclusionMap): raise TypeError(f"Maker must be of type {ExclusionMap}") - BkgMakerBase.bkg_map_maker.fset(self, maker) + self.__bkg_map_maker = maker From eab77ddb3f96dd96d61356a0aa3e7598ed80c718 Mon Sep 17 00:00:00 2001 From: Marcel Strzys Date: Wed, 31 Jul 2024 18:49:12 +0900 Subject: [PATCH 12/18] FIX: - simplified the reading of the data from the fits files with astropy.Tables follwing Max' suggestion - fixed the obstime object in the AltAz frame definition following Max' comment in the discussion of the PR --- src/pybkgmodel/data.py | 25 ++++++++++--------------- 1 file changed, 10 insertions(+), 15 deletions(-) diff --git a/src/pybkgmodel/data.py b/src/pybkgmodel/data.py index ed979dd..c12906e 100644 --- a/src/pybkgmodel/data.py +++ b/src/pybkgmodel/data.py @@ -4,12 +4,13 @@ import numpy import pandas import uproot -from astropy.io import fits -import astropy.time -import astropy.units as u from astropy.coordinates import SkyCoord, EarthLocation, AltAz from astropy.coordinates.erfa_astrom import erfa_astrom, ErfaAstromInterpolator +from astropy.io import fits +from astropy.table import Table +import astropy.time +import astropy.units as u def find_run_neighbours(target_run, run_list, time_delta, pointing_delta): @@ -645,28 +646,22 @@ def load_events(cls, file_name): try: evt_hdu = input_file["EVENTS"] evt_head = evt_hdu.header - evt_data = pandas.DataFrame(evt_hdu.data) + evt_data = Table.read(evt_hdu) event_data = {} for key, name in data_names_mapping.items(): - if key in evt_data.keys(): - event_data[name] = evt_data[key].to_numpy() - # Fits header keywords start counting from 1 - unit_name = f"TUNIT{evt_data.columns.get_loc(key)+1}" + event_data[name] = evt_data[key].quantity - try: - event_data[name] *= u.Unit(evt_head[unit_name]) - except KeyError: + if not evt_data[key].unit: event_data[name] *= u.one # Event times need to be converted from Instrument reference epoch ref_epoch = astropy.time.Time(evt_head['MJDREFI'], evt_head['MJDREFF'], format='mjd') - event_data['mjd'] = astropy.time.Time((evt_data['TIME'].quantity - + ref_epoch) - ).mjd * u.d + evt_time = evt_data['TIME'].quantity + ref_epoch + event_data['mjd'] = evt_time.mjd * u.d # TODO: current observatory location only La Palma, no mandatory header keyword obs_loc = EarthLocation(lat=28.761758*u.deg, @@ -675,7 +670,7 @@ def load_events(cls, file_name): if evt_head['OBS_MODE'] in ('POINTING', 'WOBBLE'): - alt_az_frame = AltAz(obstime=event_data['mjd'], + alt_az_frame = AltAz(obstime=evt_time, location=obs_loc) coords = SkyCoord(evt_head['RA_PNT'] *u.deg, From 98e48fc96c50fb28ddcdf2ac1f7a4fba46283ec8 Mon Sep 17 00:00:00 2001 From: Marcel Strzys Date: Wed, 31 Jul 2024 19:13:18 +0900 Subject: [PATCH 13/18] FIX: - fixed the setter function of the bkg_map_maker --- src/pybkgmodel/data.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/pybkgmodel/data.py b/src/pybkgmodel/data.py index c12906e..8234adf 100644 --- a/src/pybkgmodel/data.py +++ b/src/pybkgmodel/data.py @@ -577,8 +577,6 @@ def is_compatible(cls, file_name): ext = Path(file_name) - print("PATH", ext) - try: with fits.open(ext) as file: pass From 8966b98c8a618dca00888a0371c736769c55e1d8 Mon Sep 17 00:00:00 2001 From: Marcel Strzys Date: Wed, 31 Jul 2024 19:25:28 +0900 Subject: [PATCH 14/18] FIX: - debugging the setter for the bkg_map_maker - current version not running --- src/pybkgmodel/processing.py | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/src/pybkgmodel/processing.py b/src/pybkgmodel/processing.py index 66a5f4a..adac553 100644 --- a/src/pybkgmodel/processing.py +++ b/src/pybkgmodel/processing.py @@ -518,13 +518,13 @@ def __init__(self, @property def bkg_map_maker(self): - return self.__bkg_map_maker + return super().bkg_map_maker @bkg_map_maker.setter def bkg_map_maker(self, maker): if not isinstance(maker, WobbleMap): raise TypeError(f"Maker must be of type {WobbleMap}") - self.__bkg_map_maker = maker + super(RunwiseWobbleMap, type(self)).bkg_map_maker.__set__(self, maker) class StackedWobbleMap(Stacked): """ @@ -648,13 +648,13 @@ def __init__(self, @property def bkg_map_maker(self): - return self.__bkg_map_maker + return super().bkg_map_maker @bkg_map_maker.setter def bkg_map_maker(self, maker): if not isinstance(maker, WobbleMap): raise TypeError(f"Maker must be of type {WobbleMap}") - self.__bkg_map_maker = maker + super(StackedWobbleMap, type(self)).bkg_map_maker.__set__(self, maker) class RunwiseExclusionMap(Runwise): """ @@ -788,13 +788,17 @@ def __init__(self, @property def bkg_map_maker(self): - return self.__bkg_map_maker + print("Getter called. Current value:", self.__bkg_map_maker) + return super().bkg_map_maker @bkg_map_maker.setter def bkg_map_maker(self, maker): + print("Setter called with:", maker) if not isinstance(maker, ExclusionMap): raise TypeError(f"Maker must be of type {ExclusionMap}") - self.__bkg_map_maker = maker + super(RunwiseExclusionMap, type(self)).bkg_map_maker.__set__(self, maker) + print("Value set:", self.__bkg_map_maker) + class StackedExclusionMap(Stacked): """ @@ -928,10 +932,10 @@ def __init__(self, @property def bkg_map_maker(self): - return self.__bkg_map_maker + return super().bkg_map_maker @bkg_map_maker.setter def bkg_map_maker(self, maker): if not isinstance(maker, ExclusionMap): raise TypeError(f"Maker must be of type {ExclusionMap}") - self.__bkg_map_maker = maker + super(StackedExclusionMap, type(self)).bkg_map_maker.__set__(self, maker) From 84aeceb727bf5e5867456348d6b84ccd17fb12d3 Mon Sep 17 00:00:00 2001 From: Marcel Strzys Date: Wed, 31 Jul 2024 19:32:40 +0900 Subject: [PATCH 15/18] FIX: - intermediate commit for checking the updated setter --- src/pybkgmodel/processing.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/pybkgmodel/processing.py b/src/pybkgmodel/processing.py index adac553..23cb658 100644 --- a/src/pybkgmodel/processing.py +++ b/src/pybkgmodel/processing.py @@ -179,7 +179,7 @@ def __init__( self._bkg_maps = {} - self.__bkg_map_maker = None + self._bkg_map_maker = None @property @@ -187,12 +187,12 @@ def bkg_map_maker(self): """Getter for bkg_map_maker.""" print("This class uses the background method:", self.__bkg_map_maker.__class__.__name__) - return self.__bkg_map_maker + return self._bkg_map_maker @bkg_map_maker.setter def bkg_map_maker(self, value): """Setter for bkg_map_maker.""" - self.__bkg_map_maker = value + self._bkg_map_maker = value @property def bkg_maps(self): @@ -271,7 +271,7 @@ def generate_runwise_maps(self) -> dict: # Here the corrsponding bkg reconstruction algorith is applied # to obtain the runwise bkg map - bkg_map = self.__bkg_map_maker.get_runwise_bkg(target_run = run) + bkg_map = self._bkg_map_maker.get_runwise_bkg(target_run = run) # get corresponding names for the bkg maps under which they can # be safed @@ -788,7 +788,7 @@ def __init__(self, @property def bkg_map_maker(self): - print("Getter called. Current value:", self.__bkg_map_maker) + print("Getter called. Current value:", self._bkg_map_maker) return super().bkg_map_maker @bkg_map_maker.setter @@ -797,7 +797,7 @@ def bkg_map_maker(self, maker): if not isinstance(maker, ExclusionMap): raise TypeError(f"Maker must be of type {ExclusionMap}") super(RunwiseExclusionMap, type(self)).bkg_map_maker.__set__(self, maker) - print("Value set:", self.__bkg_map_maker) + print("Value set:", self._bkg_map_maker) class StackedExclusionMap(Stacked): From a5f98dde2b8b6fb1661dfc8aebee00a8804dc9ea Mon Sep 17 00:00:00 2001 From: Marcel Strzys Date: Wed, 31 Jul 2024 19:50:37 +0900 Subject: [PATCH 16/18] FIX: - ready for further review - setter for bkg_map_maker fixed - time scale added to ref_epoch - checked that the code runs with with LST-1 data --- src/pybkgmodel/data.py | 6 +++++- src/pybkgmodel/processing.py | 4 ---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/pybkgmodel/data.py b/src/pybkgmodel/data.py index 8234adf..661af9e 100644 --- a/src/pybkgmodel/data.py +++ b/src/pybkgmodel/data.py @@ -656,7 +656,11 @@ def load_events(cls, file_name): event_data[name] *= u.one # Event times need to be converted from Instrument reference epoch - ref_epoch = astropy.time.Time(evt_head['MJDREFI'], evt_head['MJDREFF'], format='mjd') + ref_epoch = astropy.time.Time(evt_head['MJDREFI'], + evt_head['MJDREFF'], + scale=evt_head['TIMESYS'].lower(), + format='mjd' + ) evt_time = evt_data['TIME'].quantity + ref_epoch event_data['mjd'] = evt_time.mjd * u.d diff --git a/src/pybkgmodel/processing.py b/src/pybkgmodel/processing.py index 23cb658..82515a4 100644 --- a/src/pybkgmodel/processing.py +++ b/src/pybkgmodel/processing.py @@ -788,18 +788,14 @@ def __init__(self, @property def bkg_map_maker(self): - print("Getter called. Current value:", self._bkg_map_maker) return super().bkg_map_maker @bkg_map_maker.setter def bkg_map_maker(self, maker): - print("Setter called with:", maker) if not isinstance(maker, ExclusionMap): raise TypeError(f"Maker must be of type {ExclusionMap}") super(RunwiseExclusionMap, type(self)).bkg_map_maker.__set__(self, maker) - print("Value set:", self._bkg_map_maker) - class StackedExclusionMap(Stacked): """ A class used to store the settings from the configuation file and to From af0e3a94c8618f144c6ce87c12d65f0428ebb52e Mon Sep 17 00:00:00 2001 From: Marcel Strzys Date: Mon, 5 Aug 2024 13:51:24 +0900 Subject: [PATCH 17/18] FIX: - changed the remaining numpy's to np, not covered by the merger --- src/pybkgmodel/data.py | 37 +++++++++++++++++++++---------------- 1 file changed, 21 insertions(+), 16 deletions(-) diff --git a/src/pybkgmodel/data.py b/src/pybkgmodel/data.py index 1e5d1a7..f65ec32 100644 --- a/src/pybkgmodel/data.py +++ b/src/pybkgmodel/data.py @@ -494,20 +494,25 @@ def load_events(cls, file_name, cuts): is_simulated = is_mc and 'trigger_time' in data if not is_mc or is_simulated: - event_data['mjd'] = astropy.time.Time(data['trigger_time'].to_numpy(), format='unix').mjd + event_data['mjd'] = astropy.time.Time(data['trigger_time'].to_numpy(), + format='unix').mjd lst_time = astropy.time.Time(event_data['mjd'], format='mjd') lst_loc = EarthLocation(lat=28.761758*u.deg, lon=-17.890659*u.deg, height=2200*u.m) alt_az_frame = AltAz(obstime=lst_time, location=lst_loc) if 'pointing_ra' not in event_data: - coords = SkyCoord(alt=data['alt_tel'].to_numpy()*u.rad, az=data['az_tel'].to_numpy()*u.rad, frame=alt_az_frame).icrs + coords = SkyCoord(alt=data['alt_tel'].to_numpy()*u.rad, + az=data['az_tel'].to_numpy()*u.rad, + frame=alt_az_frame).icrs event_data['pointing_ra'] = coords.ra.to(data_units['pointing_ra']).value event_data['pointing_dec'] = coords.dec.to(data_units['pointing_dec']).value if 'event_ra' not in event_data: - coords = SkyCoord(alt=data['reco_alt'].to_numpy()*u.rad, az=data['reco_az'].to_numpy()*u.rad, frame=alt_az_frame).icrs + coords = SkyCoord(alt=data['reco_alt'].to_numpy()*u.rad, + az=data['reco_az'].to_numpy()*u.rad, + frame=alt_az_frame).icrs event_data['event_ra'] = coords.ra.to(data_units['event_ra']).value event_data['event_dec'] = coords.dec.to(data_units['event_dec']).value @@ -576,14 +581,14 @@ def is_compatible(cls, file_name): """ ext = Path(file_name) - + try: with fits.open(ext) as file: pass compatible = True except OSError: compatible = False - + return compatible @classmethod @@ -656,8 +661,8 @@ def load_events(cls, file_name): event_data[name] *= u.one # Event times need to be converted from Instrument reference epoch - ref_epoch = astropy.time.Time(evt_head['MJDREFI'], - evt_head['MJDREFF'], + ref_epoch = astropy.time.Time(evt_head['MJDREFI'], + evt_head['MJDREFF'], scale=evt_head['TIMESYS'].lower(), format='mjd' ) @@ -685,11 +690,11 @@ def load_events(cls, file_name): event_data['pointing_az'] = altaz_pointing.az.to(u.deg) - event_data['pointing_ra'] = numpy.array( + event_data['pointing_ra'] = np.array( [evt_head['RA_PNT']] \ * len(event_data['pointing_zd']) ) * u.deg - event_data['pointing_dec'] = numpy.array( + event_data['pointing_dec'] = np.array( [evt_head['DEC_PNT']] \ * len(event_data['pointing_zd']) ) * u.deg @@ -716,17 +721,17 @@ def load_events(cls, file_name): Function will return zeros for the event location.\ Supported modes: POINTING, DRIFT") - event_data['pointing_zd'] = numpy.zeros(len(event_data['mjd'])) * u.deg - event_data['pointing_az'] = numpy.zeros(len(event_data['mjd'])) * u.deg - event_data['pointing_ra'] = numpy.zeros(len(event_data['mjd'])) * u.deg - event_data['pointing_dec'] = numpy.zeros(len(event_data['mjd'])) * u.deg + event_data['pointing_zd'] = np.zeros(len(event_data['mjd'])) * u.deg + event_data['pointing_az'] = np.zeros(len(event_data['mjd'])) * u.deg + event_data['pointing_ra'] = np.zeros(len(event_data['mjd'])) * u.deg + event_data['pointing_dec'] = np.zeros(len(event_data['mjd'])) * u.deg except KeyError: print(f"File {file_name} corrupted or missing the Events hdu." + "Empty arrays will be returned.") - finite = [numpy.isfinite(item) for key, item in event_data.items() if item is not None] - all_finite = numpy.prod(finite, axis=0, dtype=bool) + finite = [np.isfinite(item) for key, item in event_data.items() if item is not None] + all_finite = np.prod(finite, axis=0, dtype=bool) for key, item in event_data.items(): if item is not None: @@ -742,7 +747,7 @@ def load_events(cls, file_name): event_data['pointing_zd'], event_data['mjd'], None, - numpy.array(evt_head['LIVETIME']) * u.s + np.array(evt_head['LIVETIME']) * u.s ) return event_sample From 33ea973005bd9e2a7048aff32e0707efbe8645ed Mon Sep 17 00:00:00 2001 From: Marcel Strzys Date: Mon, 5 Aug 2024 14:59:35 +0900 Subject: [PATCH 18/18] FIX: - removed the np.bool for bool inherited from the main branch - tested with LST1 data and ensure that it runs properly --- src/pybkgmodel/camera.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pybkgmodel/camera.py b/src/pybkgmodel/camera.py index 79efe12..57b4128 100644 --- a/src/pybkgmodel/camera.py +++ b/src/pybkgmodel/camera.py @@ -169,7 +169,7 @@ def __init__(self, counts, xedges, yedges, energy_edges, center=None, mask=None, ny = yedges.size - 1 if mask is None: - mask = np.ones((nx, ny), dtype=np.bool) + mask = np.ones((nx, ny), dtype=bool) if exposure is None: exposure = np.ones((nx, ny), dtype=np.float) * u.s @@ -329,7 +329,7 @@ def mask_region(self, region): def mask_reset(self): """_summary_ """ - self.mask = np.ones((self.xedges.size - 1, self.yedges.size - 1), dtype=np.bool) + self.mask = np.ones((self.xedges.size - 1, self.yedges.size - 1), dtype=bool) def plot(self, energy_bin_id=0, ax_unit='deg', val_unit='1/s', **kwargs):