# Licensed under a 3-clause BSD style license - see LICENSE.rst"""FoV background estimation."""importloggingimportnumpyasnpfromgammapy.mapsimportMap,RegionGeomfromgammapy.modelingimportFitfromgammapy.modeling.modelsimportFoVBackgroundModel,Modelfrom..coreimportMaker__all__=["FoVBackgroundMaker"]log=logging.getLogger(__name__)
[docs]classFoVBackgroundMaker(Maker):"""Normalize template background on the whole field-of-view. The dataset background model can be simply scaled (method="scale") or fitted (method="fit") on the dataset counts. The normalization is performed outside the exclusion mask that is passed on init. This also internally takes into account the dataset fit mask. If a SkyModel is set on the input dataset its parameters are frozen during the FoV re-normalization. If the requirement (greater than) of either min_counts or min_npred_background is not satisfied, the background will not be normalised. Parameters ---------- method : {'scale', 'fit'} The normalization method to be applied. Default 'scale'. exclusion_mask : `~gammapy.maps.WcsNDMap` Exclusion mask. spectral_model : SpectralModel or str Reference norm spectral model to use for the `FoVBackgroundModel`, if none is defined on the dataset. By default, use pl-norm. min_counts : int Minimum number of counts, or residuals counts if a SkyModel is set, required outside the exclusion region. min_npred_background : float Minimum number of predicted background counts required outside the exclusion region. """tag="FoVBackgroundMaker"available_methods=["fit","scale"]def__init__(self,method="scale",exclusion_mask=None,spectral_model="pl-norm",min_counts=0,min_npred_background=0,fit=None,):self.method=methodself.exclusion_mask=exclusion_maskself.min_counts=min_countsself.min_npred_background=min_npred_backgroundifisinstance(spectral_model,str):spectral_model=Model.create(tag=spectral_model,model_type="spectral")ifnotspectral_model.is_norm_spectral_model:raiseValueError("Spectral model must be a norm spectral model")self.default_spectral_model=spectral_modeliffitisNone:fit=Fit()self.fit=fit@propertydefmethod(self):"""Method property."""returnself._method@method.setterdefmethod(self,value):"""Method setter."""ifvaluenotinself.available_methods:raiseValueError(f"Not a valid method for FoVBackgroundMaker: {value}."f" Choose from {self.available_methods}")self._method=value
[docs]defmake_default_fov_background_model(self,dataset):"""Add FoV background model to the model definition. Parameters ---------- dataset : `~gammapy.datasets.MapDataset` Input map dataset. Returns ------- dataset : `~gammapy.datasets.MapDataset` Map dataset including background model. """bkg_model=FoVBackgroundModel(dataset_name=dataset.name,spectral_model=self.default_spectral_model.copy())ifdataset.modelsisNone:dataset.models=bkg_modelelse:dataset.models=dataset.models+bkg_modelreturndataset
def_make_masked_summed_counts(self,dataset):"""Compute the sums of the counts, npred, and background maps within the mask."""npred=dataset.npred()mask=dataset.mask&~np.isnan(npred)count_tot=dataset.counts.data[mask].sum()npred_tot=npred.data[mask].sum()bkg_tot=dataset.npred_background().data[mask].sum()return{"counts":count_tot,"npred":npred_tot,"bkg":bkg_tot,}def_verify_requirements(self,dataset):"""Verify that the requirements of min_counts and min_npred_background are satisfied."""total=self._make_masked_summed_counts(dataset)not_bkg_tot=total["npred"]-total["bkg"]value=(total["counts"]-not_bkg_tot)/total["bkg"]ifnotnp.isfinite(value):log.warning(f"FoVBackgroundMaker failed. Non-finite normalisation value for {dataset.name}. ""Setting mask to False.")returnFalseeliftotal["bkg"]<=self.min_npred_background:log.warning(f"FoVBackgroundMaker failed. Only {int(total['bkg'])} background"f" counts outside exclusion mask for {dataset.name}. ""Setting mask to False.")returnFalseeliftotal["counts"]<=self.min_counts:log.warning(f"FoVBackgroundMaker failed. Only {int(total['counts'])} counts "f"outside exclusion mask for {dataset.name}. ""Setting mask to False.")returnFalseeliftotal["counts"]-not_bkg_tot<=0:log.warning(f"FoVBackgroundMaker failed. Negative residuals counts for {dataset.name}. ""Setting mask to False.")returnFalseelse:returnTrue
[docs]defrun(self,dataset,observation=None):"""Run FoV background maker. Parameters ---------- dataset : `~gammapy.datasets.MapDataset` Input map dataset. """ifisinstance(dataset.counts.geom,RegionGeom):raiseTypeError("FoVBackgroundMaker does not support region based datasets.")mask_fit=dataset.mask_fitifmask_fit:dataset.mask_fit*=self.make_exclusion_mask(dataset)else:dataset.mask_fit=self.make_exclusion_mask(dataset)ifdataset.background_modelisNone:dataset=self.make_default_fov_background_model(dataset)ifself._verify_requirements(dataset)isTrue:ifself.method=="fit":dataset=self.make_background_fit(dataset)else:# always scale the background firstdataset=self.make_background_scale(dataset)else:dataset.mask_safe.data[...]=Falsedataset.mask_fit=mask_fitreturndataset
[docs]defmake_background_fit(self,dataset):"""Fit the FoV background model on the dataset counts data. Parameters ---------- dataset : `~gammapy.datasets.MapDataset` Input dataset. Returns ------- dataset : `~gammapy.datasets.MapDataset` Map dataset with fitted background model. """# freeze all model components not related to background modelmodels=dataset.models.select(tag="sky-model")withmodels.restore_status(restore_values=False):models.select(tag="sky-model").freeze()fit_result=self.fit.run(datasets=[dataset])ifnotfit_result.success:log.warning(f"FoVBackgroundMaker failed. Fit did not converge for {dataset.name}. ""Setting mask to False.")dataset.mask_safe.data[...]=Falsereturndataset
[docs]defmake_background_scale(self,dataset):"""Fit the FoV background model on the dataset counts data. Parameters ---------- dataset : `~gammapy.datasets.MapDataset` Input dataset. Returns ------- dataset : `~gammapy.datasets.MapDataset` Map dataset with scaled background model. """total=self._make_masked_summed_counts(dataset)not_bkg_tot=total["npred"]-total["bkg"]value=(total["counts"]-not_bkg_tot)/total["bkg"]error=np.sqrt(total["counts"]-not_bkg_tot)/total["bkg"]dataset.models[f"{dataset.name}-bkg"].spectral_model.norm.value=valuedataset.models[f"{dataset.name}-bkg"].spectral_model.norm.error=errorreturndataset