progress on #14 Use NaturalRateCylinder for IceTop

icecube · Oct 23, 2023 · 6d583de · 6d583de
1 parent 447561e
commit 6d583de
Show file tree

Hide file tree

Showing 3 changed files with 5 additions and 5 deletions.
diff --git 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
@@ -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
@@ -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)