add option of changable background level

figures_of_merit.py

@@ -20,12 +20,12 @@ DIFF.fc.__doc__ = r"90% upper limit (Feldman-Cousins construction)"
 
 
 class GZK(object):
-    def __init__(self, exposures):
-        self.bundle = factory.component_bundle(exposures, self.make_components)
+    def __init__(self, exposures, **fixed):
+        self.bundle = factory.component_bundle(exposures, self.make_components, **fixed)
 
     def benchmark(self, fom, nobackground=False, **kwargs):
         components = self.bundle.get_components()
-        components["gamma"] = multillh.NuisanceParam(-12.3, 0.5, min=-12.7, max=-11.7)
+        components["gamma"] = multillh.NuisanceParam(-2.3, 0.5, min=-2.7, max=-1.7)
         components["uhe_gamma"] = multillh.NuisanceParam(-2, 0.5, min=-2.7, max=-1.7)
 
         #print(f"... components: {list(components.keys())}")
@@ -91,7 +91,7 @@ class GZK(object):
         return dict(ns=ns, nb=nb)
 
     @staticmethod
-    def make_components(aeffs):
+    def make_components(aeffs, **kwargs):
         aeff, muon_aeff = aeffs
         energy_threshold = np.inf
         atmo = diffuse.AtmosphericNu.conventional(
@@ -108,108 +108,22 @@ class GZK(object):
         uhe = diffuse.DiffuseAstro(aeff, 1.0, gamma_name="uhe_gamma")
         gzk = diffuse.ReasonableGZK(aeff, 1.0) #edited by steffen
         #return dict(atmo=atmo, prompt=prompt, astro=astro, gzk=gzk, uhe=uhe)
-        if muon_aeff is not None:
-            muon = radio_aeff_generation.MuonBackground(
-                  muon_aeff, 1
-              )#diffuse.One(muon_aeff, 1)
-            return dict(gzk=gzk, uhe=uhe, muons=muon)
-        else:
-            return dict(gzk=gzk, uhe=uhe)
-
-class GZK_mu5(object):
-    def __init__(self, exposures):
-        self.bundle = factory.component_bundle(exposures, self.make_components)
-
-    def benchmark(self, fom, nobackground=False, **kwargs):
-        components = self.bundle.get_components()
-        components["gamma"] = multillh.NuisanceParam(-12.3, 0.5, min=-12.7, max=-11.7)
-        components["uhe_gamma"] = multillh.NuisanceParam(-2, 0.5, min=-2.7, max=-1.7)
-
-        print(components.keys())
-        if nobackground:
-            components.pop("muons")
-        gzk = components.pop("gzk")
-        print(gzk)
-        uhe = components.pop("uhe")
-        print(uhe)
-        if "uhe_gamma" not in kwargs:
-            kwargs["uhe_gamma"] = -2.0
-
-        if fom == TOT.ul:
-            return pointsource.upper_limit(
-                uhe, components, baseline=100, tolerance=1e-4, gamma=-2.3, **kwargs
-            )
-        elif fom == TOT.dp:
-            return pointsource.discovery_potential(
-                uhe, components, baseline=100, tolerance=1e-4, gamma=-2.3, **kwargs
-            )
-        elif fom == TOT.fc:
-            return pointsource.fc_upper_limit(uhe, components, gamma=-2.3, **kwargs)
-        elif fom == DIFF.ul:
-            return pointsource.differential_upper_limit(
-                uhe, components, tolerance=1e-4, gamma=-2.3, **kwargs
-            )
-        elif fom == DIFF.dp:
-            return pointsource.differential_discovery_potential(
-                uhe, components, tolerance=1e-4, gamma=-2.3, **kwargs
-            )
-        elif fom == DIFF.fc:
-            return pointsource.differential_fc_upper_limit(
-                uhe, components, gamma=-2.3, **kwargs
-            )
-        else:
-            raise RuntimeError("No such fom")
 
-    def event_numbers(self, sigma=5):
-        """Returns event numbers from ahlers flux needed to reject null
-        hypothesis at sigma *sigma*
-        """
-        scale = self.benchmark(TOT.dp, sigma=sigma)
-
-        components = self.bundle.get_components()
-        components.pop("uhe")
-        components["gamma"] = multillh.NuisanceParam(-2.3, 0.5, min=-2.7, max=-1.7)
-        llh = multillh.asimov_llh(components)
-
-        bin_edges = components["gzk"].bin_edges
-
-        exes = multillh.get_expectations(llh, gzk=scale, gamma=-2.3)
-        nb = pointsource.events_above(
-            exes["atmo"], bin_edges
-        ) + pointsource.events_above(exes["astro"], bin_edges)
-        ns = pointsource.events_above(exes["gzk"], bin_edges)
-        return dict(ns=ns, nb=nb)
-
-    @staticmethod
-    def make_components(aeffs):
-        aeff, muon_aeff = aeffs
-        energy_threshold = np.inf
-        atmo = diffuse.AtmosphericNu.conventional(
-            aeff, 1.0, hard_veto_threshold=energy_threshold
-        )
-        atmo.uncertainty = 0.1
-        prompt = diffuse.AtmosphericNu.prompt(
-            aeff, 1.0, hard_veto_threshold=energy_threshold
-        )
-        prompt.min = 0.5
-        prompt.max = 3
-        astro = diffuse.DiffuseAstro(aeff, 1.0)
-        astro.seed = 2
-        uhe = diffuse.DiffuseAstro(aeff, 1.0, gamma_name="uhe_gamma")
-        gzk = diffuse.ReasonableGZK(aeff, 1.0) #edited by steffen
-        #return dict(atmo=atmo, prompt=prompt, astro=astro, gzk=gzk, uhe=uhe)
+        muon_norm = 1
+        if "muon_norm" in kwargs:
+            muon_norm = kwargs["muon_norm"]
         if muon_aeff is not None:
             muon = radio_aeff_generation.MuonBackground(
-                  muon_aeff, 5
-              )#diffuse.One(muon_aeff, 1)
+                  muon_aeff, muon_norm
+              )
             return dict(gzk=gzk, uhe=uhe, muons=muon)
         else:
             return dict(gzk=gzk, uhe=uhe)
 
 class PointSource(object):
