From 35ea099f2a8efe1d370b57496ef1ff87df6ce3fd Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Wed, 13 Sep 2023 16:05:52 +0200 Subject: [PATCH 1/2] Fix PowerLaw.from_simulation for new SimulationInfo, add test --- pyirf/io/eventdisplay.py | 3 ++- pyirf/spectral.py | 7 ++++--- pyirf/tests/test_spectral.py | 39 ++++++++++++++++++++++++++++++++++++ 3 files changed, 45 insertions(+), 4 deletions(-) diff --git a/pyirf/io/eventdisplay.py b/pyirf/io/eventdisplay.py index 5c8131923..a46054c4b 100644 --- a/pyirf/io/eventdisplay.py +++ b/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 diff --git a/pyirf/spectral.py b/pyirf/spectral.py index e4db8b8db..35bad8647 100644 --- a/pyirf/spectral.py +++ b/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 diff --git a/pyirf/tests/test_spectral.py b/pyirf/tests/test_spectral.py index 0fab23dec..fcf687120 100644 --- a/pyirf/tests/test_spectral.py +++ b/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)) From ee8ddaf630f85f6840c0c838c217cc51ea779685 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Wed, 13 Sep 2023 16:07:24 +0200 Subject: [PATCH 2/2] Add changelog --- docs/changes/258.bugfix.rst | 2 ++ 1 file changed, 2 insertions(+) create mode 100644 docs/changes/258.bugfix.rst diff --git a/docs/changes/258.bugfix.rst b/docs/changes/258.bugfix.rst new file mode 100644 index 000000000..b4c32a6d4 --- /dev/null +++ b/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``.