Merge pull request #258 from cta-observatory/fix_from_simulation

Fix PowerLaw.from_simulation for new SimulationInfo, add test

docs/changes/258.bugfix.rst

@@ -0,0 +1,2 @@
+Fix ``PowerLaw.from_simulation`` for the new format of ``SimulatedEventsInformation``,
+it was broken since splitting the single ``viewcone`` into ``viewcone_min`` and ``viewcone_max``.

pyirf/io/eventdisplay.py

@@ -77,7 +77,8 @@ def read_eventdisplay_fits(infile, use_histogram=True):
         energy_max=u.Quantity(run_header["E_range"][1], u.TeV),
         max_impact=u.Quantity(run_header["core_range"][1], u.m),
         spectral_index=run_header["spectral_index"],
-        viewcone=u.Quantity(run_header["viewcone"][1], u.deg),
+        viewcone_min=u.Quantity(run_header["viewcone"][0], u.deg),
+        viewcone_max=u.Quantity(run_header["viewcone"][1], u.deg),
     )
 
     return events, sim_info

pyirf/spectral.py

@@ -116,10 +116,11 @@ def from_simulation(cls, simulated_event_info, obstime, e_ref=1 * u.TeV):
         e_max = simulated_event_info.energy_max
         index = simulated_event_info.spectral_index
         n_showers = simulated_event_info.n_showers
-        viewcone = simulated_event_info.viewcone
+        viewcone_min = simulated_event_info.viewcone_min
+        viewcone_max = simulated_event_info.viewcone_max
 
-        if viewcone.value > 0:
-            solid_angle = 2 * np.pi * (1 - np.cos(viewcone)) * u.sr
+        if (viewcone_max - viewcone_min).value > 0:
+            solid_angle = cone_solid_angle(viewcone_max) - cone_solid_angle(viewcone_min)
             unit = DIFFUSE_FLUX_UNIT
         else:
             solid_angle = 1

pyirf/tests/test_spectral.py

@@ -73,3 +73,42 @@ def test_powerlaw():
     power_law = PowerLaw(1e-10 * unit, -2.65)
     assert power_law(1 * u.TeV).unit == unit
     assert power_law(1 * u.GeV).unit == unit
+
+
+def test_powerlaw_from_simulations():
+    from pyirf.simulations import SimulatedEventsInfo
+    from pyirf.spectral import PowerLaw
+
+    # calculate sensitivity between 1 and 2 degrees offset from fov center
+    obstime = 50 * u.hour
+
+    simulated_events = SimulatedEventsInfo(
+        n_showers=int(1e6),
+        energy_min=10 * u.GeV,
+        energy_max=100 * u.TeV,
+        max_impact=1 * u.km,
+        spectral_index=-2,
+        viewcone_min=0 * u.deg,
+        viewcone_max=0 * u.deg,
+    )
+
+    powerlaw = PowerLaw.from_simulation(simulated_events, obstime=obstime)
+    assert powerlaw.index == -2
+    # regression test, maybe better come up with an easy to analytically verify parameter combination?
+    assert u.isclose(powerlaw.normalization, 1.76856511e-08 / (u.TeV * u.m**2 * u.s))
+
+
+    simulated_events = SimulatedEventsInfo(
+        n_showers=int(1e6),
+        energy_min=10 * u.GeV,
+        energy_max=100 * u.TeV,
+        max_impact=1 * u.km,
+        spectral_index=-2,
+        viewcone_min=5 * u.deg,
+        viewcone_max=10 * u.deg,
+    )
+
+    powerlaw = PowerLaw.from_simulation(simulated_events, obstime=obstime)
+    assert powerlaw.index == -2
+    # regression test, maybe better come up with an easy to analytically verify parameter combination?
+    assert u.isclose(powerlaw.normalization, 2.471917427911683e-07 / (u.TeV * u.m**2 * u.s * u.sr))