Skip to content

Commit

Permalink
Adds a test to make sure the HillasReconstructor fails if
Browse files Browse the repository at this point in the history
- 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.
  • Loading branch information
LukasNickel committed May 22, 2019
1 parent a771f56 commit 7058a23
Showing 1 changed file with 70 additions and 1 deletion.
71 changes: 70 additions & 1 deletion ctapipe/reco/tests/test_HillasReconstructor.py
Original file line number Diff line number Diff line change
@@ -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

Expand Down Expand Up @@ -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)

0 comments on commit 7058a23

Please sign in to comment.