Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Update gammapy stat #318

Merged
merged 8 commits into from
Apr 3, 2020
Merged

Update gammapy stat #318

Show file tree
Hide file tree
Changes from 7 commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
8511eb3
Modify excess_matching_significance from gammapy
rlopezcoto Apr 2, 2020
72fff95
correct also ebin function
rlopezcoto Apr 2, 2020
ed52454
correct some tests
rlopezcoto Apr 2, 2020
8bda455
solve test conflicts
rlopezcoto Apr 2, 2020
b9b82b0
Add tests for calculate_sensitivity_lima_ebin
rlopezcoto Apr 2, 2020
90a315e
fix calls for arrays
rlopezcoto Apr 2, 2020
1c08ca5
fix tests
rlopezcoto Apr 2, 2020
9d6f695
implement comments
rlopezcoto Apr 3, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
42 changes: 29 additions & 13 deletions lstchain/mc/sensitivity.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -176,9 +176,12 @@ 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):
Expand All @@ -204,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
---------
Expand All @@ -215,8 +218,19 @@ 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')

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
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

raise ValueError with that message, not print and return

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe if any(len(a) != n_bins_energy for a in (n_excesses, n_background, alpha)) as a more concise formulation

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks, already implemented. The only problem with this (also present with my formulation) is that (n_excesses, n_background, alpha) must be input as lists/arrays, and an error is raised if only single floats are given as input

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Then apply np.asanyarray(n_on, ndmin=1) to the input. If you also want to return a single number when the input were numbers, you can do sometehing like this:

scalar = np.isscalar(n_on)
n_on = np.asanyarray(n_on, ndmin=1)
n_off = ...
alpha = ...

...

if scalar:
     return result[0]
return result

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But let's get this merged first so the other PRs can run again


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,
Expand All @@ -226,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

Expand Down Expand Up @@ -439,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
Expand Down Expand Up @@ -661,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
Expand Down
19 changes: 16 additions & 3 deletions lstchain/mc/tests/test_senstivity.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
read_sim_par,
calculate_sensitivity,
calculate_sensitivity_lima,
calculate_sensitivity_lima_ebin,
bin_definition,
ring_containment,
)
Expand Down Expand Up @@ -45,15 +46,27 @@ 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),
[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)
Expand Down