From 399ec788ad6cf0db11d7117865c35932905630a6 Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Thu, 28 Feb 2019 13:40:28 +0100 Subject: [PATCH 1/6] only counts telescopes with a width that is not np.NaN now --- ctapipe/reco/HillasReconstructor.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/ctapipe/reco/HillasReconstructor.py b/ctapipe/reco/HillasReconstructor.py index a5a66e79ec3..736fe0657dd 100644 --- a/ctapipe/reco/HillasReconstructor.py +++ b/ctapipe/reco/HillasReconstructor.py @@ -134,8 +134,10 @@ class are set to np.nan # filter warnings for missing obs time. this is needed because MC data has no obs time warnings.filterwarnings(action='ignore', category=MissingFrameAttributeWarning) - # stereoscopy needs at least two telescopes - if len(hillas_dict) < 2: + # stereoscopy needs at least two telescopes with a valid width of the hillas ellipse + valid_telescopes = sum([1 if hillas_dict[x]['width'].value is not np.NaN else 0 for x in hillas_dict]) + + if valid_telescopes < 2: raise TooFewTelescopesException( "need at least two telescopes, have {}" .format(len(hillas_dict))) From 83b7e8c29ccfd66d7eb4edcff730f1ce4be42ceb Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Thu, 28 Feb 2019 13:49:47 +0100 Subject: [PATCH 2/6] using np.isnan to check for nan --- ctapipe/reco/HillasReconstructor.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ctapipe/reco/HillasReconstructor.py b/ctapipe/reco/HillasReconstructor.py index 736fe0657dd..d0d8e60a5b6 100644 --- a/ctapipe/reco/HillasReconstructor.py +++ b/ctapipe/reco/HillasReconstructor.py @@ -135,7 +135,7 @@ class are set to np.nan warnings.filterwarnings(action='ignore', category=MissingFrameAttributeWarning) # stereoscopy needs at least two telescopes with a valid width of the hillas ellipse - valid_telescopes = sum([1 if hillas_dict[x]['width'].value is not np.NaN else 0 for x in hillas_dict]) + valid_telescopes = sum([1 if not np.isnan(hillas_dict[x]['width'].value) else 0 for x in hillas_dict]) if valid_telescopes < 2: raise TooFewTelescopesException( From 9d16b53b2ba25304a3ab8a40a0f8110018b8113f Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Thu, 28 Feb 2019 14:11:09 +0100 Subject: [PATCH 3/6] Exception had the wrong telescope number --- ctapipe/reco/HillasReconstructor.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/ctapipe/reco/HillasReconstructor.py b/ctapipe/reco/HillasReconstructor.py index d0d8e60a5b6..4d75f8dc969 100644 --- a/ctapipe/reco/HillasReconstructor.py +++ b/ctapipe/reco/HillasReconstructor.py @@ -133,14 +133,16 @@ class are set to np.nan # filter warnings for missing obs time. this is needed because MC data has no obs time warnings.filterwarnings(action='ignore', category=MissingFrameAttributeWarning) - - # stereoscopy needs at least two telescopes with a valid width of the hillas ellipse - valid_telescopes = sum([1 if not np.isnan(hillas_dict[x]['width'].value) else 0 for x in hillas_dict]) + + # stereoscopy needs at least two telescopes with valid widths + valid_telescopes = sum([1 if not np.isnan(hillas_dict[x]['width'].value) + else 0 + for x in hillas_dict]) if valid_telescopes < 2: raise TooFewTelescopesException( "need at least two telescopes, have {}" - .format(len(hillas_dict))) + .format(valid_telescopes)) self.initialize_hillas_planes( hillas_dict, From a771f56d18ca01e969050541dcb3dbd4695056b6 Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Mon, 11 Mar 2019 14:04:28 +0100 Subject: [PATCH 4/6] Raise Exceptions for invalid widths - Hillas Reconstroctur now fails in any of these cases: - len(hillas_dict) < 2 - any width is np.nan - any width is 0 --- ctapipe/reco/HillasReconstructor.py | 27 +++++++++++++++++++-------- 1 file changed, 19 insertions(+), 8 deletions(-) diff --git a/ctapipe/reco/HillasReconstructor.py b/ctapipe/reco/HillasReconstructor.py index 4d75f8dc969..138c27325c8 100644 --- a/ctapipe/reco/HillasReconstructor.py +++ b/ctapipe/reco/HillasReconstructor.py @@ -29,13 +29,17 @@ from astropy import units as u -__all__ = ['HillasReconstructor', 'TooFewTelescopesException', 'HillasPlane'] +__all__ = ['HillasReconstructor', 'TooFewTelescopesException', 'HillasPlane', 'InvalidWidthException'] class TooFewTelescopesException(Exception): pass +class InvalidWidthException(Exception): + pass + + def angle(v1, v2): """ computes the angle between two vectors assuming carthesian coordinates @@ -129,20 +133,27 @@ class are set to np.nan ------ TooFewTelescopesException if len(hillas_dict) < 2 + InvalidWidthException + if any width is np.nan or 0 ''' # filter warnings for missing obs time. this is needed because MC data has no obs time warnings.filterwarnings(action='ignore', category=MissingFrameAttributeWarning) - # stereoscopy needs at least two telescopes with valid widths - valid_telescopes = sum([1 if not np.isnan(hillas_dict[x]['width'].value) - else 0 - for x in hillas_dict]) - - if valid_telescopes < 2: + # stereoscopy needs at least two telescopes + if len(hillas_dict) < 2: raise TooFewTelescopesException( "need at least two telescopes, have {}" - .format(valid_telescopes)) + .format(len(hillas_dict))) + + # check for np.nan or 0 width's as these screw up weights + if any([np.isnan(hillas_dict[tel]['width'].value) for tel in hillas_dict]): + raise InvalidWidthException( + "A HillasContainer contains an ellipse of width==np.nan") + + if any([hillas_dict[tel]['width'].value == 0 for tel in hillas_dict]): + raise InvalidWidthException( + "A HillasContainer contains an ellipse of width==0") self.initialize_hillas_planes( hillas_dict, From 7058a2356a2e72f18d4975a03709d074e174421b Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Wed, 22 May 2019 12:10:00 +0200 Subject: [PATCH 5/6] Adds a test to make sure the HillasReconstructor fails if - there are less than 2 telescopes for the given event - any width is 0 - any width is nan Implementation is mostly copied from the test_reconstruction test. --- .../reco/tests/test_HillasReconstructor.py | 71 ++++++++++++++++++- 1 file changed, 70 insertions(+), 1 deletion(-) diff --git a/ctapipe/reco/tests/test_HillasReconstructor.py b/ctapipe/reco/tests/test_HillasReconstructor.py index 4dc1828c13a..236638bf0f4 100644 --- a/ctapipe/reco/tests/test_HillasReconstructor.py +++ b/ctapipe/reco/tests/test_HillasReconstructor.py @@ -1,10 +1,12 @@ import numpy as np from astropy import units as u +import pytest from ctapipe.image.cleaning import tailcuts_clean from ctapipe.image.hillas import hillas_parameters, HillasParameterizationError from ctapipe.io import event_source -from ctapipe.reco.HillasReconstructor import HillasReconstructor, HillasPlane +from ctapipe.reco.HillasReconstructor import ( + HillasReconstructor, HillasPlane, InvalidWidthException, TooFewTelescopesException) from ctapipe.utils import get_dataset_path from astropy.coordinates import SkyCoord, AltAz @@ -132,3 +134,70 @@ def test_reconstruction(): fit_result.az.to(u.deg) fit_result.core_x.to(u.m) assert fit_result.is_valid + + +def test_invalid_events(): + """ + The HillasReconstructor is supposed to fail + in these cases: + - less than two teleskopes + - any width is NaN + - any width is 0 + + This test uses the same sample simtel file as + test_reconstruction(). As there are no invalid events in this + file, multiple hillas_dicts are constructed to make sure + Exceptions get thrown in the mentioned edge cases. + + Test will fail if no Exception or another Exception gets thrown.""" + + filename = get_dataset_path("gamma_test.simtel.gz") + + fit = HillasReconstructor() + + tel_azimuth = {} + tel_altitude = {} + + source = event_source(filename) + + for event in source: + + hillas_dict = {} + for tel_id in event.dl0.tels_with_data: + + geom = event.inst.subarray.tel[tel_id].camera + tel_azimuth[tel_id] = event.mc.tel[tel_id].azimuth_raw * u.rad + tel_altitude[tel_id] = event.mc.tel[tel_id].altitude_raw * u.rad + + pmt_signal = event.r0.tel[tel_id].image[0] + + mask = tailcuts_clean(geom, pmt_signal, + picture_thresh=10., boundary_thresh=5.) + pmt_signal[mask == 0] = 0 + + try: + moments = hillas_parameters(geom, pmt_signal) + hillas_dict[tel_id] = moments + except HillasParameterizationError as e: + continue + + # construct a dict only containing the last telescope events + # (#telescopes < 2) + hillas_dict_only_one_tel = dict() + hillas_dict_only_one_tel[tel_id] = hillas_dict[tel_id] + with pytest.raises(TooFewTelescopesException): + fit.predict(hillas_dict_only_one_tel, event.inst, tel_azimuth, tel_altitude) + + # construct a hillas dict with the width of the last event set to 0 + # (any width == 0) + hillas_dict_zero_width = hillas_dict.copy() + hillas_dict_zero_width[tel_id]['width'] = 0 * u.m + with pytest.raises(InvalidWidthException): + fit.predict(hillas_dict_zero_width, event.inst, tel_azimuth, tel_altitude) + + # construct a hillas dict with the width of the last event set to np.nan + # (any width == nan) + hillas_dict_nan_width = hillas_dict.copy() + hillas_dict_zero_width[tel_id]['width'] = np.nan * u.m + with pytest.raises(InvalidWidthException): + fit.predict(hillas_dict_nan_width, event.inst, tel_azimuth, tel_altitude) From 8450731e4dedd938d975a4167ffc2533388831a9 Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Wed, 22 May 2019 13:27:38 +0200 Subject: [PATCH 6/6] adapting the test to account for changes in the master branch (r0 container: image->waveform) --- ctapipe/reco/tests/test_HillasReconstructor.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/ctapipe/reco/tests/test_HillasReconstructor.py b/ctapipe/reco/tests/test_HillasReconstructor.py index 236638bf0f4..97fc81bee33 100644 --- a/ctapipe/reco/tests/test_HillasReconstructor.py +++ b/ctapipe/reco/tests/test_HillasReconstructor.py @@ -151,14 +151,14 @@ def test_invalid_events(): Test will fail if no Exception or another Exception gets thrown.""" - filename = get_dataset_path("gamma_test.simtel.gz") + filename = get_dataset_path("gamma_test_large.simtel.gz") fit = HillasReconstructor() tel_azimuth = {} tel_altitude = {} - source = event_source(filename) + source = event_source(filename, max_events=10) for event in source: @@ -169,7 +169,7 @@ def test_invalid_events(): tel_azimuth[tel_id] = event.mc.tel[tel_id].azimuth_raw * u.rad tel_altitude[tel_id] = event.mc.tel[tel_id].altitude_raw * u.rad - pmt_signal = event.r0.tel[tel_id].image[0] + pmt_signal = event.r0.tel[tel_id].waveform[0].sum(axis=1) mask = tailcuts_clean(geom, pmt_signal, picture_thresh=10., boundary_thresh=5.)