Add method integrating sim. events in FOV bins

pyirf/irf/effective_area.py

@@ -45,6 +45,6 @@ def point_like_effective_area(selected_events, simulation_info, true_energy_bins
     hist_selected = create_histogram_table(
         selected_events, true_energy_bins, "true_energy"
     )
-    hist_simulated = simulation_info.calculate_n_showers(true_energy_bins)
+    hist_simulated = simulation_info.calculate_n_showers_per_energy(true_energy_bins)
 
     return effective_area(hist_selected["n"], hist_simulated, area)

pyirf/simulations.py

@@ -1,4 +1,5 @@
 import astropy.units as u
+import numpy as np
 
 
 class SimulatedEventsInfo:
@@ -48,9 +49,12 @@ def __init__(
             raise ValueError("spectral index must be <= -1")
 
     @u.quantity_input(energy_bins=u.TeV)
-    def calculate_n_showers(self, energy_bins):
+    def calculate_n_showers_per_energy(self, energy_bins):
         """
-        Calculate number of showers that were simulated in the given interval
+        Calculate number of showers that were simulated in the given energy intervals
+
+        This assumes the events were generated and from a powerlaw
+        like CORSIKA simulates events.
 
         Parameters
         ----------
@@ -65,18 +69,97 @@ def calculate_n_showers(self, energy_bins):
             The actual numbers will follow a poissionian distribution around this
             expected value.
         """
-        bins = energy_bins.to_value(u.TeV)
+        bins = energy_bins
         e_low = bins[:-1]
         e_high = bins[1:]
+        e_min = self.energy_min
+        e_max = self.energy_max
+
+        integral = _powerlaw_pdf_integral(
+            self.spectral_index, e_low, e_high, e_min, e_max
+        )
+
+        integral[e_high <= e_min] = 0
+        integral[e_low >= e_max] = 0
+
+        mask = (e_high > e_max) & (e_low < e_max)
+        integral[mask] = _powerlaw_pdf_integral(
+            self.spectral_index, e_low[mask], e_max, e_min, e_max
+        )
+
+        mask = (e_high > e_min) & (e_low < e_min)
+        integral[mask] = _powerlaw_pdf_integral(
+            self.spectral_index, e_min, e_high[mask], e_min, e_max
+        )
+
+        return self.n_showers * integral
+
+    def calculate_n_showers_per_fov(self, fov_bins):
+        """
+        Calculate number of showers that were simulated in the given fov bins.
+
+        This assumes the events were generated uniformly distributed per solid angle,
+        like CORSIKA simulates events with the VIEWCONE option.
+
+        Parameters
+        ----------
+        fov_bins: astropy.units.Quantity[angle]
+            The FOV bin edges for which to calculate the number of simulated showers
+
+        Returns
+        -------
+        n_showers: numpy.ndarray(ndim=2)
+            The expected number of events inside each of the ``fov_bins``.
+            This is a floating point number.
+            The actual numbers will follow a poissionian distribution around this
+            expected value.
+        """
+        fov_bins = fov_bins
+        fov_low = fov_bins[:-1]
+        fov_high = fov_bins[1:]
+
+        fov_integral = _viewcone_pdf_integral(self.viewcone, fov_low, fov_high)
+        viewcone = self.viewcone
+        # check if any of the bins are outside the max viewcone
+        fov_integral = np.where(fov_high <= viewcone, fov_integral, 0)
+
+        # identify the bin with the maximum viewcone inside
+        mask = (viewcone > fov_low) & (viewcone < fov_high)
+        fov_integral[mask] = _viewcone_pdf_integral(viewcone, fov_low[mask], viewcone)
 
-        int_index = self.spectral_index + 1
-        e_min = self.energy_min.to_value(u.TeV)
-        e_max = self.energy_max.to_value(u.TeV)
+        return self.n_showers * fov_integral
 
-        e_term = e_low ** int_index - e_high ** int_index
-        normalization = int_index / (e_max ** int_index - e_min ** int_index)
+    @u.quantity_input(energy_bins=u.TeV, fov_bins=u.deg)
+    def calculate_n_showers_per_energy_and_fov(self, energy_bins, fov_bins):
+        """
+        Calculate number of showers that were simulated in the given
+        energy and fov bins.
 
-        return self.n_showers * normalization * e_term
+        This assumes the events were generated uniformly distributed per solid angle,
+        and from a powerlaw in energy like CORSIKA simulates events.
+
+        Parameters
+        ----------
+        energy_bins: astropy.units.Quantity[energy]
+            The energy bin edges for which to calculate the number of simulated showers
+        fov_bins: astropy.units.Quantity[angle]
+            The FOV bin edges for which to calculate the number of simulated showers
+
+        Returns
+        -------
+        n_showers: numpy.ndarray(ndim=2)
+            The expected number of events inside each of the
+            ``energy_bins`` and ``fov_bins``.
+            Dimension (n_energy_bins, n_fov_bins)
+            This is a floating point number.
+            The actual numbers will follow a poissionian distribution around this
+            expected value.
+        """
+        # energy distribution and fov distribution are independent in CORSIKA,
+        # so just multiply both distributions.
+        e_integral = self.calculate_n_showers_per_energy(energy_bins)
+        fov_integral = self.calculate_n_showers_per_fov(fov_bins)
+        return e_integral[:, np.newaxis] * fov_integral / self.n_showers
 
     def __repr__(self):
         return (
@@ -89,3 +172,27 @@ def __repr__(self):
             f"viewcone={self.viewcone}"
             ")"
         )
+
+
+def _powerlaw_pdf_integral(index, e_low, e_high, e_min, e_max):
+    # strip units, make sure all in the same unit
+    e_low = e_low.to_value(u.TeV)
+    e_high = e_high.to_value(u.TeV)
+    e_min = e_min.to_value(u.TeV)
+    e_max = e_max.to_value(u.TeV)
+
+    int_index = index + 1
+    e_term = e_low ** int_index - e_high ** int_index
+    normalization = int_index / (e_max ** int_index - e_min ** int_index)
+    return e_term * normalization
+
+
+def _viewcone_pdf_integral(viewcone, fov_low, fov_high):
+    if viewcone.value == 0:
+        raise ValueError("Only supported for diffuse simulations")
+    else:
+        norm = 1 / (1 - np.cos(viewcone))
+
+    integral = np.cos(fov_low) - np.cos(fov_high)
+
+    return norm * integral

pyirf/tests/test_simulations.py

@@ -0,0 +1,82 @@
+import numpy as np
+from pyirf.io import read_eventdisplay_fits
+import astropy.units as u
+
+
+def test_integrate_energy():
+    from pyirf.simulations import SimulatedEventsInfo
+
+    proton_info = SimulatedEventsInfo(
+        n_showers=int(1e6),
+        energy_min=100 * u.GeV,
+        energy_max=10 * u.TeV,
+        max_impact=500 * u.m,
+        spectral_index=-2,
+        viewcone=10 * u.deg,
+    )
+    # simplest case, no bins  outside e_min, e_max
+    energy_bins = np.geomspace(proton_info.energy_min, proton_info.energy_max, 20)
+    assert np.isclose(
+        proton_info.calculate_n_showers_per_energy(energy_bins).sum(),
+        proton_info.n_showers,
+    )
+
+    # simple case, e_min and e_max on bin edges
+    energy_bins = np.geomspace(proton_info.energy_min, proton_info.energy_max, 20)
+    energy_bins = np.append(0.9 * proton_info.energy_min, energy_bins)
+    energy_bins = np.append(energy_bins, 1.1 * proton_info.energy_max)
+
+    events = proton_info.calculate_n_showers_per_energy(energy_bins)
+    assert np.isclose(proton_info.n_showers, events.sum())
+
+    assert events[0] == 0
+    assert events[-1] == 0
+
+    # complex case, e_min and e_max inside bins
+    energy_bins = np.geomspace(
+        0.5 * proton_info.energy_min, 2 * proton_info.energy_max, 5
+    )
+    events = proton_info.calculate_n_showers_per_energy(energy_bins)
+    assert np.isclose(proton_info.n_showers, events.sum())
+    assert events[0] > 0
+    assert events[-1] > 0
+
+
+def test_integrate_energy_fov():
+    from pyirf.simulations import SimulatedEventsInfo
+
+    # simple case, max viewcone on bin edge
+    proton_info = SimulatedEventsInfo(
+        n_showers=int(1e6),
+        energy_min=100 * u.GeV,
+        energy_max=10 * u.TeV,
+        max_impact=500 * u.m,
+        spectral_index=-2,
+        viewcone=10 * u.deg,
+    )
+
+    fov_bins = [0, 10, 20] * u.deg
+    energy_bins = np.geomspace(proton_info.energy_min, proton_info.energy_max, 20)
+
+    n_events = proton_info.calculate_n_showers_per_energy_and_fov(energy_bins, fov_bins)
+
+    assert np.all(n_events[:, 1:] == 0)
+    assert np.isclose(np.sum(n_events), int(1e6))
+
+    # viewcone inside of bin
+    proton_info = SimulatedEventsInfo(
+        n_showers=int(1e6),
+        energy_min=100 * u.GeV,
+        energy_max=10 * u.TeV,
+        max_impact=500 * u.m,
+        spectral_index=-2,
+        viewcone=10 * u.deg,
+    )
+
+    fov_bins = [0, 9, 11, 20] * u.deg
+    energy_bins = np.geomspace(proton_info.energy_min, proton_info.energy_max, 20)
+    n_events = proton_info.calculate_n_showers_per_energy_and_fov(energy_bins, fov_bins)
+
+    assert np.all(n_events[:, 1:2] > 0)
+    assert np.all(n_events[:, 2:] == 0)
+    assert np.isclose(np.sum(n_events), int(1e6))