Second try of fixing #14 started from a clean branch

src/simweights/__init__.py

@@ -18,6 +18,7 @@
     "CircleInjector",
     "CorsikaWeighter",
     "NaturalRateCylinder",
+    "IceTopNaturalRateCylinder",
     "UniformSolidAngleCylinder",
     "TIG1996",
     "FixedFractionFlux",
@@ -63,6 +64,11 @@
 from ._nugen_weighter import NuGenWeighter
 from ._pdgcode import PDGCode
 from ._powerlaw import PowerLaw
-from ._spatial import CircleInjector, NaturalRateCylinder, UniformSolidAngleCylinder
+from ._spatial import (
+    CircleInjector,
+    IceTopNaturalRateCylinder,
+    NaturalRateCylinder,
+    UniformSolidAngleCylinder,
+)
 from ._utils import corsika_to_pdg
 from ._weighter import Weighter

src/simweights/_icetop_weighter.py

@@ -9,7 +9,7 @@
 
 from ._generation_surface import GenerationSurface, generation_surface
 from ._powerlaw import PowerLaw
-from ._spatial import CircleInjector
+from ._spatial import IceTopNaturalRateCylinder
 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 = CircleInjector(
+        spatial = IceTopNaturalRateCylinder(
+            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]),

src/simweights/_spatial.py

@@ -126,6 +126,29 @@ def pdf(self, cos_zen: ArrayLike) -> NDArray[np.float64]:
         )
 
 
+class IceTopNaturalRateCylinder(CylinderBase):
+    r"""Angular distribution when zenith distribution matched the natural rate of an isotropic source.
+
+    For a given zenith angle the intensity of particles thrown was proportional to the cross-sectional area
+    perpendicular to the direction of the particle.
+    """
+
+    def __init__(self, length: float, radius: float, cos_zen_min: float, cos_zen_max: float) -> None:
+        super().__init__(length, radius, cos_zen_min, cos_zen_max)
+        self.etendue = (
+            2 * np.pi * (self.cos_zen_max - self.cos_zen_min) * self._cap
+        )  # This is needed to make the tests work
+        self._normalization = 1 / self.etendue
+
+    def pdf(self, cos_zen: ArrayLike) -> NDArray[np.float64]:
+        cosz = np.asfarray(cos_zen)
+        return np.piecewise(
+            cosz,
+            [(cosz >= self.cos_zen_min) & (cosz <= self.cos_zen_max)],
+            [self._normalization],
+        )
+
+
 class CircleInjector:
     """Particles are generated on a circle perpendicular to the direction of the particle.
 

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.IceTopNaturalRateCylinder(0, 300, 0, 1)
         p1 = simweights.PowerLaw(0, 1e3, 1e4)
 
         weight = np.zeros(nevents, dtype=particle_dtype)

tests/test_spatial.py

@@ -7,7 +7,11 @@
 import unittest
 
 import numpy as np
-from simweights import CircleInjector, NaturalRateCylinder, UniformSolidAngleCylinder
+from simweights import (
+    CircleInjector,
+    NaturalRateCylinder,
+    UniformSolidAngleCylinder,
+)
 from simweights._spatial import CylinderBase