From e14118b7f0d1c0dc644da62f77d0d9958649fd25 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Mon, 30 Mar 2026 18:05:41 +0200 Subject: [PATCH 1/2] Intermediary output for models --- baccmod/base_acceptance_map_creator.py | 69 +++++++++++++++++++++--- baccmod/grid3d_acceptance_map_creator.py | 15 +++++- baccmod/radial_acceptance_map_creator.py | 14 ++++- test/integration_test.py | 18 +++---- 4 files changed, 97 insertions(+), 19 deletions(-) diff --git a/baccmod/base_acceptance_map_creator.py b/baccmod/base_acceptance_map_creator.py index f15fe4e..4083045 100644 --- a/baccmod/base_acceptance_map_creator.py +++ b/baccmod/base_acceptance_map_creator.py @@ -180,7 +180,7 @@ def __init__(self, self.mini_irf_time_resolution = mini_irf_time_resolution @abstractmethod - def create_model(self, observations: Observations) -> BackgroundIRF: + def create_model(self, observations: Observations, add_intermediary_output: bool = False) -> BackgroundIRF: """ Abstract method to calculate an acceptance map from a list of observations. @@ -190,6 +190,8 @@ def create_model(self, observations: Observations) -> BackgroundIRF: ---------- observations : gammapy.data.observations.Observations The collection of observations used to create the acceptance map. + add_intermediary_output: bool + Add for each model a baccmod field in meta with all the intermediary state of the computation Returns ------- @@ -732,7 +734,8 @@ def _cluster_observations(self, return observations_dict, split_obs def _create_model_cos_zenith_binned(self, - observations: Observations + observations: Observations, + add_intermediary_output: bool = False, ) -> BackgroundCollectionZenith: """ Calculate a model for each cos zenith bin @@ -741,6 +744,8 @@ def _create_model_cos_zenith_binned(self, ---------- observations : gammapy.data.observations.Observations The collection of observations used to make the background model + add_intermediary_output: bool + Add for each model a baccmod field in meta with all the intermediary state of the computation Returns ------- @@ -877,7 +882,7 @@ def _create_model_cos_zenith_binned(self, binned_model = [] for i, binned_obs in enumerate(binned_observations): logger.info(f"Creating model for the bin at cos zenith = {np.round(bin_center[i], 2)}°") - binned_model.append(self.create_model(binned_obs)) + binned_model.append(self.create_model(binned_obs, add_intermediary_output)) dict_binned_model = {} for i in range(len(binned_model)): @@ -894,7 +899,8 @@ def _create_model_cos_zenith_binned(self, def _create_acceptance_model(self, off_observations: dict[str, Observations], zenith_binning: bool = False, - zenith_interpolation: bool = False + zenith_interpolation: bool = False, + add_intermediary_output: bool = False, ) -> Union[BackgroundIRF,BackgroundCollectionZenith]: """ A function to create a combined background model for observations previously separated in subsets. @@ -903,6 +909,8 @@ def _create_acceptance_model(self, off_observations: gammapy.data.observations.Observations zenith_binning: bool zenith_interpolation: bool + add_intermediary_output: bool + Add for each model a baccmod field in meta with all the intermediary state of the computation Returns ------- @@ -911,9 +919,9 @@ def _create_acceptance_model(self, models={} for key in off_observations.keys(): if zenith_interpolation or zenith_binning: - models[key] = self._create_model_cos_zenith_binned(observations=off_observations[key]) + models[key] = self._create_model_cos_zenith_binned(observations=off_observations[key], add_intermediary_output=add_intermediary_output) else: - models[key] = self.create_model(observations=off_observations[key]) + models[key] = self.create_model(observations=off_observations[key], add_intermediary_output=add_intermediary_output) if self.azimuth_east_west_splitting: model=BackgroundCollectionZenithSplitAzimuth(bkg_east=models['east'], bkg_west=models['west'], interpolation_type=self.interpolation_type, @@ -1258,7 +1266,8 @@ def _normalise_model_per_run(self, def create_acceptance_model(self, off_observations: Observations, - zenith_binning: bool = False + zenith_binning: bool = False, + add_intermediary_output: bool = False, ) -> Union[BackgroundIRF,BackgroundCollectionZenith]: """ A high level function to create a background model from observations. @@ -1266,6 +1275,8 @@ def create_acceptance_model(self, ---------- off_observations: gammapy.data.observations.Observations zenith_binning: bool + add_intermediary_output: bool + Add for each model a baccmod field in meta with all the intermediary state of the computation Returns ------- @@ -1273,7 +1284,9 @@ def create_acceptance_model(self, """ off_observations_sets, _ = self._cluster_observations(off_observations) - return self._create_acceptance_model(off_observations_sets, zenith_binning) + return self._create_acceptance_model(off_observations_sets, + zenith_binning, + add_intermediary_output=add_intermediary_output) def create_acceptance_map_per_observation(self, @@ -1338,3 +1351,43 @@ def create_acceptance_map_per_observation(self, acceptance_map = self._normalise_model_per_run(observations, acceptance_map) return acceptance_map + + @staticmethod + def meta_observations(observations: Observations) -> List[Dict]: + """ + Return the meta block with the information detail of which observation where used to create the model + + Parameters + ---------- + observations : gammapy.data.observations.Observations + The collection of observations used to make the acceptance map + + Returns + ------- + meta : List of Dict + """ + + meta = [] + + for obs in observations: + meta.append({ + "obs_id": obs.obs_id, + "obs_tstart": obs.tstart, + "obs_tstop": obs.tstop, + "obs_mid_alt_az": obs.get_pointing_altaz(obs.tmid), + "obs_pointing": obs.get_pointing_icrs(obs.tmid) + }) + return meta + + def meta_exclusion_region(self) -> List[Dict]: + """ + Return the meta block with the information detail of which exclusion regions where used to create the model + + Returns + ------- + meta : List of Dict + """ + + if self.exclude_regions is None: + return [] + return copy.deepcopy(self.exclude_regions) \ No newline at end of file diff --git a/baccmod/grid3d_acceptance_map_creator.py b/baccmod/grid3d_acceptance_map_creator.py index 1ca57d3..996d207 100644 --- a/baccmod/grid3d_acceptance_map_creator.py +++ b/baccmod/grid3d_acceptance_map_creator.py @@ -165,7 +165,7 @@ def f(*args): return fnc(x, y, **m.values.to_dict()) - def create_model(self, observations: Observations) -> Background3D: + def create_model(self, observations: Observations, add_intermediary_output: bool = False) -> Background3D: """ Calculate a 3D grid acceptance map @@ -173,6 +173,8 @@ def create_model(self, observations: Observations) -> Background3D: ---------- observations : gammapy.data.observations.Observations The collection of observations used to make the acceptance map + add_intermediary_output: bool + Add for each model a meta_baccmod field with all the intermediary state of the computation Returns ------- @@ -227,6 +229,17 @@ def create_model(self, observations: Observations) -> Background3D: data=data_background.to(u.Unit('s-1 MeV-1 sr-1')), fov_alignment=FoVAlignment.ALTAZ) + + if add_intermediary_output: + meta = {'observations': self.meta_observations(observations), + 'exclusion_region': self.meta_exclusion_region(), + 'count_map': count_background, + 'corrected_exposure': exp_map_background, + 'total_exposure': exp_map_background_total, + 'energy_axis_computation': energy_axis_computation, + 'background_model_computation': data_background} + acceptance_map.meta['baccmod'] = meta + return acceptance_map def _create_base_computation_map(self, observations: Observations) -> Tuple[np.ndarray, WcsNDMap, WcsNDMap, u.Quantity, MapAxis]: diff --git a/baccmod/radial_acceptance_map_creator.py b/baccmod/radial_acceptance_map_creator.py index 761f66d..9230e4b 100644 --- a/baccmod/radial_acceptance_map_creator.py +++ b/baccmod/radial_acceptance_map_creator.py @@ -65,7 +65,7 @@ def __init__(self, exclude_regions=exclude_regions, **kwargs) - def create_model(self, observations: Observations) -> Background2D: + def create_model(self, observations: Observations, add_intermediary_output: bool = False) -> Background2D: """ Calculate a radial acceptance map @@ -73,6 +73,8 @@ def create_model(self, observations: Observations) -> Background2D: ---------- observations : Observations The collection of observations used to make the acceptance map + add_intermediary_output: bool + Add for each model a baccmod field in meta with all the intermediary state of the computation Returns ------- @@ -104,6 +106,16 @@ def create_model(self, observations: Observations) -> Background2D: acceptance_map = Background2D(axes=[self.energy_axis, self.offset_axis], data=self._interpolate_bkg_to_energy_axis(data_background, energy_axis_computation)) + if add_intermediary_output: + meta = {'observations': self.meta_observations(observations), + 'exclusion_region': self.meta_exclusion_region(), + 'count_map': count_map_background, + 'corrected_exposure': exp_map_background, + 'total_exposure': exp_map_background_total, + 'energy_axis_computation': energy_axis_computation, + 'background_model_computation': data_background} + acceptance_map.meta['baccmod'] = meta + return acceptance_map def _create_base_computation_map(self, observations: Observations) -> Tuple[np.ndarray, WcsNDMap, WcsNDMap, u.Quantity, MapAxis]: diff --git a/test/integration_test.py b/test/integration_test.py index 13572d9..2097fe2 100644 --- a/test/integration_test.py +++ b/test/integration_test.py @@ -63,7 +63,7 @@ def test_integration_3D(self): offset_axis=self.offset_axis, oversample_map=5, exclude_regions=self.exclude_region_PKS_2155) - background_model = bkg_maker.create_model(observations=self.obs_collection_pks_2155) + background_model = bkg_maker.create_model(observations=self.obs_collection_pks_2155, add_intermediary_output=True) assert type(background_model) is Background3D reference = Background3D.read('ressource/test_data/reference_model/pks_2155_3D.fits') self._print_model_precision(background_model, reference) @@ -77,7 +77,7 @@ def test_integration_spatial_fit(self): oversample_map=5, exclude_regions=self.exclude_region_PKS_2155, method='fit') - background_model = bkg_maker.create_model(observations=self.obs_collection_pks_2155) + background_model = bkg_maker.create_model(observations=self.obs_collection_pks_2155, add_intermediary_output=True) assert type(background_model) is Background3D reference = Background3D.read('ressource/test_data/reference_model/pks_2155_spatial_fit_bkg.fits') self._print_model_precision(background_model, reference) @@ -90,7 +90,7 @@ def test_integration_2D(self): offset_axis=self.offset_axis, oversample_map=5, exclude_regions=self.exclude_region_PKS_2155) - background_model = bkg_maker.create_model(observations=self.obs_collection_pks_2155) + background_model = bkg_maker.create_model(observations=self.obs_collection_pks_2155, add_intermediary_output=True) assert type(background_model) is Background2D reference = Background2D.read('ressource/test_data/reference_model/pks_2155_2D.fits') self._print_model_precision(background_model, reference) @@ -104,7 +104,7 @@ def test_integration_3D_irregular_computation_axis(self): offset_axis=self.offset_axis, oversample_map=5, exclude_regions=self.exclude_region_PKS_2155) - background_model = bkg_maker.create_model(observations=self.obs_collection_pks_2155) + background_model = bkg_maker.create_model(observations=self.obs_collection_pks_2155, add_intermediary_output=True) assert type(background_model) is Background3D reference = Background3D.read('ressource/test_data/reference_model/pks_2155_3D_bkg_irregular_energy.fits') self._print_model_precision(background_model, reference) @@ -119,7 +119,7 @@ def test_integration_spatial_fit_irregular_computation_axis(self): oversample_map=5, exclude_regions=self.exclude_region_PKS_2155, method='fit') - background_model = bkg_maker.create_model(observations=self.obs_collection_pks_2155) + background_model = bkg_maker.create_model(observations=self.obs_collection_pks_2155, add_intermediary_output=True) assert type(background_model) is Background3D reference = Background3D.read('ressource/test_data/reference_model/pks_2155_spatial_fit_bkg_irregular_energy.fits') self._print_model_precision(background_model, reference) @@ -133,7 +133,7 @@ def test_integration_2D_irregular_computation_axis(self): offset_axis=self.offset_axis, oversample_map=5, exclude_regions=self.exclude_region_PKS_2155) - background_model = bkg_maker.create_model(observations=self.obs_collection_pks_2155) + background_model = bkg_maker.create_model(observations=self.obs_collection_pks_2155, add_intermediary_output=True) assert type(background_model) is Background2D reference = Background2D.read('ressource/test_data/reference_model/pks_2155_2D_bkg_irregular_energy.fits') self._print_model_precision(background_model, reference) @@ -151,7 +151,7 @@ def test_integration_3D_dynamic_irregular_computation_axis(self): dynamic_energy_axis_target_statistics=100, oversample_map=5, exclude_regions=self.exclude_region_PKS_2155) - background_model = bkg_maker.create_model(observations=self.obs_collection_pks_2155) + background_model = bkg_maker.create_model(observations=self.obs_collection_pks_2155, add_intermediary_output=True) assert type(background_model) is Background3D reference = Background3D.read('ressource/test_data/reference_model/pks_2155_3D_bkg_dynamic_irregular_energy.fits') self._print_model_precision(background_model, reference) @@ -168,7 +168,7 @@ def test_integration_spatial_fit_dynamic_irregular_computation_axis(self): oversample_map=5, exclude_regions=self.exclude_region_PKS_2155, method='fit') - background_model = bkg_maker.create_model(observations=self.obs_collection_pks_2155) + background_model = bkg_maker.create_model(observations=self.obs_collection_pks_2155, add_intermediary_output=True) assert type(background_model) is Background3D reference = Background3D.read('ressource/test_data/reference_model/pks_2155_spatial_fit_bkg_dynamic_irregular_energy.fits') self._print_model_precision(background_model, reference) @@ -184,7 +184,7 @@ def test_integration_2D_dynamic_irregular_computation_axis(self): dynamic_energy_axis_target_statistics=500, oversample_map=5, exclude_regions=self.exclude_region_PKS_2155) - background_model = bkg_maker.create_model(observations=self.obs_collection_pks_2155) + background_model = bkg_maker.create_model(observations=self.obs_collection_pks_2155, add_intermediary_output=True) assert type(background_model) is Background2D reference = Background2D.read('ressource/test_data/reference_model/pks_2155_2D_bkg_dynamic_irregular_energy.fits') self._print_model_precision(background_model, reference) From 288c6f8f97aee5325a42cf3018f089e48525b286 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Tue, 21 Apr 2026 11:19:47 +0200 Subject: [PATCH 2/2] Fix --- baccmod/grid3d_acceptance_map_creator.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/baccmod/grid3d_acceptance_map_creator.py b/baccmod/grid3d_acceptance_map_creator.py index 996d207..68bc333 100644 --- a/baccmod/grid3d_acceptance_map_creator.py +++ b/baccmod/grid3d_acceptance_map_creator.py @@ -223,10 +223,10 @@ def create_model(self, observations: Observations, add_intermediary_output: bool data_background = corrected_counts / solid_angle[np.newaxis, :, :] / energy_axis_computation.bin_width[:, np.newaxis, np.newaxis] / livetime - data_background = self._interpolate_bkg_to_energy_axis(data_background, energy_axis_computation) + data_background_interpolated = self._interpolate_bkg_to_energy_axis(data_background, energy_axis_computation) acceptance_map = Background3D(axes=[self.energy_axis, extended_offset_axis_x, extended_offset_axis_y], - data=data_background.to(u.Unit('s-1 MeV-1 sr-1')), + data=data_background_interpolated.to(u.Unit('s-1 MeV-1 sr-1')), fov_alignment=FoVAlignment.ALTAZ)