From 9ab133147753c2203dc2471b14e25ac132bb913c Mon Sep 17 00:00:00 2001
From: Matthias Plum <mplum@icecube.wisc.edu>
Date: Tue, 31 Oct 2023 11:29:35 -0400
Subject: [PATCH] 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                      | 8 ++++++++
 docs/icetop_tutorial.rst           | 3 +++
 src/simweights/_icetop_weighter.py | 5 +++--
 tests/test_icetop_datasets.py      | 2 +-
 tests/test_icetop_weighter.py      | 2 +-
 tests/test_spatial.py              | 6 +++++-
 6 files changed, 21 insertions(+), 5 deletions(-)

diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 451c3c4..5f10ac6 100644
--- a/CHANGELOG.rst
+++ b/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
 ---------------------
 
diff --git a/docs/icetop_tutorial.rst b/docs/icetop_tutorial.rst
index 8ccbc2f..cc18f8c 100644
--- a/docs/icetop_tutorial.rst
+++ b/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).
diff --git a/src/simweights/_icetop_weighter.py b/src/simweights/_icetop_weighter.py
index 34ae24a..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 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]),
diff --git a/tests/test_icetop_datasets.py b/tests/test_icetop_datasets.py
index 0428f10..4a14f7c 100755
--- a/tests/test_icetop_datasets.py
+++ b/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"]
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)
diff --git a/tests/test_spatial.py b/tests/test_spatial.py
index 0f2770f..98cbcf3 100755
--- a/tests/test_spatial.py
+++ b/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