From 8511eb3b5e88f4598f7d02ff15ff0cab5dba4473 Mon Sep 17 00:00:00 2001 From: rlopezcoto Date: Thu, 2 Apr 2020 12:45:52 +0200 Subject: [PATCH 1/8] Modify excess_matching_significance from gammapy --- lstchain/mc/sensitivity.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/lstchain/mc/sensitivity.py b/lstchain/mc/sensitivity.py index 14fa40b81c..29b8a89897 100644 --- a/lstchain/mc/sensitivity.py +++ b/lstchain/mc/sensitivity.py @@ -5,7 +5,7 @@ from .mc import rate, weight from lstchain.spectra.crab import crab_hegra,crab_magic from lstchain.spectra.proton import proton_bess -from gammapy.stats.poisson import excess_matching_significance_on_off +from gammapy.stats import WStatCountsStatistic from lstchain.reco.utils import reco_source_position_sky from astropy.coordinates.angle_utilities import angular_separation from lstchain.io import read_simu_info_merged_hdf5 @@ -176,9 +176,13 @@ def calculate_sensitivity_lima(n_excesses, n_background, alpha, n_bins_energy, n """ - n_excesses_5sigma = excess_matching_significance_on_off( \ - n_off=n_background, alpha=alpha, significance=5, method='lima') + stat = WStatCountsStatistic( + ) + n_on=np.ones_like(n_background), + n_off=n_background, + alpha=alpha) + n_excesses_5sigma = stat.excess_matching_significance(5) for i in range(0, n_bins_energy): for j in range(0, n_bins_gammaness): From 72fff95dd736457d708028fe2320e8279304e6f9 Mon Sep 17 00:00:00 2001 From: rlopezcoto Date: Thu, 2 Apr 2020 13:11:34 +0200 Subject: [PATCH 2/8] correct also ebin function --- lstchain/mc/sensitivity.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/lstchain/mc/sensitivity.py b/lstchain/mc/sensitivity.py index 29b8a89897..3fa047ea1d 100644 --- a/lstchain/mc/sensitivity.py +++ b/lstchain/mc/sensitivity.py @@ -177,10 +177,9 @@ def calculate_sensitivity_lima(n_excesses, n_background, alpha, n_bins_energy, n """ stat = WStatCountsStatistic( - ) - n_on=np.ones_like(n_background), - n_off=n_background, - alpha=alpha) + n_on=np.ones_like(n_background), + n_off=n_background, + alpha=alpha) n_excesses_5sigma = stat.excess_matching_significance(5) @@ -219,8 +218,12 @@ def calculate_sensitivity_lima_ebin(n_excesses, n_background, alpha, n_bins_ener a 5 sigma significance """ - n_excesses_5sigma = excess_matching_significance_on_off(\ - n_off = n_background, alpha = alpha, significance = 5, method = 'lima') + stat = WStatCountsStatistic( + n_on=np.ones_like(n_background), + n_off=n_background, + alpha=alpha) + + n_excesses_5sigma = stat.excess_matching_significance(5) for i in range(0, n_bins_energy): # If the excess needed to get 5 sigma is less than 10, From ed52454e3bde42dba7e5a473f9367c42df3859c2 Mon Sep 17 00:00:00 2001 From: rlopezcoto Date: Thu, 2 Apr 2020 13:36:59 +0200 Subject: [PATCH 3/8] correct some tests --- lstchain/mc/tests/test_senstivity.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lstchain/mc/tests/test_senstivity.py b/lstchain/mc/tests/test_senstivity.py index e4dd690fae..b3d448e5a1 100644 --- a/lstchain/mc/tests/test_senstivity.py +++ b/lstchain/mc/tests/test_senstivity.py @@ -45,10 +45,10 @@ def test_calculate_sensitivity_lima(): np.testing.assert_allclose(calculate_sensitivity_lima( 50, 10, 0.2, 1, 0, 0), - ([13.48],[26.97]), rtol = 1.e-3) + ([13.48, 26.97]), rtol = 1.e-3) np.testing.assert_allclose(calculate_sensitivity_lima( 200, 50, 1, 0, 1, 0), - ([63.00],[31.5]), rtol = 1.e-3) + ([63.00, 31.5]), rtol = 1.e-3) # Testing an array np.testing.assert_allclose(calculate_sensitivity_lima( [10, 100], [50,100], 1, 1, 1, 0), From 8bda4554e7bacf9bad3cb8c5c6c521251435f26d Mon Sep 17 00:00:00 2001 From: rlopezcoto Date: Thu, 2 Apr 2020 15:14:41 +0200 Subject: [PATCH 4/8] solve test conflicts --- lstchain/mc/tests/test_senstivity.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lstchain/mc/tests/test_senstivity.py b/lstchain/mc/tests/test_senstivity.py index b3d448e5a1..86ba2a5d38 100644 --- a/lstchain/mc/tests/test_senstivity.py +++ b/lstchain/mc/tests/test_senstivity.py @@ -51,7 +51,7 @@ def test_calculate_sensitivity_lima(): ([63.00, 31.5]), rtol = 1.e-3) # Testing an array np.testing.assert_allclose(calculate_sensitivity_lima( - [10, 100], [50,100], 1, 1, 1, 0), + [10, 100], [50,100], [1, 1], 1, 1, 0), ([63.00, 83.57],[630.07, 83.57]), rtol = 1.e-3) def test_bin_definition(): From b9b82b0d8888d5848680533ae5640d3bce2ab4c3 Mon Sep 17 00:00:00 2001 From: rlopezcoto Date: Thu, 2 Apr 2020 16:11:47 +0200 Subject: [PATCH 5/8] Add tests for calculate_sensitivity_lima_ebin --- lstchain/mc/tests/test_senstivity.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/lstchain/mc/tests/test_senstivity.py b/lstchain/mc/tests/test_senstivity.py index 86ba2a5d38..2d1f4173b9 100644 --- a/lstchain/mc/tests/test_senstivity.py +++ b/lstchain/mc/tests/test_senstivity.py @@ -5,6 +5,7 @@ read_sim_par, calculate_sensitivity, calculate_sensitivity_lima, + calculate_sensitivity_lima_ebin, bin_definition, ring_containment, ) @@ -54,6 +55,18 @@ def test_calculate_sensitivity_lima(): [10, 100], [50,100], [1, 1], 1, 1, 0), ([63.00, 83.57],[630.07, 83.57]), rtol = 1.e-3) +def test_calculate_sensitivity_lima_ebin(): + + np.testing.assert_allclose(calculate_sensitivity_lima_ebin( + [50], [10], [0.2], 1), ([13.48], [26.97]), + rtol = 1.e-3) + + np.testing.assert_allclose(calculate_sensitivity_lima_ebin( + [50, 20, 10], [10, 10, 10], [0.2, 0.2, 0.2], 3), + (([13.48, 13.48, 13.48]), + [ 26.97208396, 67.43020989, 134.86041979]), + rtol = 1.e-3) + def test_bin_definition(): gammaness_bins, theta2_bins = bin_definition(11,10) From 90a315e6aa359825031e5bd913feefa789cd6b3d Mon Sep 17 00:00:00 2001 From: rlopezcoto Date: Thu, 2 Apr 2020 16:12:14 +0200 Subject: [PATCH 6/8] fix calls for arrays --- lstchain/mc/sensitivity.py | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/lstchain/mc/sensitivity.py b/lstchain/mc/sensitivity.py index 3fa047ea1d..7a91581ccf 100644 --- a/lstchain/mc/sensitivity.py +++ b/lstchain/mc/sensitivity.py @@ -207,9 +207,9 @@ def calculate_sensitivity_lima_ebin(n_excesses, n_background, alpha, n_bins_ener Parameters --------- n_excesses: `numpy.ndarray` number of excess events in the signal region - n_background: `numpy.ndarray` number of events in the background region - alpha: `float` inverse of the number of off positions - n_bins_energy: `int` number of bins in energy + n_background: `numpy.ndarray` number of events in the background region + alpha: `float` inverse of the number of off positions + n_bins_energy:`int` number of bins in energy Returns --------- @@ -218,6 +218,13 @@ def calculate_sensitivity_lima_ebin(n_excesses, n_background, alpha, n_bins_ener a 5 sigma significance """ + + if not (len(n_excesses) == n_bins_energy and + len(n_background) == n_bins_energy and + len(alpha) == n_bins_energy): + print("Excess, background and alpha arrays must have the same length") + return + stat = WStatCountsStatistic( n_on=np.ones_like(n_background), n_off=n_background, @@ -233,8 +240,8 @@ def calculate_sensitivity_lima_ebin(n_excesses, n_background, alpha, n_bins_ener # If the excess needed to get 5 sigma is less than 5% # of the background, we force it to be at least 5% of # the background - if n_excesses_5sigma[i] < 0.05 * n_background[i] * alpha: - n_excesses_5sigma[i] = 0.05 * n_background[i] * alpha + if n_excesses_5sigma[i] < 0.05 * n_background[i] * alpha[i]: + n_excesses_5sigma[i] = 0.05 * n_background[i] * alpha[i] sensitivity = n_excesses_5sigma / n_excesses * 100 # percentage of Crab From 1c08ca5c7f7c45b17bc48412d91e04d2e726a186 Mon Sep 17 00:00:00 2001 From: rlopezcoto Date: Thu, 2 Apr 2020 16:25:11 +0200 Subject: [PATCH 7/8] fix tests --- lstchain/mc/sensitivity.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/lstchain/mc/sensitivity.py b/lstchain/mc/sensitivity.py index 7a91581ccf..e877a3173d 100644 --- a/lstchain/mc/sensitivity.py +++ b/lstchain/mc/sensitivity.py @@ -453,8 +453,10 @@ def find_best_cuts_sensitivity(simtelfile_gammas, simtelfile_protons, ngamma_per_ebin[i] = np.sum(rate_weighted_g[(e_reco_g < energy[i+1]) & (e_reco_g > energy[i])]) * obstime nhadron_per_ebin[i] = np.sum(rate_weighted_p[(e_reco_p < energy[i+1]) & (e_reco_p > energy[i])]) * obstime - n_excesses_5sigma, sensitivity_3Darray = calculate_sensitivity_lima(final_gamma, final_hadrons * noff, 1/noff, - n_bins_energy, n_bins_gammaness, n_bins_theta2) + n_excesses_5sigma, sensitivity_3Darray = calculate_sensitivity_lima(final_gamma, final_hadrons * noff, + 1/noff * np.ones(len(final_gamma)), + n_bins_energy, n_bins_gammaness, + n_bins_theta2) # Avoid bins which are empty or have too few events: min_num_events = 10 @@ -675,7 +677,7 @@ def sensitivity(simtelfile_gammas, simtelfile_protons, * obstime.to(u.s).value n_excesses_5sigma, sensitivity_3Darray = calculate_sensitivity_lima_ebin(final_gamma, final_hadrons * noff, - 1 / noff, + 1 / noff * np.ones(len(final_gamma)), n_bins_energy) # Avoid bins which are empty or have too few events: min_num_events = 5 From 9d6f695c4e1b24a0a586cba9b8c842d6d655886d Mon Sep 17 00:00:00 2001 From: rlopezcoto Date: Fri, 3 Apr 2020 09:34:01 +0200 Subject: [PATCH 8/8] implement comments --- lstchain/mc/sensitivity.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/lstchain/mc/sensitivity.py b/lstchain/mc/sensitivity.py index e877a3173d..a8fe9763bc 100644 --- a/lstchain/mc/sensitivity.py +++ b/lstchain/mc/sensitivity.py @@ -219,11 +219,9 @@ def calculate_sensitivity_lima_ebin(n_excesses, n_background, alpha, n_bins_ener """ - if not (len(n_excesses) == n_bins_energy and - len(n_background) == n_bins_energy and - len(alpha) == n_bins_energy): - print("Excess, background and alpha arrays must have the same length") - return + if any(len(a) != n_bins_energy for a in (n_excesses, n_background, alpha)): + raise ValueError( + 'Excess, background and alpha arrays must have the same length') stat = WStatCountsStatistic( n_on=np.ones_like(n_background),