From f5ffd0898d576a4f3fffbdba3e388c2b865c4174 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maximilian=20N=C3=B6the?= Date: Tue, 29 Sep 2020 12:36:56 +0200 Subject: [PATCH 1/3] Add method integrating sim. events in FOV bins MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Michele Peresano Co-authored-by: Thomas Villaume Co-authored-by: Léa Jouvin --- pyirf/irf/effective_area.py | 2 +- pyirf/simulations.py | 125 +++++++++++++++++++++++++++++--- pyirf/tests/test_simulations.py | 82 +++++++++++++++++++++ 3 files changed, 199 insertions(+), 10 deletions(-) create mode 100644 pyirf/tests/test_simulations.py diff --git a/pyirf/irf/effective_area.py b/pyirf/irf/effective_area.py index 017240540..33341c0cc 100644 --- a/pyirf/irf/effective_area.py +++ b/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) diff --git a/pyirf/simulations.py b/pyirf/simulations.py index 3850e5b92..46e80ec18 100644 --- a/pyirf/simulations.py +++ b/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 diff --git a/pyirf/tests/test_simulations.py b/pyirf/tests/test_simulations.py new file mode 100644 index 000000000..bedce85e3 --- /dev/null +++ b/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)) From dd68ad07f6cad39a683a2d48aee956b25204a46c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maximilian=20N=C3=B6the?= Date: Tue, 29 Sep 2020 13:17:15 +0200 Subject: [PATCH 2/3] Remove unused import --- pyirf/tests/test_simulations.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pyirf/tests/test_simulations.py b/pyirf/tests/test_simulations.py index bedce85e3..3639670dc 100644 --- a/pyirf/tests/test_simulations.py +++ b/pyirf/tests/test_simulations.py @@ -1,5 +1,4 @@ import numpy as np -from pyirf.io import read_eventdisplay_fits import astropy.units as u From cbf2f0745ef6875a04cf749405aacf6b35b94aea Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maximilian=20N=C3=B6the?= Date: Tue, 29 Sep 2020 13:53:00 +0200 Subject: [PATCH 3/3] Add test for error --- pyirf/tests/test_simulations.py | 57 ++++++++++++++++++++++----------- 1 file changed, 39 insertions(+), 18 deletions(-) diff --git a/pyirf/tests/test_simulations.py b/pyirf/tests/test_simulations.py index 3639670dc..c7b14fa68 100644 --- a/pyirf/tests/test_simulations.py +++ b/pyirf/tests/test_simulations.py @@ -1,11 +1,12 @@ import numpy as np import astropy.units as u +import pytest def test_integrate_energy(): from pyirf.simulations import SimulatedEventsInfo - proton_info = SimulatedEventsInfo( + info = SimulatedEventsInfo( n_showers=int(1e6), energy_min=100 * u.GeV, energy_max=10 * u.TeV, @@ -14,29 +15,29 @@ def test_integrate_energy(): 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) + energy_bins = np.geomspace(info.energy_min, info.energy_max, 20) assert np.isclose( - proton_info.calculate_n_showers_per_energy(energy_bins).sum(), - proton_info.n_showers, + info.calculate_n_showers_per_energy(energy_bins).sum(), + 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) + energy_bins = np.geomspace(info.energy_min, info.energy_max, 20) + energy_bins = np.append(0.9 * info.energy_min, energy_bins) + energy_bins = np.append(energy_bins, 1.1 * info.energy_max) - events = proton_info.calculate_n_showers_per_energy(energy_bins) - assert np.isclose(proton_info.n_showers, events.sum()) + events = info.calculate_n_showers_per_energy(energy_bins) + assert np.isclose(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 + 0.5 * info.energy_min, 2 * info.energy_max, 5 ) - events = proton_info.calculate_n_showers_per_energy(energy_bins) - assert np.isclose(proton_info.n_showers, events.sum()) + events = info.calculate_n_showers_per_energy(energy_bins) + assert np.isclose(info.n_showers, events.sum()) assert events[0] > 0 assert events[-1] > 0 @@ -45,7 +46,7 @@ def test_integrate_energy_fov(): from pyirf.simulations import SimulatedEventsInfo # simple case, max viewcone on bin edge - proton_info = SimulatedEventsInfo( + info = SimulatedEventsInfo( n_showers=int(1e6), energy_min=100 * u.GeV, energy_max=10 * u.TeV, @@ -55,15 +56,15 @@ def test_integrate_energy_fov(): ) fov_bins = [0, 10, 20] * u.deg - energy_bins = np.geomspace(proton_info.energy_min, proton_info.energy_max, 20) + energy_bins = np.geomspace(info.energy_min, info.energy_max, 20) - n_events = proton_info.calculate_n_showers_per_energy_and_fov(energy_bins, fov_bins) + n_events = 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( + info = SimulatedEventsInfo( n_showers=int(1e6), energy_min=100 * u.GeV, energy_max=10 * u.TeV, @@ -73,9 +74,29 @@ def test_integrate_energy_fov(): ) 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) + energy_bins = np.geomspace(info.energy_min, info.energy_max, 20) + n_events = 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)) + + +def test_integrate_energy_fov_pointlike(): + from pyirf.simulations import SimulatedEventsInfo + + 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=0 * u.deg, + ) + + fov_bins = [0, 9, 11, 20] * u.deg + energy_bins = np.geomspace(info.energy_min, info.energy_max, 20) + + # make sure we raise an error on invalid input + with pytest.raises(ValueError): + info.calculate_n_showers_per_energy_and_fov(energy_bins, fov_bins)