Skip to content

Commit

Permalink
Test the volscale line, even if we're only assuming a volume scale of…
Browse files Browse the repository at this point in the history
… 1.0.
  • Loading branch information
mjlarson committed Aug 29, 2024
1 parent 5398f42 commit a14a719
Showing 1 changed file with 41 additions and 34 deletions.
75 changes: 41 additions & 34 deletions tests/test_genie_weighter.py
Original file line number Diff line number Diff line change
Expand Up @@ -33,42 +33,49 @@ def test_genie(self):
c1 = simweights.CircleInjector(300, 0, 1)
p1 = simweights.PowerLaw(0, 1e3, 1e4)

weight = np.zeros(nevents, dtype=result_dtype)
weight["neu"] = pdgid
weight["pzv"] = coszen
weight["Ev"] = p1.ppf(np.linspace(0, 1, nevents))

for event_weight in [1e-6, 1e-3, 1]:
weight["wght"] = event_weight

for nfiles in [1, 5, 50]:
rows = nfiles * [
(
pdgid,
nevents,
1,
c1.radius,
np.arccos(c1.cos_zen_max),
np.arccos(c1.cos_zen_min),
p1.a,
p1.b,
p1.g,
),
]
info = np.array(rows, dtype=info_dtype)
d = {"I3GenieResult": weight, "I3GenieInfo": info}

for flux in [0.1, 1, 10]:
wobj = simweights.GenieWeighter(d)
w = wobj.get_weights(flux)
np.testing.assert_allclose(
w.sum(),
flux * event_weight * c1.etendue * p1.integral / nfiles,
)
E = d["I3GenieResult"]["Ev"]
y, x = np.histogram(E, weights=w, bins=51, range=[p1.a, p1.b])
Ewidth = np.ediff1d(x)
np.testing.assert_allclose(y, flux * event_weight * Ewidth * c1.etendue / nfiles, 5e-3)
for include_volscale in [True, False]:
result_dtype = [("neu", np.int32), ("pzv", np.float64), ("Ev", np.float64), ("wght", np.float64)]
if include_volscale:
result_dtype.append(("volscale", np.float64))

weight = np.zeros(nevents, dtype=result_dtype)
weight["neu"] = pdgid
weight["pzv"] = coszen
weight["Ev"] = p1.ppf(np.linspace(0, 1, nevents))
weight["wght"] = event_weight

if include_volscale:
weight["volscale"] = 1

rows = nfiles * [
(
pdgid,
nevents,
1,
c1.radius,
np.arccos(c1.cos_zen_max),
np.arccos(c1.cos_zen_min),
p1.a,
p1.b,
p1.g,
),
]
info = np.array(rows, dtype=info_dtype)
d = {"I3GenieResult": weight, "I3GenieInfo": info}

for flux in [0.1, 1, 10]:
wobj = simweights.GenieWeighter(d)
w = wobj.get_weights(flux)
np.testing.assert_allclose(
w.sum(),
flux * event_weight * c1.etendue * p1.integral / nfiles,
)
E = d["I3GenieResult"]["Ev"]
y, x = np.histogram(E, weights=w, bins=51, range=[p1.a, p1.b])
Ewidth = np.ediff1d(x)
np.testing.assert_allclose(y, flux * event_weight * Ewidth * c1.etendue / nfiles, 5e-3)

with self.assertRaises(TypeError):
simweights.GenieWeighter(d, nfiles=10)
Expand Down

0 comments on commit a14a719

Please sign in to comment.