spectral asimov

titrate/datasets.py

@@ -1,5 +1,5 @@
 import numpy as np
-from gammapy.datasets import MapDataset
+from gammapy.datasets import MapDataset, SpectrumDataset
 
 from titrate.utils import copy_models_to_dataset
 
@@ -40,3 +40,40 @@ def from_MapDataset(self, dataset):
         asimov_dataset.fake()
 
         return asimov_dataset
+
+
+class AsimovSpectralDataset(SpectrumDataset):
+    """
+    AsimovMapDataset is a subclass of a gammapy MapDataset
+    and provides asimov-like fake method.
+    """
+
+    def fake(self):
+        """
+        This method generates Asimov like counts,
+        i.e. the counts are not drawn from a poisson distribution.
+        """
+        npred = self.npred()
+        data = np.nan_to_num(npred.data, copy=True, nan=0.0, posinf=0.0, neginf=0.0)
+        npred.data = data
+        self.counts = npred
+
+    @classmethod
+    def from_SpectralDataset(self, dataset):
+        dataset_dict = dataset.__dict__.copy()
+
+        delete_keys = [key for key in dataset_dict.keys() if key.startswith("_")]
+
+        deleted_entries = {}
+        for key in delete_keys:
+            deleted_entries[key] = dataset_dict.pop(key)
+
+        asimov_dataset = AsimovSpectralDataset(**dataset_dict)
+        # copy IRFs separately
+        asimov_dataset._background_parameters_cached = deleted_entries[
+            "_background_parameters_cached"
+        ]
+        copy_models_to_dataset(dataset.models, asimov_dataset)
+        asimov_dataset.fake()
+
+        return asimov_dataset

titrate/statistics.py

@@ -58,7 +58,7 @@ def __init__(self, dataset, poi_name=""):
             )
 
         self.fit = Fit()
-        self.fit_result = self.fit.run(datasets=[self.dataset])
+        self.fit_result = self.fit.run(datasets=self.dataset)
         self.poi_best = self.fit_result.parameters[self.poi_name].value
         self.likelihood_minimum = self.dataset.stat_sum()
 

titrate/upperlimits.py

@@ -4,6 +4,7 @@
 import numpy as np
 from astropy.table import QTable
 from gammapy.astro.darkmatter import DarkMatterAnnihilationSpectralModel
+from gammapy.datasets import Datasets
 from gammapy.modeling import Fit
 from gammapy.modeling.models import (
     FoVBackgroundModel,
@@ -15,7 +16,7 @@
 from scipy.optimize import brentq
 from scipy.stats import norm
 
-from titrate.datasets import AsimovMapDataset
+from titrate.datasets import AsimovMapDataset, AsimovSpectralDataset
 from titrate.statistics import QMuTestStatistic, QTildeMuTestStatistic
 from titrate.utils import copy_models_to_dataset
 
@@ -67,6 +68,12 @@ def compute(self):
             ]
         ).ravel()
         interpolation = interp1d(interp_ul_points, interp_pvalues - 1 + self.cl)
+        if interpolation(interp_ul_points[0]) < 0:
+            print(
+                interpolation(interp_ul_points[0]),
+                interpolation(interp_ul_points[-1]),
+                self.measurement_dataset,
+            )
         poi_ul = brentq(interpolation, poi_ul / 2, poi_ul)
 
         return poi_ul
@@ -113,12 +120,18 @@ def pvalue(self, poi_ul, cl_type):
 
     def expected_uls(self):
         # Create asimov dataset
-        asimov_dataset = AsimovMapDataset.from_MapDataset(self.measurement_dataset)
+        # asimov_dataset = AsimovMapDataset.from_MapDataset(self.measurement_dataset)
+        asimov_dataset = []
+        for dataset in self.measurement_dataset:
+            asimov_dataset.append(AsimovSpectralDataset.from_SpectralDataset(dataset))
+        asimov_dataset = Datasets(asimov_dataset)
+        # asimov_dataset =
+        # AsimovSpectralDataset.from_SpectralDatasets(self.measurement_dataset)
         asimov_dataset.models.parameters[self.poi_name].value = 0
-        asimov_dataset.fake()
+        # asimov_dataset.fake()
 
         fit = Fit()
