From 6d583deeb297c755f4e933866de54e4b13c55693 Mon Sep 17 00:00:00 2001 From: Matthias Plum Date: Mon, 23 Oct 2023 09:44:09 -0400 Subject: [PATCH] progress on #14 Use NaturalRateCylinder for IceTop --- src/simweights/_icetop_weighter.py | 5 +++-- src/simweights/_spatial.py | 3 +-- tests/test_icetop_weighter.py | 2 +- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/simweights/_icetop_weighter.py b/src/simweights/_icetop_weighter.py index 6ec197b..730fa92 100644 --- a/src/simweights/_icetop_weighter.py +++ b/src/simweights/_icetop_weighter.py @@ -9,7 +9,7 @@ from ._generation_surface import GenerationSurface, generation_surface from ._powerlaw import PowerLaw -from ._spatial import SineSquaredThetaCircleInjector +from ._spatial import NaturalRateCylinder from ._utils import get_column, get_table from ._weighter import Weighter @@ -25,7 +25,8 @@ def sframe_icetop_surface(table: Any) -> GenerationSurface: get_column(table, "min_energy")[i], get_column(table, "max_energy")[i], ) - spatial = SineSquaredThetaCircleInjector( + spatial = NaturalRateCylinder( + 0, # set cylinder height to 0 to get simple surface plane get_column(table, "sampling_radius")[i], np.cos(get_column(table, "max_zenith")[i]), np.cos(get_column(table, "min_zenith")[i]), diff --git a/src/simweights/_spatial.py b/src/simweights/_spatial.py index 438a661..1bfb04f 100644 --- a/src/simweights/_spatial.py +++ b/src/simweights/_spatial.py @@ -184,7 +184,6 @@ def __init__(self, radius: float, cos_zen_min: float, cos_zen_max: float) -> Non self._cap = 1e4 * np.pi * self.radius**2 self.etendue = 2 * np.pi * (self.cos_zen_max**2 - self.cos_zen_min**2) * self._cap self._normalization = 1 / self.etendue - self._normalization_pdf = 1 / (self.cos_zen_max**2 - self.cos_zen_min**2) def projected_area(self, cos_zen: float) -> float: # noqa: ARG002 """Returns the cross sectional area of the injection area in cm^2.""" @@ -193,7 +192,7 @@ def projected_area(self, cos_zen: float) -> float: # noqa: ARG002 def _pdf(self, cosz: NDArray[np.float64]) -> NDArray[np.float64]: sin2_dist = 2 * cosz - return sin2_dist * self._normalization * self._normalization_pdf + return sin2_dist * self._normalization def pdf(self, cos_zen: ArrayLike) -> NDArray[np.float64]: """Returns diff --git a/tests/test_icetop_weighter.py b/tests/test_icetop_weighter.py index 9873125..4b610bb 100755 --- a/tests/test_icetop_weighter.py +++ b/tests/test_icetop_weighter.py @@ -27,7 +27,7 @@ class TestIceTopWeighter(unittest.TestCase): def test_icetop_corsika(self): nevents = 10000 pdgid = 12 - c1 = simweights.CircleInjector(300, 0, 1) + c1 = simweights.NaturalRateCylinder(0, 300, 0, 1) p1 = simweights.PowerLaw(0, 1e3, 1e4) weight = np.zeros(nevents, dtype=particle_dtype)