From f4b2f360cc6bafb3089bee2f306bb21993a55233 Mon Sep 17 00:00:00 2001 From: Gabriel Emery Date: Tue, 28 Apr 2026 15:59:06 +0200 Subject: [PATCH 1/2] Implement center point integration in the fitting procedure. Fix support for positionnal arguments in fitting class. --- baccmod/base_acceptance_map_creator.py | 2 +- baccmod/fit_acceptance_map_creator.py | 53 +++++++++++++++----------- baccmod/fitting.py | 10 +++-- 3 files changed, 39 insertions(+), 26 deletions(-) diff --git a/baccmod/base_acceptance_map_creator.py b/baccmod/base_acceptance_map_creator.py index 93ffc3c..6208efb 100644 --- a/baccmod/base_acceptance_map_creator.py +++ b/baccmod/base_acceptance_map_creator.py @@ -1135,7 +1135,7 @@ def _interpolate_bkg_to_energy_axis(self, data_bkg: u.Quantity, energy_axis_comp """ # Return the provided bkg data if energy axis are matching - if len(energy_axis_computation.edges) == len(self.energy_axis.edges) and np.all(energy_axis_computation.edges == self.energy_axis.edges): + if energy_axis_computation.nbin == self.energy_axis.nbin and np.all(energy_axis_computation.edges == self.energy_axis.edges): logger.info('Identical computation energy axis and model energy axis, no interpolation required') return data_bkg diff --git a/baccmod/fit_acceptance_map_creator.py b/baccmod/fit_acceptance_map_creator.py index ffd936b..47956a3 100644 --- a/baccmod/fit_acceptance_map_creator.py +++ b/baccmod/fit_acceptance_map_creator.py @@ -31,6 +31,8 @@ class FitAcceptanceMapCreator(Grid3DAcceptanceMapCreator): def __init__( self, + energy_axis: MapAxis, + offset_axis: MapAxis, model_to_fit: FittableModel = Gaussian2D(), list_name_normalisation_parameter: List[str] = ['amplitude'], maxiter: int = 1000, @@ -56,12 +58,10 @@ def __init__( """ # Call the “stack”‐only constructor in Grid3D, to set up geometry, offset axes, etc. - super().__init__(**kwargs) + super().__init__(energy_axis, offset_axis, **kwargs) self.sq_rel_residuals = {"mean": [], "std": []} self.model_to_fit = model_to_fit - if model_to_fit.n_inputs !=2: #TODO remove once 1D and 3D models are tested - logger.warning('Only 2D models, fitted per energy bin, were thoroughly tested.') self.list_name_normalisation_parameter = list_name_normalisation_parameter self.maxiter = maxiter @@ -89,8 +89,19 @@ def create_model(self, observations) -> Background3D: exp_ds = exp_map_background.downsample(self.oversample_map, preserve_counts=True) exp_total_ds = exp_map_background_total.downsample(self.oversample_map, preserve_counts=True) + # 3) build final offset axes (matching Grid3DAcceptanceMapCreator) + edges = self.offset_axis.edges + extended_edges = np.concatenate((-np.flip(edges), edges[1:]), axis=None) + extended_offset_axis_x = MapAxis.from_edges(extended_edges, name="fov_lon") + extended_offset_axis_y = MapAxis.from_edges(extended_edges, name="fov_lat") + + bin_width_x = np.repeat(extended_offset_axis_x.bin_width[:, np.newaxis], extended_offset_axis_x.nbin, axis=1) + bin_width_y = np.repeat(extended_offset_axis_y.bin_width[np.newaxis, :], extended_offset_axis_y.nbin, axis=0) + solid_angle = 4.0 * (np.sin(bin_width_x / 2) * np.sin(bin_width_y / 2)) * u.steradian + + energy_bin_width = energy_axis_computation.bin_width - # 3) fit function on counts → “corrected counts” + # 4) fit function on counts → “corrected counts” predicted_counts = np.empty(count_background.shape) self.sq_rel_residuals = {"mean": [], "std": []} @@ -108,6 +119,7 @@ def create_model(self, observations) -> Background3D: # perform the fit, looping over non fitted axes if self.model_to_fit.n_inputs == 3: + bin_size = (solid_angle[np.newaxis, :, :] * energy_bin_width[:, np.newaxis, np.newaxis]).value logger.info( "Fitting background with a 3D model." ) @@ -117,8 +129,10 @@ def create_model(self, observations) -> Background3D: count_map=count_background.astype(int), exp_map_total=exp_total_ds.data, exp_map=exp_ds.data, + bin_size=bin_size ) elif self.model_to_fit.n_inputs == 2: + bin_size = solid_angle.value logger.info( "Fitting background per enery bin" ) @@ -134,8 +148,10 @@ def create_model(self, observations) -> Background3D: count_map=count_background[e].astype(int), exp_map_total=exp_total_ds.data[e], exp_map=exp_ds.data[e], + bin_size=bin_size ) elif self.model_to_fit.n_inputs == 1: + bin_size = energy_bin_width.value logger.info( "Fitting background per spatial bin" ) @@ -150,6 +166,7 @@ def create_model(self, observations) -> Background3D: count_map=count_background[:,x,y].astype(int), exp_map_total=exp_total_ds.data[:,x,y], exp_map=exp_ds.data[:,x,y], + bin_size=bin_size ) else: raise RuntimeError(f"The provided model dimension is incorrect : {self.model_to_fit.n_inputs}") @@ -160,16 +177,6 @@ def create_model(self, observations) -> Background3D: np.array_str(np.round(self.sq_rel_residuals['std'], 2)) ) - # 4) build final offset axes (matching Grid3DAcceptanceMapCreator) - edges = self.offset_axis.edges - extended_edges = np.concatenate((-np.flip(edges), edges[1:]), axis=None) - extended_offset_axis_x = MapAxis.from_edges(extended_edges, name="fov_lon") - extended_offset_axis_y = MapAxis.from_edges(extended_edges, name="fov_lat") - - bin_width_x = np.repeat(extended_offset_axis_x.bin_width[:, np.newaxis], extended_offset_axis_x.nbin, axis=1) - bin_width_y = np.repeat(extended_offset_axis_y.bin_width[np.newaxis, :], extended_offset_axis_y.nbin, axis=0) - solid_angle = 4.0 * (np.sin(bin_width_x / 2) * np.sin(bin_width_y / 2)) * u.steradian - # 5) normalize to flux units data_background = ( predicted_counts @@ -186,7 +193,7 @@ def create_model(self, observations) -> Background3D: return acceptance_map def _fit_background(self, model: Model, *coords: np.ndarray, count_map: np.ndarray, - exp_map_total: np.ndarray, exp_map: np.ndarray) -> np.ndarray: + exp_map_total: np.ndarray, exp_map: np.ndarray, bin_size: np.ndarray) -> np.ndarray: """ Perform a Poisson fit on the given map (could be 1D, 2D or 3D), given the fine‐binned exposure (exp_map) and total‐exposure (exp_map_total), @@ -199,16 +206,18 @@ def _fit_background(self, model: Model, *coords: np.ndarray, count_map: np.ndarr coords : np.ndarray the list of coordinates count_map : np.ndarray - counts in each pixel, integer + counts in each bin, integer exp_map_total : np.ndarray exposure WITHOUT exclusion regions exp_map : np.ndarray exposure CORRECTED for exclusion regions + bin_size : np.ndarray + integral size of each bin Returns ------- - fitted_counts_model : 2D np.ndarray - The best‐fit model of counts (on the same pixel grid as count_map). + fitted_counts_model : np.ndarray + The best‐fit model of counts (on the same bin grid as count_map). """ total_counts = np.sum(count_map) # Skip the fit if there are no events @@ -230,7 +239,7 @@ def _fit_background(self, model: Model, *coords: np.ndarray, count_map: np.ndarr # Correct normalisation of the model if self.list_name_normalisation_parameter is not None and len(self.list_name_normalisation_parameter) > 0: # Compute correction - init_count_model = np.sum(model_init(*coords)*exp_correction) + init_count_model = np.sum(model_init(*coords)*exp_correction*bin_size) correction_norm = total_counts/init_count_model # build list of free parameters @@ -248,10 +257,10 @@ def _fit_background(self, model: Model, *coords: np.ndarray, count_map: np.ndarr # Fit the model best_model = poisson_fitter(model_init, *coords, data=count_map, exposure_correction=exp_correction, - mask=mask, maxiter=self.maxiter) + mask=mask, maxiter=self.maxiter, bin_size=bin_size) # Collect results & (optionally) track residuals if logger.getEffectiveLevel() <= logging.INFO: - fitted_model = best_model(*coords) * exp_correction + fitted_model = best_model(*coords) * exp_correction * bin_size fitted_model[exp_map == 0] = 1.0 rel_resid = 100 * (count_map - fitted_model) / fitted_model sq_rel_resid = (count_map - fitted_model) / np.sqrt(fitted_model) @@ -263,4 +272,4 @@ def _fit_background(self, model: Model, *coords: np.ndarray, count_map: np.ndarr "Avg rel residual: %.1f, Std = %.2f\n", np.mean(rel_resid), np.std(rel_resid) ) - return best_model(*coords) + return best_model(*coords) * bin_size diff --git a/baccmod/fitting.py b/baccmod/fitting.py index da4f31d..23fc43e 100644 --- a/baccmod/fitting.py +++ b/baccmod/fitting.py @@ -21,7 +21,8 @@ def poisson_fitter(model: Model, data: np.ndarray, exposure_correction: np.ndarray = None, mask: np.ndarray = None, - maxiter: int = 1000) -> Model: + maxiter: int = 1000, + bin_size: np.ndarray = 1) -> Model: """ Fit the model to the data following poisson likelihood statistics and using iminuit @@ -40,6 +41,8 @@ def poisson_fitter(model: Model, maxiter : int maximum number of iteration for fitting, as the fitting is performed in two steps, could be the double of this value in practice + bin_size : np.ndarray + integral size of each bin Returns ------- @@ -88,7 +91,7 @@ def apply_params_and_tied(pars): # negative log‑likelihood def neg_logL(**pars): apply_params_and_tied(pars) - mu = model_copy(*coords) * exposure_correction + mu = model_copy(*coords) * exposure_correction * bin_size return -np.sum(log_poisson(mu[mask], data[mask], log_fact[mask])) # wrapper to accept either positional arguments or keyword arguments (exclusives) @@ -139,7 +142,8 @@ def log_factorial(x: np.ndarray) -> np.ndarray: def log_poisson(mu: np.ndarray, x: np.ndarray, log_factorial_x: np.ndarray) -> np.ndarray: + """ Poisson pdf in log scale. """ if np.any(mu <= 0): logger.warning('log poisson received a zero or negative value.') - mu[mu <= 0] = np.min(mu[mu > 0])/2 + mu[mu <= 0] = np.finfo(mu.dtype).tiny return -mu + x * np.log(mu) - log_factorial_x From 29682989218cb08a99a26c9d2c77dc34186a6940 Mon Sep 17 00:00:00 2001 From: Gabriel Emery Date: Mon, 15 Jun 2026 15:21:47 +0200 Subject: [PATCH 2/2] Refactor spatial axis manipulation in Grid3DAcceptanceMapCreator --- baccmod/fit_acceptance_map_creator.py | 12 ++---------- baccmod/grid3d_acceptance_map_creator.py | 23 +++++++++++++++-------- 2 files changed, 17 insertions(+), 18 deletions(-) diff --git a/baccmod/fit_acceptance_map_creator.py b/baccmod/fit_acceptance_map_creator.py index 47956a3..7c5b235 100644 --- a/baccmod/fit_acceptance_map_creator.py +++ b/baccmod/fit_acceptance_map_creator.py @@ -89,16 +89,8 @@ def create_model(self, observations) -> Background3D: exp_ds = exp_map_background.downsample(self.oversample_map, preserve_counts=True) exp_total_ds = exp_map_background_total.downsample(self.oversample_map, preserve_counts=True) - # 3) build final offset axes (matching Grid3DAcceptanceMapCreator) - edges = self.offset_axis.edges - extended_edges = np.concatenate((-np.flip(edges), edges[1:]), axis=None) - extended_offset_axis_x = MapAxis.from_edges(extended_edges, name="fov_lon") - extended_offset_axis_y = MapAxis.from_edges(extended_edges, name="fov_lat") - - bin_width_x = np.repeat(extended_offset_axis_x.bin_width[:, np.newaxis], extended_offset_axis_x.nbin, axis=1) - bin_width_y = np.repeat(extended_offset_axis_y.bin_width[np.newaxis, :], extended_offset_axis_y.nbin, axis=0) - solid_angle = 4.0 * (np.sin(bin_width_x / 2) * np.sin(bin_width_y / 2)) * u.steradian - + # 3) Get spatial axis for the bkg model and the bins size + extended_offset_axis_x, extended_offset_axis_y, solid_angle = self._get_ext_axis_and_solid_angle() energy_bin_width = energy_axis_computation.bin_width # 4) fit function on counts → “corrected counts” diff --git a/baccmod/grid3d_acceptance_map_creator.py b/baccmod/grid3d_acceptance_map_creator.py index b305410..5d6d7f1 100644 --- a/baccmod/grid3d_acceptance_map_creator.py +++ b/baccmod/grid3d_acceptance_map_creator.py @@ -81,6 +81,19 @@ def __init__(self, exclude_regions=exclude_regions, **kwargs) + def _get_ext_axis_and_solid_angle(self): + """ Compute the solid angle of spatial bins """ + edges = self.offset_axis.edges + extended_edges = np.concatenate((-np.flip(edges), edges[1:]), axis=None) + extended_offset_axis_x = MapAxis.from_edges(extended_edges, name="fov_lon") + extended_offset_axis_y = MapAxis.from_edges(extended_edges, name="fov_lat") + + bin_width_x = np.repeat(extended_offset_axis_x.bin_width[:, np.newaxis], extended_offset_axis_x.nbin, axis=1) + bin_width_y = np.repeat(extended_offset_axis_y.bin_width[np.newaxis, :], extended_offset_axis_y.nbin, axis=0) + solid_angle = 4.0 * (np.sin(bin_width_x / 2) * np.sin(bin_width_y / 2)) * u.steradian + + return extended_offset_axis_x, extended_offset_axis_y, solid_angle + def create_model(self, observations: Observations) -> Background3D: """ Calculate a 3D grid acceptance map @@ -104,18 +117,12 @@ def create_model(self, observations: Observations) -> Background3D: exp_map_background_total_downsample = exp_map_background_total.downsample(self.oversample_map, preserve_counts=True) - # Create axis for bkg model - edges = self.offset_axis.edges - extended_edges = np.concatenate((-np.flip(edges), edges[1:]), axis=None) - extended_offset_axis_x = MapAxis.from_edges(extended_edges, name='fov_lon') - bin_width_x = np.repeat(extended_offset_axis_x.bin_width[:, np.newaxis], extended_offset_axis_x.nbin, axis=1) - extended_offset_axis_y = MapAxis.from_edges(extended_edges, name='fov_lat') - bin_width_y = np.repeat(extended_offset_axis_y.bin_width[np.newaxis, :], extended_offset_axis_y.nbin, axis=0) + # Create spatial axis for the bkg model and the bin solid angle + extended_offset_axis_x, extended_offset_axis_y, solid_angle = self._get_ext_axis_and_solid_angle() # Compute acceptance_map corrected_counts = count_background * (exp_map_background_total_downsample.data / exp_map_background_downsample.data) - solid_angle = 4. * (np.sin(bin_width_x / 2.) * np.sin(bin_width_y / 2.)) * u.steradian data_background = (corrected_counts / solid_angle[np.newaxis, :, :] / energy_axis_computation.bin_width[:, np.newaxis, np.newaxis] /