Refactor scripts/analaysis.py
@@ -4,6 +4,7 @@ import copy
import logging
from collections import defaultdict
from pathlib import Path
+import numpy as np
from astropy import units as u
from astropy.coordinates import Angle, SkyCoord
from regions import CircleSkyRegion
@@ -18,8 +19,9 @@ from gammapy.modeling.serialize import dict_to_models
from gammapy.spectrum import (
FluxPointsDataset,
FluxPointsEstimator,
- ReflectedRegionsBackgroundEstimator,
- SpectrumExtraction,
+ SpectrumDatasetMaker,
+ SafeMaskMaker,
+ ReflectedRegionsBackgroundMaker
)
from gammapy.utils.scripts import make_path, read_yaml
@@ -303,40 +305,51 @@ class Analysis:
on_lat = Angle(region["center"][1])
on_center = SkyCoord(on_lon, on_lat, frame=region["frame"])
on_region = CircleSkyRegion(on_center, Angle(region["radius"]))
- background_params = {"on_region": on_region}
+
+ maker_config = {}
+ if "containment_correction" in self.settings["datasets"]:
+ maker_config["containment_correction"] = self.settings["datasets"][
+ "containment_correction"
+ ]
+ params = self.settings["datasets"]["geom"]["axes"][0]
+ e_reco = MapAxis.from_bounds(**params).edges
+ maker_config["e_reco"] = e_reco
+
+ # TODO: remove hard-coded e_true and make it configurable
+ maker_config["e_true"] = np.logspace(-2, 2.5, 109) * u.TeV
+ maker_config["region"] = on_region
+
+ dataset_maker = SpectrumDatasetMaker(**maker_config)
+ bkg_maker_config = {"region": on_region}
background = self.settings["datasets"]["background"]
+
if "exclusion_mask" in background:
map_hdu = {}
filename = background["exclusion_mask"]["filename"]
if "hdu" in background["exclusion_mask"]:
map_hdu = {"hdu": background["exclusion_mask"]["hdu"]}
exclusion_region = Map.read(filename, **map_hdu)
- background_params["exclusion_mask"] = exclusion_region
+ bkg_maker_config["exclusion_mask"] = exclusion_region
if background["background_estimator"] == "reflected":
- self.background_estimator = ReflectedRegionsBackgroundEstimator(
- observations=self.observations, **background_params
+ reflected_bkg_maker = ReflectedRegionsBackgroundMaker(
+ **bkg_maker_config
)
- self.background_estimator.run()
else:
# TODO: raise error?
log.info("Background estimation only for reflected regions method.")
- extraction_params = {}
- if "containment_correction" in self.settings["datasets"]:
- extraction_params["containment_correction"] = self.settings["datasets"][
- "containment_correction"
- ]
- params = self.settings["datasets"]["geom"]["axes"][0]
- e_reco = MapAxis.from_bounds(**params).edges
- extraction_params["e_reco"] = e_reco
- extraction_params["e_true"] = None
- self.extraction = SpectrumExtraction(
- observations=self.observations,
- bkg_estimate=self.background_estimator.result,
- **extraction_params,
- )
- self.extraction.run()
- self.datasets = Datasets(self.extraction.spectrum_observations)
+ safe_mask_maker = SafeMaskMaker(methods=["aeff-default", "aeff-max"])
+
+ datasets = []
+ for obs in self.observations:
+ selection = ["counts", "aeff", "edisp"]
+ dataset = dataset_maker.run(obs, selection=selection)
+ dataset = reflected_bkg_maker.run(dataset, obs)
+ dataset = safe_mask_maker.run(dataset, obs)
+ datasets.append(dataset)
+
+ self.datasets = Datasets(datasets)
+
if self.settings["datasets"]["stack-datasets"]:
stacked = self.datasets.stack_reduce()
stacked.name = "stacked"
