Fix HillasReconstructor in case of too few telescopes (#994)

* only counts telescopes with a width that is not np.NaN now

* using np.isnan to check for nan

* Exception had the wrong telescope number

* 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

* 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.

* adapting the test to account for changes in the master branch (r0 container: image->waveform)

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,6 +133,8 @@ 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
@@ -140,6 +146,15 @@ class are set to np.nan
                 "need at least two telescopes, have {}"
                 .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,
             inst.subarray,

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_large.simtel.gz")
+
+    fit = HillasReconstructor()
+
+    tel_azimuth = {}
+    tel_altitude = {}
+
+    source = event_source(filename, max_events=10)
+
+    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].waveform[0].sum(axis=1)
+
+            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)