Icetop effective area fix 2try (#22)

* Second try of fixing #14 started from a clean branch

* possible fix for the failing tests and removal of unnecessary object

* Update changelog with additional information for IceTop fix

* Update IceTop documentation with effective area information

CHANGELOG.rst

@@ -5,6 +5,14 @@
 Changelog
 =========
 
+main - 2023-10-30
+---------------------
+Update the calculation of IceTop weights to fix an issue with the effective area.
+The fix replaces the CircularInjector with NaturalRateCylinder,
+where the cylinder height is set to zero to mimic the flat surface detector.
+The calculated effective area is the projected area, which can be corrected by
+dividing by cos(zenith).
+
 `0.1.0`_ - 2023-01-24
 ---------------------
 

docs/icetop_tutorial.rst

@@ -58,3 +58,6 @@ combining low- and high-energy samples from each nuclear type (for instance, dat
 protons), so that it spans the full energy range from 10^5 -> 10^9.6 GeV.
 
 .. figure:: icetop_eprimary_combinedLEHE_2012_H4a.svg
+
+Note that the calculated effective area by simweights is the projected effective area,
+which can be corrected by dividing by cos(zenith).

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 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 = CircleInjector(
+        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]),

tests/test_icetop_datasets.py

@@ -31,7 +31,7 @@ def cmp_dataset(self, fname):
         assert len(reffile["I3TopInjectorInfo"]) == 1
         si = reffile["I3TopInjectorInfo"][0]
         pri = reffile["MCPrimary"]
-        solid_angle = 2 * np.pi * (np.cos(si["min_zenith"]) - np.cos(si["max_zenith"]))
+        solid_angle = np.pi * (np.cos(si["min_zenith"]) ** 2 - np.cos(si["max_zenith"]) ** 2)
         injection_area = np.pi * (si["sampling_radius"] * 1e2) ** 2
         energy_integral = np.log(si["max_energy"] / si["min_energy"])  # assuming E^-1
         energy_factor = energy_integral * pri["energy"]

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)

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