-    def __init__(self, exposures, zi):
+    def __init__(self, exposures, zi, **fixed):
         self.bundle = factory.component_bundle(
-            exposures, partial(self.make_components, zi)
+            exposures, partial(self.make_components, zi), **fixed
         )
 
     def benchmark(self, fom, nobackground=False, gamma=-2.0, diff_gamma=-3.17, **kwargs):
@@ -282,16 +196,16 @@ class PointSource(object):
             raise RuntimeError("No such fom")
 
     def background_optimised_benchmark(self, fom, nobackground=False, **kwargs):
-        energies, result = self.benchmark(fom, nobackground=nobackground, **kwargs)
+        result = self.benchmark(fom, nobackground=nobackground, **kwargs)
         if nobackground==False:
             E_cut_scan = np.logspace(6,9,30)
             for E_cut in E_cut_scan:
-                theE, theRes = self.benchmark(fom, nobackground=nobackground, ecutoff=E_cut, **kwargs)
-                result = np.nanmin(np.array([result, theRes]), axis=0)
-        return (energies, result)
+                theRes = self.benchmark(fom, nobackground=nobackground, ecutoff=E_cut, **kwargs)
+                result = np.nanmin([result, theRes])
+        return result
 
     @staticmethod
-    def make_components(zi, aeffs):
+    def make_components(zi, aeffs, **kwargs):
         aeff, muon_aeff = aeffs
         atmo = diffuse.AtmosphericNu.conventional(aeff, 1.0, veto_threshold=None)
         atmo.uncertainty = 0.1
@@ -308,12 +222,16 @@ class PointSource(object):
         gzk_bkg = gzk.point_source_background(zenith_index=zi)
 
         components = dict(atmo=atmo_bkg, prompt=prompt_bkg, astro=astro_bkg, ps=ps, gzk=gzk_bkg)
+
+        muon_norm = 1
+        if "muon_norm" in kwargs:
+              muon_norm = kwargs["muon_norm"]
         if muon_aeff is not None:
             components["muon"] = radio_aeff_generation.MuonBackground(
-                muon_aeff, 1
+                muon_aeff, muon_norm
             ).point_source_background(zenith_index=zi, psi_bins=aeff.bin_edges[-1][:-1])
         #components.pop("astro")
         components.pop("atmo")
         components.pop("prompt")
-        components.pop("muon")
+        #components.pop("muon")
         return components