-        fit_result = fit.run(datasets=[asimov_dataset])
+        fit_result = fit.run(datasets=asimov_dataset)
 
         sigma = np.sqrt(fit_result.covariance_result.matrix[0, 0])
         med_ul = self.compute_band(sigma, 0, self.cl_type)
@@ -150,7 +163,8 @@ def __init__(
         mass_min,
         mass_max,
         n_steps,
-        jfactor_map,
+        jfactor,
+        analysis="3d",
         cl_type="s",
         cl=0.95,
         max_workers=None,
@@ -161,7 +175,8 @@ def __init__(
         self.masses = np.geomspace(
             mass_min.to_value("TeV"), mass_max.to_value("TeV"), n_steps
         )
-        self.jfactor_map = jfactor_map
+        self.jfactor = jfactor
+        self.analysis = analysis
         self.cl_type = cl_type
         self.cl = cl
         self.kwargs = kwargs
@@ -173,19 +188,34 @@ def __init__(
 
     def setup_models(self):
         models = []
-        spatial_model = TemplateSpatialModel(self.jfactor_map, normalize=False)
-        bkg_model = FoVBackgroundModel(dataset_name="foo")
-        for channel in self.channels:
-            for mass in self.masses:
-                spectral_model = DarkMatterAnnihilationSpectralModel(
-                    mass=mass * u.TeV, channel=channel
-                )
-                sky_model = SkyModel(
-                    spatial_model=spatial_model,
-                    spectral_model=spectral_model,
-                    name=f"mass_{mass:.2f}TeV_channel_{channel}",
-                )
-                models.append(Models([sky_model, bkg_model]))
+        if self.analysis == "3d":
+            spatial_model = TemplateSpatialModel(self.jfactor, normalize=False)
+
+            bkg_model = FoVBackgroundModel(dataset_name="foo")
+            for channel in self.channels:
+                for mass in self.masses:
+                    spectral_model = DarkMatterAnnihilationSpectralModel(
+                        mass=mass * u.TeV, channel=channel
+                    )
+                    sky_model = SkyModel(
+                        spatial_model=spatial_model,
+                        spectral_model=spectral_model,
+                        name=f"mass_{mass:.2f}TeV_channel_{channel}",
+                    )
+                    models.append(Models([sky_model, bkg_model]))
+
+        elif self.analysis == "1d":
+            bkg_model = FoVBackgroundModel(dataset_name="foo")
+            for channel in self.channels:
+                for mass in self.masses:
+                    spectral_model = DarkMatterAnnihilationSpectralModel(
+                        mass=mass * u.TeV, channel=channel, jfactor=self.jfactor
+                    )
+                    sky_model = SkyModel(
+                        spectral_model=spectral_model,
+                        name=f"mass_{mass:.2f}TeV_channel_{channel}",
+                    )
+                    models.append(Models([sky_model, bkg_model]))
 
         return models
 

titrate/utils.py

@@ -1,3 +1,6 @@
+from gammapy.datasets import Datasets
+
+
 def CopyModelError(AttributeError):
     """Something went wrong during copying a model."""
 
@@ -27,15 +30,20 @@ def copy_dataset_with_models(dataset):
 def copy_models_to_dataset(models, dataset):
     """Copies models and assigns them to dataset."""
     model_copies = models.copy()
-    for model in model_copies:
-        if hasattr(model, "_name"):
-            model._name = f"{dataset.name}-{model._name}"
-        elif hasattr(model, "datasets_names"):
-            model.datasets_names = [f"{dataset.name}"]
-        else:
-            raise CopyModelError(
-                f"{model.__class__.__name__} doesn't provided `._name`"
-                f"nor `.datasets_names`."
-            )
-
-    dataset.models = model_copies
+    if isinstance(dataset, Datasets):
+        for d in dataset:
+            b = model_copies[1].copy()
+            b.datasets_names = [d.name]
+            d.models = [model_copies[0], b]
+    else:
+        for model in model_copies:
+            if hasattr(model, "_name"):
+                model._name = f"{dataset.name}-{model._name}"
+            elif hasattr(model, "datasets_names"):
+                model.datasets_names = [dataset.name]
+            else:
+                raise CopyModelError(
+                    f"{model.__class__.__name__} doesn't provided `._name`"
+                    f"nor `.datasets_names`."
+                )
+        dataset.models = model_copies