diff --git a/figures_of_merit.py b/figures_of_merit.py
index 19f56c0ae4f3a157d383a90cc285c1154d2d74fb..bc44483195e6eb27b30af097681f150017daf2e4 100644
--- a/figures_of_merit.py
+++ b/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