diff --git a/ctapipe/coordinates/__init__.py b/ctapipe/coordinates/__init__.py index c5ba24ed797..2d9ab2a1de3 100644 --- a/ctapipe/coordinates/__init__.py +++ b/ctapipe/coordinates/__init__.py @@ -2,13 +2,6 @@ """ Coordinates. """ - - -from .coordinate_transformations import (pixel_position_to_direction, - alt_to_theta, az_to_phi, - transform_pixel_position) - - from .angular_frames import * from .ground_frames import * diff --git a/ctapipe/coordinates/angular_frames.py b/ctapipe/coordinates/angular_frames.py index 512a332ecdd..a772a816be2 100644 --- a/ctapipe/coordinates/angular_frames.py +++ b/ctapipe/coordinates/angular_frames.py @@ -215,8 +215,12 @@ def offset_to_altaz(x_off, y_off, azimuth, altitude): offset = sqrt(x_off * x_off + y_off * y_off) pos = np.where(offset == 0) # find offset 0 positions + if len(pos[0]) > 0: - offset[pos] = 1e-12 # add a very small offset to prevent math errors + if np.isscalar(offset): + offset = 1e-14 + else: + offset[pos] = 1e-14 # add a very small offset to prevent math errors atan_off = arctan(offset) @@ -236,8 +240,12 @@ def offset_to_altaz(x_off, y_off, azimuth, altitude): obj_azimuth = arctan2(yp0, -xp0) + azimuth if len(pos[0]) > 0: - obj_altitude[pos] = altitude - obj_azimuth[pos] = azimuth + if np.isscalar(obj_altitude): + obj_altitude = altitude + obj_azimuth = azimuth + else: + obj_altitude[pos] = altitude + obj_azimuth[pos] = azimuth obj_altitude = obj_altitude * u.rad obj_azimuth = obj_azimuth * u.rad diff --git a/ctapipe/coordinates/coordinate_transformations.py b/ctapipe/coordinates/coordinate_transformations.py deleted file mode 100644 index cbc2691b6e1..00000000000 --- a/ctapipe/coordinates/coordinate_transformations.py +++ /dev/null @@ -1,81 +0,0 @@ -import os - -import numpy as np -from astropy.utils.decorators import deprecated -from astropy import units as u - -from ctapipe import utils -from ctapipe.utils import linalg - - -@deprecated(0.1, "will be replaced with proper coord transform") -def pixel_position_to_direction(pix_x, pix_y, tel_phi, tel_theta, tel_foclen): - """ - TODO replace with proper implementation - calculates the direction vector of corresponding to a - (x,y) position on the camera - - beta is the pixel's angular distance to the centre - according to beta / tel_view = r / maxR - alpha is the polar angle between the y-axis and the pixel - to find the direction the pixel is looking at: - - - the pixel direction is set to the telescope direction - - offset by beta towards up - - rotated around the telescope direction by the angle alpha - - - Parameters - ----------- - pix_x, pix_y : ndarray - lists of x and y positions on the camera - tel_phi, tel_theta: astropy quantities - two angles that describe the orientation of the telescope - tel_foclen : astropy quantity - focal length of the telescope - - Returns - ------- - pix_dirs : ndarray - shape (n,3) list of "direction vectors" - corresponding to a position on the camera - """ - pix_alpha = np.arctan2(-pix_y, pix_x) - - pix_rho = (pix_x ** 2 + pix_y ** 2) ** .5 - - pix_beta = pix_rho / tel_foclen * u.rad - - tel_dir = linalg.set_phi_theta(tel_phi, tel_theta) - - pix_dirs = [] - for a, b in zip(pix_alpha, pix_beta): - pix_dir = linalg.set_phi_theta(tel_phi, tel_theta - b) - - pix_dir = linalg.rotate_around_axis(pix_dir, tel_dir, a) - pix_dirs.append(pix_dir * u.dimless) - - return pix_dirs - - -def alt_to_theta(alt): - """transforms altitude (angle from the horizon upwards) to theta (angle from z-axis) - for simtel array coordinate systems - """ - return (90 * u.deg - alt).to(alt.unit) - - -def az_to_phi(az): - """transforms azimuth (angle from north towards east) - to phi (angle from x-axis towards y-axis) - for simtel array coordinate systems - """ - return -az - - -def transform_pixel_position(x, y): - """transforms the x and y coordinates on the camera plane so that they correspond to - the view as if looked along the pointing direction of the telescope, i.e. y->up and - x->right - """ - return x, -y diff --git a/ctapipe/reco/HillasReconstructor.py b/ctapipe/reco/HillasReconstructor.py index e4996caf68e..3db7b6bc4c7 100644 --- a/ctapipe/reco/HillasReconstructor.py +++ b/ctapipe/reco/HillasReconstructor.py @@ -5,11 +5,10 @@ """ -from ctapipe.utils import linalg -from ctapipe.coordinates import coordinate_transformations as trafo from ctapipe.reco.reco_algorithms import Reconstructor from ctapipe.io.containers import ReconstructedShowerContainer - +from ctapipe.coordinates import TiltedGroundFrame, HorizonFrame, CameraFrame +from astropy.coordinates import SkyCoord, spherical_to_cartesian, cartesian_to_spherical from itertools import combinations import numpy as np @@ -17,115 +16,48 @@ from scipy.optimize import minimize from astropy import units as u -u.dimless = u.dimensionless_unscaled -__all__ = ['HillasReconstructor', - 'TooFewTelescopes', - 'dist_to_traces', 'MEst', 'GreatCircle'] +__all__ = ['HillasReconstructor', 'TooFewTelescopesException', 'HillasPlane'] -class TooFewTelescopes(Exception): +class TooFewTelescopesException(Exception): pass -def dist_to_traces(core, circles): - """This function calculates the M-Estimator from the distances of the - suggested core position to all traces of the given GreatCircles. - The distance of the core to the trace line is the length of the - vector between the core and an arbitrary point on the trace - projected perpendicular to the trace. - - This is implemented as the scalar product of the connecting vector - between the core and the position of the telescope and { trace[1], - -trace[0] } as the normal vector of the trace. - - Notes - ----- - uses the M-Estimator of the distance instead of the distance itself: - - .. math:: - - M_\\text{Est} = \sum_i{ 2 \cdot \sqrt{1 + d_i^2} - 2} - - - """ - - mest = 0. - for circle in circles.values(): - - # the distanece of the core - D = core - circle.pos[:2] / u.m - dist = D[0] * circle.trace[1] - D[1] * circle.trace[0] - - # summing up the M-Estimator with the given circle weights - mest += (2 * np.sqrt(1 + (dist ** 2)) - 2) * circle.weight - return mest - - -def MEst(origin, circles, weights): - """calculates the M-Estimator: a modified χ² that becomes - asymptotically linear for high values and is therefore less - sensitive to outliers - - the test is performed to maximise the angles between the fit - direction and the all the normal vectors of the great circles - - .. math:: - - M_\\text{Est} = \sum_i{ 2 \cdot \sqrt{1 + d_i^2} - 2} - - - Notes - ----- - seemingly inferior to negative sum of sin(angle)... - +def angle(v1, v2): + """ computes the angle between two vectors + assuming carthesian coordinates Parameters - ----------- - origin : length-3 array - direction vector of the gamma's origin used as seed - circles : GreatCircle array - collection of great circles created from the camera images - weights : array - list of weights for each image/great circle + ---------- + v1 : numpy array + v2 : numpy array Returns ------- - MEstimator : float - - + the angle between v1 and v2 as a dimensioned astropy quantity """ + norm = np.linalg.norm(v1) * np.linalg.norm(v2) + return np.arccos(np.clip(v1.dot(v2) / norm, -1.0, 1.0)) - sin_ang = np.array([linalg.length(np.cross(origin, circ.norm)) - for circ in circles.values()]) - return -2 * np.sum(weights * np.sqrt((1 + np.square(sin_ang))) - 2) - -def neg_angle_sum(origin, circles, weights): - """calculates the negative sum of the angle between the fit direction - and all the normal vectors of the great circles +def normalise(vec): + """ Sets the length of the vector to 1 + without changing its direction Parameters - ----------- - origin : length-3 array - direction vector of the gamma's origin used as seed - circles : GreatCircle array - collection of great circles created from the camera images - weights : array - list of weights for each image/great circle + ---------- + vec : numpy array Returns - -------- - n_sum_angles : float - negative of the sum of the angles between the test direction - and all normal vectors of the given great circles - + ------- + numpy array with the same direction but length of 1 """ - - sin_ang = np.array([linalg.length(np.cross(origin, circ.norm)) - for circ in circles.values()]) - return -np.sum(weights * sin_ang) + try: + return vec / np.linalg.norm(vec) + except ZeroDivisionError: + return vec class HillasReconstructor(Reconstructor): @@ -142,9 +74,9 @@ class that reconstructs the direction of an atmospheric shower def __init__(self, config=None, tool=None, **kwargs): super().__init__(config=config, parent=tool, **kwargs) - self.circles = {} + self.hillas_planes = {} - def predict(self, hillas_dict, inst, tel_phi, tel_theta, seed_pos=(0, 0)): + def predict(self, hillas_dict, inst, pointing_alt, pointing_az, seed_pos=(0, 0)): """The function you want to call for the reconstruction of the event. It takes care of setting up the event and consecutively calls the functions for the direction and core position @@ -158,8 +90,8 @@ class are set to np.nan MomentParameters instances as values inst : ctapipe.io.InstrumentContainer instrumental description - tel_phi: - tel_theta: + pointing_alt: + pointing_az: seed_pos : python tuple shape (2) tuple with a possible seed for the core position fit (e.g. CoG of all telescope images) @@ -177,27 +109,28 @@ class are set to np.nan "need at least two telescopes, have {}" .format(len(hillas_dict))) - self.get_great_circles(hillas_dict, inst.subarray, tel_phi, tel_theta) + self.inititialize_hillas_planes( + hillas_dict, + inst.subarray, + pointing_alt, + pointing_az + ) # algebraic direction estimate - direction, err_est_dir = self.fit_origin_crosses() + direction, err_est_dir = self.estimate_direction() # core position estimate using a geometric approach - pos, err_est_pos = self.fit_core_crosses() - - # numerical minimisations do not really improve the fit - # direction estimate using numerical minimisation - # dir = self.fit_origin_minimise(dir) - # - # core position estimate using numerical minimisation - # pos = self.fit_core_minimise(seed_pos) + pos, err_est_pos = self.estimate_core_position() # container class for reconstructed showers result = ReconstructedShowerContainer() - phi, theta = linalg.get_phi_theta(direction).to(u.deg) + _, lat, lon = cartesian_to_spherical(*direction) + + # estimate max height of shower + h_max = self.estimate_h_max(hillas_dict, inst.subarray, pointing_alt, pointing_az) - # TODO fix coordinates! - result.alt, result.az = 90 * u.deg - theta, -phi + + result.alt, result.az = lat, lon result.core_x = pos[0] result.core_y = pos[1] result.core_uncert = err_est_pos @@ -208,15 +141,23 @@ class are set to np.nan result.alt_uncert = err_est_dir result.az_uncert = np.nan - result.h_max = self.fit_h_max(hillas_dict, inst.subarray, tel_phi, tel_theta) + + result.h_max = h_max result.h_max_uncert = np.nan + result.goodness_of_fit = np.nan return result - def get_great_circles(self, hillas_dict, subarray, tel_phi, tel_theta): + def inititialize_hillas_planes( + self, + hillas_dict, + subarray, + pointing_alt, + pointing_az + ): """ - creates a dictionary of :class:`.GreatCircle` from a dictionary of + creates a dictionary of :class:`.HillasPlane` from a dictionary of hillas parameters @@ -231,38 +172,56 @@ def get_great_circles(self, hillas_dict, subarray, tel_phi, tel_theta): needs to contain at least the same keys as in `hillas_dict` """ - self.circles = {} + self.hillas_planes = {} for tel_id, moments in hillas_dict.items(): + p2_x = moments.cen_x + 0.1 * u.m * np.cos(moments.psi) + p2_y = moments.cen_y + 0.1 * u.m * np.sin(moments.psi) + focal_length = subarray.tel[tel_id].optics.equivalent_focal_length + + pointing = SkyCoord( + alt=pointing_alt[tel_id], + az=pointing_az[tel_id], + frame='altaz' + ) + + hf = HorizonFrame( + array_direction=pointing, + pointing_direction=pointing + ) + cf = CameraFrame( + focal_length=focal_length, + array_direction=pointing, + pointing_direction=pointing + ) + + cog_coord = SkyCoord(x=moments.cen_x, y=moments.cen_y, frame=cf) + cog_coord = cog_coord.transform_to(hf) + + p2_coord = SkyCoord(x=p2_x, y=p2_y, frame=cf) + p2_coord = p2_coord.transform_to(hf) - # NOTE this is correct: +cos(psi) ; +sin(psi) - p2_x = moments.cen_x + moments.length * np.cos(moments.psi) - p2_y = moments.cen_y + moments.length * np.sin(moments.psi) - foclen = subarray.tel[tel_id].optics.equivalent_focal_length - - circle = GreatCircle( - trafo.pixel_position_to_direction( - np.array([moments.cen_x / u.m, p2_x / u.m]) * u.m, - np.array([moments.cen_y / u.m, p2_y / u.m]) * u.m, - tel_phi[tel_id], tel_theta[tel_id], foclen), - moments.size * (moments.length / moments.width) + circle = HillasPlane( + p1=cog_coord, + p2=p2_coord, + telescope_position=subarray.positions[tel_id], + weight=moments.size * (moments.length / moments.width), ) - circle.pos = subarray.positions[tel_id] - self.circles[tel_id] = circle + self.hillas_planes[tel_id] = circle - def fit_origin_crosses(self): + def estimate_direction(self): """calculates the origin of the gamma as the weighted average - direction of the intersections of all great circles + direction of the intersections of all hillas planes Returns ------- gamma : shape (3) numpy array direction of origin of the reconstructed shower as a 3D vector crossings : shape (n,3) list - list of all the crossings of the `GreatCircle` list + an error esimate """ crossings = [] - for perm in combinations(self.circles.values(), 2): + for perm in combinations(self.hillas_planes.values(), 2): n1, n2 = perm[0].norm, perm[1].norm # cross product automatically weighs in the angle between # the two vectors: narrower angles have less impact, @@ -277,65 +236,25 @@ def fit_origin_crosses(self): crossing *= -1 crossings.append(crossing * perm[0].weight * perm[1].weight) - result = linalg.normalise(np.sum(crossings, axis=0)) * u.dimless - off_angles = [linalg.angle(result, cross) / u.rad for cross in crossings] + result = normalise(np.sum(crossings, axis=0)) + off_angles = [angle(result, cross) for cross in crossings] * u.rad + err_est_dir = np.average( off_angles, weights=[len(cross) for cross in crossings] - ) * u.rad + ) - # averaging over the solutions of all permutations return result, err_est_dir - def fit_origin_minimise(self, seed=(0, 0, 1), test_function=neg_angle_sum): - """ Fits the origin of the gamma with a minimisation procedure this - function expects that `get_great_circles` - has been run already. A seed should be given otherwise it defaults to - "straight up" supperted functions to minimise are an M-estimator and the - negative sum of the angles to all normal vectors of the great - circles - - Parameters - ----------- - seed : length-3 array - starting point of the minimisation - test_function : function object, optional (default: neg_angle_sum) - either neg_angle_sum or MEst (or own implementation...) - neg_angle_sum seemingly superior to MEst - - Returns - ------- - direction : length-3 numpy array as dimensionless quantity - best fit for the origin of the gamma from - the minimisation process - - """ - - # using the sum of the cosines of each direction with every - # other direction as weight; don't use the product -- with many - # steep angles, the product will become too small and the - # weight (and the whole fit) useless - weights = [np.sum([linalg.length(np.cross(A.norm, B.norm)) - for A in self.circles.values()] - ) * B.weight - for B in self.circles.values()] - - # minimising the test function - self.fit_result_origin = minimize(test_function, seed, - args=(self.circles, weights), - method='BFGS', - options={'disp': False} - ) - return np.array(linalg.normalise(self.fit_result_origin.x)) * u.dimless - def fit_core_crosses(self): + def estimate_core_position(self): r"""calculates the core position as the least linear square solution of an (over-constrained) equation system Notes ----- - The basis is the "trace" of each telescope's `GreatCircle` which + The basis is the "trace" of each telescope's `HillasPlane` which can be determined by the telescope's position P=(Px, Py) and the circle's normal vector, projected to the ground n=(nx, ny), so that for every r=(x, y) on the trace @@ -400,9 +319,9 @@ def fit_core_crosses(self): """ - A = np.zeros((len(self.circles), 2)) - D = np.zeros(len(self.circles)) - for i, circle in enumerate(self.circles.values()): + A = np.zeros((len(self.hillas_planes), 2)) + D = np.zeros(len(self.hillas_planes)) + for i, circle in enumerate(self.hillas_planes.values()): # apply weight from circle and from the tilt of the circle # towards the horizontal plane: simply projecting # circle.norm to the ground gives higher weight to planes @@ -421,78 +340,57 @@ def fit_core_crosses(self): # instead used directly the numpy implementation # speed is the same, just handles already "SingularMatrixError" - pos = np.linalg.lstsq(A, D)[0] * u.m + if np.all(np.isfinite(A)) and np.all(np.isfinite(D)): + # note that NaN values create a value error with MKL + # installations but not otherwise. + pos = np.linalg.lstsq(A, D, rcond=None)[0] * u.m + else: + return [np.nan, np.nan], [np.nan, np.nan] weighted_sum_dist = np.sum([np.dot(pos[:2] - c.pos[:2], c.norm[:2]) * c.weight - for c in self.circles.values()]) * pos.unit - norm_sum_dist = np.sum([c.weight * linalg.length(c.norm[:2]) - for c in self.circles.values()]) + for c in self.hillas_planes.values()]) * pos.unit + norm_sum_dist = np.sum([c.weight * np.linalg.norm(c.norm[:2]) + for c in self.hillas_planes.values()]) pos_uncert = abs(weighted_sum_dist / norm_sum_dist) return pos, pos_uncert - def fit_core_minimise(self, seed=(0, 0), test_function=dist_to_traces): - """ - reconstructs the shower core position from the already set up - great circles - Notes - ----- - The core of the shower lies on the cross section of the great - circle with the horizontal plane. The direction of this cross - section is the cross-product of the circle's normal vector and - the horizontal plane. Here, we only care about the direction; - not the orientation... + def estimate_h_max(self, hillas_dict, subarray, pointing_alt, pointing_az): + weights = [] + tels = [] + dirs = [] - Parameters - ---------- - seed : tuple - shape (2) tuple with optional starting coordinates - tuple of floats or astropy.length -- if floats, assume metres - test_function : function object, optional (default: dist_to_traces) - function to be used by the minimiser - - """ - - if type(seed) == u.Quantity: - unit = seed.unit - else: - unit = u.m + for tel_id, moments in hillas_dict.items(): - # the direction of the cross section of the great circle with - # the horizontal frame is the cross product of the great - # circle's normal vector with the z-axis: - # n × z = (n[1], -n[0], 0) - for circle in self.circles.values(): - circle.trace = linalg.normalise(np.array([circle.norm[1], - -circle.norm[0], 0])) + focal_length = subarray.tel[tel_id].optics.equivalent_focal_length - # minimising the test function (note: minimize strips seed off its - # unit) - self.fit_result_core = minimize(test_function, seed[:2], - args=(self.circles), - method='BFGS', options={'disp': False}) + pointing = SkyCoord( + alt=pointing_alt[tel_id], + az=pointing_az[tel_id], + frame='altaz' + ) - if not self.fit_result_core.success: - print("fit_core: fit no success") + hf = HorizonFrame( + array_direction=pointing, + pointing_direction=pointing + ) + cf = CameraFrame( + focal_length=focal_length, + array_direction=pointing, + pointing_direction=pointing + ) - return np.array(self.fit_result_core.x) * unit + cog_coord = SkyCoord(x=moments.cen_x, y=moments.cen_y, frame=cf) + cog_coord = cog_coord.transform_to(hf) - def fit_h_max(self, hillas_dict, subarray, tel_phi, tel_theta): + cog_direction = spherical_to_cartesian(1, cog_coord.alt, cog_coord.az) + cog_direction = np.array(cog_direction).ravel() - weights = [] - tels = [] - dirs = [] - for tel_id, hillas in hillas_dict.items(): - foclen = subarray.tel[tel_id].optics.equivalent_focal_length - max_dir, = trafo.pixel_position_to_direction( - np.array([hillas.cen_x / u.m]) * u.m, - np.array([hillas.cen_y / u.m]) * u.m, - tel_phi[tel_id], tel_theta[tel_id], foclen) - weights.append(self.circles[tel_id].weight) - tels.append(self.circles[tel_id].pos) - dirs.append(max_dir) + weights.append(self.hillas_planes[tel_id].weight) + tels.append(self.hillas_planes[tel_id].pos) + dirs.append(cog_direction) # minimising the test function pos_max = minimize(dist_to_line3d, np.array([0, 0, 10000]), @@ -509,23 +407,34 @@ def dist_to_line3d(pos, tels, dirs, weights): return result -class GreatCircle: +class HillasPlane: """ a tiny helper class to collect some parameters for each great great circle + + Stores some vectors a, b, and c + + These vectors are eucledian [x, y, z] where positive z values point towards the sky + and x and y are parallel to the ground. """ - def __init__(self, dirs, weight=1): - """the constructor takes two directions on the circle and creates the - normal vector belonging to that plane + def __init__(self, p1, p2, telescope_position, weight=1): + """The constructor takes two coordinates in the horizontal + frame (alt, az) which define a plane perpedicular + to the camera. Parameters ----------- - dirs : shape (2,3) ndarray - contains two 3D direction-vectors + p1: astropy.coordinates.SkyCoord + One of two direction vectors which define the plane. + This coordinate has to be defined in the ctapipe.coordinates.HorizonFrame + p2: astropy.coordinates.SkyCoord + One of two direction vectors which define the plane. + This coordinate has to be defined in the ctapipe.coordinates.HorizonFrame + telescope_position: np.array(3) + Position of the telescope on the ground weight : float, optional - weight of this telescope for later use during the - reconstruction + weight of this plane for later use during the reconstruction Notes ----- @@ -539,15 +448,17 @@ def __init__(self, dirs, weight=1): perpendicular to a, b and c """ - self.a = dirs[0] - self.b = dirs[1] + self.pos = telescope_position + + self.a = np.array(spherical_to_cartesian(1, p1.alt, p1.az)).ravel() + self.b = np.array(spherical_to_cartesian(1, p2.alt, p2.az)).ravel() # a and c form an orthogonal basis for the great circle # not really necessary since the norm can be calculated # with a and b just as well self.c = np.cross(np.cross(self.a, self.b), self.a) # normal vector for the plane defined by the great circle - self.norm = linalg.normalise(np.cross(self.a, self.c)) + self.norm = normalise(np.cross(self.a, self.c)) # some weight for this circle # (put e.g. uncertainty on the Hillas parameters # or number of PE in here) diff --git a/ctapipe/reco/tests/test_HillasReconstructor.py b/ctapipe/reco/tests/test_HillasReconstructor.py index 8f2a4275ac8..2a579fade70 100644 --- a/ctapipe/reco/tests/test_HillasReconstructor.py +++ b/ctapipe/reco/tests/test_HillasReconstructor.py @@ -3,94 +3,52 @@ from ctapipe.image.cleaning import tailcuts_clean from ctapipe.image.hillas import hillas_parameters, HillasParameterizationError -from ctapipe.io.hessio import hessio_event_source -from ctapipe.reco.HillasReconstructor import HillasReconstructor, GreatCircle +from ctapipe.io.eventsourcefactory import EventSourceFactory +from ctapipe.reco.HillasReconstructor import HillasReconstructor, HillasPlane from ctapipe.utils import get_dataset_path +from astropy.coordinates import SkyCoord -def test_fit_core(): +def test_estimator_results(): """ - creating some great circles pointing in different directions (two - north-south, - two east-west) and that have a slight position errors (+- 0.1 m in one of - the four - cardinal directions """ - circle1 = GreatCircle([[1, 0, 0], [0, 0, 1]]) - circle1.pos = [0, 0.1] * u.m - circle1.trace = [1, 0, 0] - - circle2 = GreatCircle([[0, 1, 0], [0, 0, 1]]) - circle2.pos = [0.1, 0] * u.m - circle2.trace = [0, 1, 0] - - circle3 = GreatCircle([[1, 0, 0], [0, 0, 1]]) - circle3.pos = [0, -.1] * u.m - circle3.trace = [1, 0, 0] - - circle4 = GreatCircle([[0, 1, 0], [0, 0, 1]]) - circle4.pos = [-.1, 0] * u.m - circle4.trace = [0, 1, 0] - - # creating the fit class and setting the the great circle member - fit = HillasReconstructor() - fit.circles = {1: circle1, 2: circle2, 3: circle3, 4: circle4} - - # performing the position fit with the minimisation algorithm - # and a seed that is quite far away - pos_fit_minimise = fit.fit_core_minimise([100, 1000] * u.m) - print("position fit test minimise:", pos_fit_minimise) - print() - - # performing the position fit with the geometric algorithm - pos_fit_crosses, err_est_pos_fit_crosses = fit.fit_core_crosses() - print("position fit test crosses:", pos_fit_crosses) - print("error estimate:", err_est_pos_fit_crosses) - print() - - # the results should be close to the origin of the coordinate system - np.testing.assert_allclose(pos_fit_minimise / u.m, [0, 0], atol=1e-3) - np.testing.assert_allclose(pos_fit_crosses / u.m, [0, 0], atol=1e-3) - - -def test_fit_origin(): - """ - creating some great circles pointing in different directions (two + creating some planes pointing in different directions (two north-south, two east-west) and that have a slight position errors (+- 0.1 m in one of the four cardinal directions """ - circle1 = GreatCircle([[1, 0, 0], [0, 0, 1]]) - circle1.pos = [0, 0.1] * u.m - circle1.trace = [1, 0, 0] - circle2 = GreatCircle([[0, 1, 0], [0, 0, 1]]) - circle2.pos = [0.1, 0] * u.m - circle2.trace = [0, 1, 0] + p1 = SkyCoord(alt=43 * u.deg, az=45 * u.deg, frame='altaz') + p2 = SkyCoord(alt=47 * u.deg, az=45 * u.deg, frame='altaz') + circle1 = HillasPlane(p1=p1, p2=p2, telescope_position=[0, 1, 0] * u.m) + + p1 = SkyCoord(alt=44 * u.deg, az=90 * u.deg, frame='altaz') + p2 = SkyCoord(alt=46 * u.deg, az=90 * u.deg, frame='altaz') + circle2 = HillasPlane(p1=p1, p2=p2, telescope_position=[1, 0, 0] * u.m) - circle3 = GreatCircle([[1, 0, 0], [0, 0, 1]]) - circle3.pos = [0, -.1] * u.m - circle3.trace = [1, 0, 0] + p1 = SkyCoord(alt=44.5 * u.deg, az=45 * u.deg, frame='altaz') + p2 = SkyCoord(alt=46.5 * u.deg, az=45 * u.deg, frame='altaz') + circle3 = HillasPlane(p1=p1, p2=p2, telescope_position=[0, -1, 0] * u.m) - circle4 = GreatCircle([[0, 1, 0], [0, 0, 1]]) - circle4.pos = [-.1, 0] * u.m - circle4.trace = [0, 1, 0] + p1 = SkyCoord(alt=43.5 * u.deg, az=90 * u.deg, frame='altaz') + p2 = SkyCoord(alt=45.5 * u.deg, az=90 * u.deg, frame='altaz') + circle4 = HillasPlane(p1=p1, p2=p2, telescope_position=[-1, 0, 0] * u.m) # creating the fit class and setting the the great circle member fit = HillasReconstructor() - fit.circles = {1: circle1, 2: circle2, 3: circle3, 4: circle4} + fit.hillas_planes = {1: circle1, 2: circle2, 3: circle3, 4: circle4} # performing the direction fit with the minimisation algorithm # and a seed that is perpendicular to the up direction - dir_fit_minimise = fit.fit_origin_minimise((0.1, 0.1, 1)) + dir_fit_minimise, _ = fit.estimate_direction() print("direction fit test minimise:", dir_fit_minimise) print() # performing the direction fit with the geometric algorithm - dir_fit_crosses = fit.fit_origin_crosses()[0] - print("direction fit test crosses:", dir_fit_crosses) + fitted_core_position, _ = fit.estimate_core_position() + print("direction fit test core position:", fitted_core_position) print() # the results should be close to the direction straight up - # np.testing.assert_allclose(dir_fit_minimise, [0, 0, 1], atol=1e-1) - np.testing.assert_allclose(dir_fit_crosses, [0, 0, 1], atol=1e-3) + np.testing.assert_allclose(dir_fit_minimise, [0, 0, 1], atol=1e-3) + np.testing.assert_allclose(fitted_core_position.value, [0, 0], atol=1e-3) def test_reconstruction(): @@ -98,7 +56,7 @@ def test_reconstruction(): a test of the complete fit procedure on one event including: • tailcut cleaning • hillas parametrisation - • GreatCircle creation + • HillasPlane creation • direction fit • position fit @@ -108,10 +66,13 @@ def test_reconstruction(): fit = HillasReconstructor() - tel_phi = {} - tel_theta = {} + tel_azimuth = {} + tel_altitude = {} - source = hessio_event_source(filename) + source = EventSourceFactory.produce( + input_url=filename, + product='HESSIOEventSource', + ) for event in source: @@ -119,9 +80,8 @@ def test_reconstruction(): for tel_id in event.dl0.tels_with_data: geom = event.inst.subarray.tel[tel_id].camera - tel_phi[tel_id] = event.mc.tel[tel_id].azimuth_raw * u.rad - tel_theta[tel_id] = (np.pi / 2 - - event.mc.tel[tel_id].altitude_raw) * u.rad + 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] @@ -139,7 +99,7 @@ def test_reconstruction(): if len(hillas_dict) < 2: continue - fit_result = fit.predict(hillas_dict, event.inst, tel_phi, tel_theta) + fit_result = fit.predict(hillas_dict, event.inst, tel_azimuth, tel_altitude) print(fit_result) fit_result.alt.to(u.deg) diff --git a/ctapipe/utils/linalg.py b/ctapipe/utils/linalg.py index 84b6410e623..30ed088a7de 100644 --- a/ctapipe/utils/linalg.py +++ b/ctapipe/utils/linalg.py @@ -4,8 +4,7 @@ import numpy as np from numpy import cos, sin, arctan2 as atan2, arccos as acos -__all__ = ['rotate_around_axis', 'rotation_matrix_2d', 'length', 'normalise', - 'angle', 'set_phi_theta', 'set_phi_theta_r'] +__all__ = ['rotation_matrix_2d', 'length', 'normalise', 'angle'] def rotation_matrix_2d(angle): @@ -18,35 +17,6 @@ def rotation_matrix_2d(angle): [sin(psi), cos(psi)]]) -def rotate_around_axis(vec, axis, angle): - """ rotates a vector around an axis by an angle - creates a rotation matrix and multiplies - the initial vector with it - - Parameters - ---------- - vec : length-3 numpy array - 3D vector to be rotated - axis : length-3 numpy array - axis around which the rotation is performed - angle : astropy angle quantity or float - angle (in rad if float) by which vec is rotated around axis - - Returns - ------- - rotated numpy array - """ - - axis = normalise(axis) - a = cos(angle / 2.0) - b, c, d = -axis * sin(angle / 2.0) - aa, bb, cc, dd = a * a, b * b, c * c, d * d - bc, ad, ac, ab, bd, cd = b * c, a * d, a * c, a * b, b * d, c * d - rot_matrix = np.array([[aa + bb - cc - dd, 2 * (bc + ad), 2 * (bd - ac)], - [2 * (bc - ad), aa + cc - bb - dd, 2 * (cd + ab)], - [2 * (bd + ac), 2 * (cd - ab), aa + dd - bb - cc]]) - return vec.dot(rot_matrix) - def length(vec): """ returns the length/norm of a numpy array @@ -86,37 +56,5 @@ def angle(v1, v2): ------- the angle between v1 and v2 as a dimensioned astropy quantity """ - return acos(np.clip(v1.dot(v2) / (length(v1) * length(v2)), -1.0, 1.0)) - - -def set_phi_theta_r(phi, theta, r=1): - """ sets a 3D vector according to the given angles - - Parameters - ---------- - phi : astropy.Quantity - theta : astropy.Quantity - r : (optional) - the length of the vector - can have a unit, doesn't have to - - Returns - ------- - numpy array with the given direction and length - """ - return np.array([sin(theta) * cos(phi), - sin(theta) * sin(phi), - cos(theta)]) * r - - -""" simple alias for set_phi_theta_r """ -set_phi_theta = set_phi_theta_r - - -def get_phi_theta(vec): - """ returns a tupel of the phi and theta angles of the given vector - """ - try: - return (atan2(vec[1], vec[0]), acos(np.clip(vec[2] / length(vec), -1, 1))) * u.rad - except ValueError: - return (0, 0) + norm = np.linalg.norm(v1) * np.linalg.norm(v2) + return np.arccos(np.clip(v1.dot(v2) / norm, -1.0, 1.0)) diff --git a/examples/notebooks/Theta square plot.ipynb b/examples/notebooks/Theta square plot.ipynb index 8f42ddeff13..8729f5431ae 100644 --- a/examples/notebooks/Theta square plot.ipynb +++ b/examples/notebooks/Theta square plot.ipynb @@ -4,22 +4,32 @@ "cell_type": "code", "execution_count": 1, "metadata": { + "ExecuteTime": { + "end_time": "2018-06-15T12:49:35.515499Z", + "start_time": "2018-06-15T12:49:34.968051Z" + }, "collapsed": true }, "outputs": [], "source": [ - "%matplotlib notebook" + "%matplotlib inline" ] }, { "cell_type": "code", - "execution_count": 1, + "execution_count": 2, "metadata": { + "ExecuteTime": { + "end_time": "2018-06-15T12:49:37.807612Z", + "start_time": "2018-06-15T12:49:35.520552Z" + }, "collapsed": true }, "outputs": [], "source": [ "from astropy import units as u\n", + "from astropy.coordinates.angle_utilities import angular_separation\n", + "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", @@ -36,26 +46,14 @@ }, { "cell_type": "code", - "execution_count": 2, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "env: CTAPIPE_SVC_PATH=/Users/jer/DATA/MAGIC/prod3blp_sample\n" - ] - } - ], - "source": [ - "# path to MC dataset\n", - "%env CTAPIPE_SVC_PATH=/Users/jer/DATA/MAGIC/prod3blp_sample" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": {}, + "execution_count": 3, + "metadata": { + "ExecuteTime": { + "end_time": "2018-06-15T12:49:37.815579Z", + "start_time": "2018-06-15T12:49:37.810814Z" + }, + "collapsed": true + }, "outputs": [], "source": [ "# MC dataset\n", @@ -64,139 +62,131 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 4, "metadata": { - "collapsed": true + "ExecuteTime": { + "end_time": "2018-06-15T12:49:37.887391Z", + "start_time": "2018-06-15T12:49:37.818824Z" + } }, "outputs": [], "source": [ "# get source events in MC dataset\n", - "source = event_source(filename, allowed_tels={1, 2, 3, 4})" + "source = event_source(filename, allowed_tels={1, 2, 3, 4})\n", + "reco = HillasReconstructor()\n", + "calib = CameraCalibrator(r1_product='HESSIOR1Calibrator')" ] }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 5, "metadata": { - "collapsed": true + "ExecuteTime": { + "end_time": "2018-06-15T12:49:47.500199Z", + "start_time": "2018-06-15T12:49:37.893169Z" + } }, - "outputs": [], - "source": [ - "# for each event\n", - "off_angles = [] \n", - "\n", - "reco = HillasReconstructor()\n", - "calib = CameraCalibrator()" - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ - "WARNING: AstropyDeprecationWarning: will be replaced with proper coord transform [ctapipe.reco.HillasReconstructor]\n", - "WARNING:astropy:AstropyDeprecationWarning: will be replaced with proper coord transform\n", - "/Users/kosack/Projects/CTA/Working/ctapipe/ctapipe/image/hillas.py:582: RuntimeWarning: invalid value encountered in sqrt\n", + "/Users/mackaiver/Development/ctapipe/ctapipe/image/hillas.py:571: RuntimeWarning: invalid value encountered in sqrt\n", + " width = np.sqrt((vy2 + vx2 - z) / 2.0)\n", + "/Users/mackaiver/Development/ctapipe/ctapipe/image/hillas.py:571: RuntimeWarning: invalid value encountered in sqrt\n", " width = np.sqrt((vy2 + vx2 - z) / 2.0)\n" ] } ], "source": [ + "off_angles = []\n", "for event in source:\n", "\n", - " # shower direction\n", - " # converting MC shower angular parameters to 3D spatial vector\n", - " shower_azimuth = event.mc.az # same as in MC file\n", - " shower_altitude = np.pi * u.rad / 2 - event.mc.alt # altitude = 90 - mc.alt\n", - " shower_direction = linalg.set_phi_theta(shower_azimuth, shower_altitude)\n", - "\n", " # calibrating the event\n", " calib.calibrate(event)\n", - "\n", - " # for each telescope and event\n", + " \n", + " hillas_params = {}\n", + " # pointing direction of the telescopes\n", " point_azimuth = {}\n", " point_altitude = {}\n", - " hillas_params = {}\n", - "\n", + " \n", + " subarray = event.inst.subarray\n", + " \n", " # get hillas params for each event in different telescopes\n", " for tel_id in event.dl0.tels_with_data:\n", "\n", - " # telescope pointing direction extracted from MC dataset\n", + " # telescope pointing direction\n", " point_azimuth[tel_id] = event.mc.tel[tel_id].azimuth_raw * u.rad\n", - " point_altitude[tel_id] = (\n", - " np.pi / 2 - event.mc.tel[tel_id].altitude_raw) * u.rad\n", + " point_altitude[tel_id] = event.mc.tel[tel_id].altitude_raw * u.rad\n", + " # print(point_azimuth,point_altitude)\n", "\n", - " # camera geometry required for hillas parametrization\n", - " camgeom = event.inst.subarray.tel[tel_id].camera\n", + " # Camera Geometry required for hillas parametrization\n", + " camgeom = subarray.tel[tel_id].camera\n", "\n", - " # note that [0] is for channel 0 which is high-gain channel\n", + " # note the [0] is for channel 0 which is high-gain channel\n", " image = event.dl1.tel[tel_id].image[0]\n", "\n", - " # cleaning the image\n", + " # Cleaning of the image\n", " cleaned_image = image\n", " # create a clean mask of pixels above the threshold\n", " cleanmask = tailcuts_clean(\n", - " camgeom, image, picture_thresh=10, boundary_thresh=5)\n", + " camgeom, image, picture_thresh=10, boundary_thresh=5\n", + " )\n", " # set all rejected pixels to zero\n", " cleaned_image[~cleanmask] = 0\n", "\n", - " # calulate hillas parameters\n", - " # it fails for empty pixels\n", + " # Calulate hillas parameters\n", + " # It fails for empty pixels\n", " try:\n", " hillas_params[tel_id] = hillas_parameters(camgeom, cleaned_image)\n", " except:\n", " pass\n", "\n", - " # not a stereo event -> next event\n", " if len(hillas_params) < 2:\n", " continue\n", "\n", - " # fit stereo event direction\n", - " reco.get_great_circles(\n", - " hillas_params, event.inst.subarray, point_azimuth, point_altitude)\n", - " \n", - " # reconstruct direction (3 components) with errors on the values\n", - " reco_direction, reco_dir_err = reco.fit_origin_crosses()\n", + " reco_result = reco.predict(hillas_params, event.inst, point_altitude, point_azimuth)\n", "\n", - " # in case fit reconstruct fails to get any real value -> next event\n", - " if np.isnan(reco_direction).any():\n", - " continue\n", + " # get angular offset between reconstructed shower direction and MC\n", + " # generated shower direction\n", + " off_angle = angular_separation(event.mc.az, event.mc.alt, reco_result.az, reco_result.alt)\n", "\n", - " # get angular offset between reconstructed event direction and MC shower direction\n", - " off_angle = linalg.angle(reco_direction, shower_direction)\n", - "\n", - " # appending all estimated off angles of each event\n", + " # Appending all estimated off angles\n", " off_angles.append(off_angle.to(u.deg).value)" ] }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 6, "metadata": { + "ExecuteTime": { + "end_time": "2018-06-15T12:49:47.507369Z", + "start_time": "2018-06-15T12:49:47.502642Z" + }, "collapsed": true }, "outputs": [], "source": [ - "# square the angles\n", - "thetasq = []\n", - "for i in off_angles:\n", - " thetasq.append(i**2)" + "# calculate theta square for angles which are not nan\n", + "off_angles = np.array(off_angles)\n", + "thetasquare = off_angles[np.isfinite(off_angles)]**2" ] }, { "cell_type": "code", - "execution_count": 11, - "metadata": {}, + "execution_count": 7, + "metadata": { + "ExecuteTime": { + "end_time": "2018-06-15T12:49:48.264122Z", + "start_time": "2018-06-15T12:49:47.511172Z" + } + }, "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEQCAYAAACjnUNyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFwFJREFUeJzt3X+wnmWd3/H3xxh0RVZXckA3cAysqbuhA1TPRlbpAF1l\nAxXRFisUceqPprrSrnbdFa3Fne50pq0zO5WKZjOaZe1U6LYSjEsQxDqLitgkDIYfApNGHJLahl+i\n+GMh+O0fz53t4+GcnCvJuc9zTs77NfNMnvu6rvt+vrcP5pP7x3PdqSokSZrJs0ZdgCRpYTAwJElN\nDAxJUhMDQ5LUxMCQJDUxMCRJTQwMSVITA0OS1MTAkCQ1MTAkSU2ePeoCZtOyZctqxYoVoy5DkhaM\nbdu2PVxVYy1jD6vAWLFiBVu3bh11GZK0YCT5XutYT0lJkpoYGJKkJgaGJKmJgSFJatJbYCQ5PslX\nk9yT5O4kvzfFmCS5IsmOJNuTvGKob02S+7q+y/qqU5LUps8jjL3A71fVKuA04L1JVk0acw6wsnut\nBT4FkGQJcGXXvwq4aIp1JUlzqLfAqKrvV9Xt3fsfAd8Blk8adj7w2Rq4DXhhkpcAq4EdVbWzqp4E\nrunGSpJGZE6uYSRZAfwd4FuTupYDDw4t7+rapmuXJI1I74GR5PnA54H3VdUPe9j+2iRbk2x96KGH\nZnvzh2zFZdePugRJmhW9BkaSpQzC4r9U1bVTDNkNHD+0fFzXNl37M1TV+qqaqKqJsbGmX7dLkg5C\nn3dJBfgM8J2q+pNphm0C3tbdLXUa8HhVfR/YAqxMckKSI4ALu7GSpBHpcy6p1wCXAHcmuaNr+zAw\nDlBV64DNwLnADuAnwNu7vr1JLgVuBJYAG6rq7h5rlSTNoLfAqKqvA5lhTAHvnaZvM4NAkSTNA/7S\nW5LUxMCQJDUxMCRJTQwMSVITA0OS1MTAkCQ1MTAkSU0MDElSEwNDktTEwJAkNTEwJElNDAxJUhMD\nQ5LUxMCQJDUxMCRJTQwMSVITA0OS1KS3J+4l2QC8HthTVX97iv4/AC4equM3gLGqejTJA8CPgKeB\nvVU10VedkqQ2fR5hXAWsma6zqj5WVadW1anAh4C/qqpHh4ac1fUbFpI0D/QWGFV1C/DojAMHLgKu\n7qsWSdKhG/k1jCTPY3Ak8vmh5gJuTrItydrRVCZJGtbbNYwDcB7wjUmno06vqt1JjgG+nOTe7ojl\nGbpAWQswPj7ef7WStEiN/AgDuJBJp6Oqanf35x5gI7B6upWran1VTVTVxNjYWK+FStJiNtLASPIC\n4AzgC0NtRyY5at974GzgrtFUKEnap8/baq8GzgSWJdkFfBRYClBV67phbwJuqqofD616LLAxyb76\nPldVX+qrTklSm94Co6ouahhzFYPbb4fbdgKn9FOVJOlgzYdrGJKkBcDAkCQ1MTAkSU0MDElSEwND\nktTEwJAkNTEwJElNDAxJUhMDQ5LUxMCQJDUxMCRJTQwMSVITA0OS1MTAkCQ1MTAkSU0MDElSEwND\nktSkt8BIsiHJniRTPo87yZlJHk9yR/e6fKhvTZL7kuxIcllfNUqS2vV5hHEVsGaGMV+rqlO7178B\nSLIEuBI4B1gFXJRkVY91SpIa9BYYVXUL8OhBrLoa2FFVO6vqSeAa4PxZLU6SdMBGfQ3j1Um2J7kh\nyUld23LgwaExu7o2SdIIPXuEn307MF5VTyQ5F7gOWHmgG0myFlgLMD4+PrsVSpL+xsiOMKrqh1X1\nRPd+M7A0yTJgN3D80NDjurbptrO+qiaqamJsbKzXmiVpMRtZYCR5cZJ071d3tTwCbAFWJjkhyRHA\nhcCmUdUpSRro7ZRUkquBM4FlSXYBHwWWAlTVOuAC4D1J9gI/BS6sqgL2JrkUuBFYAmyoqrv7qlOS\n1Ka3wKiqi2bo/wTwiWn6NgOb+6hLknRwRn2XlCRpgTAwJElNDAxJUhMDQ5LUxMCQJDUxMCRJTQwM\nSVITA0OS1MTAkCQ1MTAkSU0MDElSEwNDktTEwJAkNZkxMJIcmeRZ3fu/leQNSZb2X5okaT5pOcK4\nBXhukuXATcAlwFV9FiVJmn9aAiNV9RPgHwCfrKo3Ayf1W5Ykab5pCowkvwVcDFzftS3pryRJ0nzU\nEhi/B3wI2FhVdyc5EfjqTCsl2ZBkT5K7pum/OMn2JHcmuTXJKUN9D3TtdyTZ2rozkqT+tDyi9diq\nesO+harameRrDetdxeARrJ+dpv+7wBlV9ViSc4D1wKuG+s+qqocbPkeSNAdajjA+1Nj2C6rqFuDR\n/fTfWlWPdYu3Acc11CJJGpFpjzC6f/WfCyxPcsVQ1y8De2e5jncCNwwtF3BzkqeBP62q9bP8eZKk\nA7S/U1L/G9gKvAHYNtT+I+D9s1VAkrMYBMbpQ82nV9XuJMcAX05yb3fEMtX6a4G1AOPj47NVliRp\nkmkDo6q+DXw7yeeq6qk+PjzJycCngXOq6pGhz97d/bknyUZgNYPfg0xV53oG1z+YmJioPuqUJLVd\nw1id5MtJ7k+yM8l3k+w81A9OMg5cC1xSVfcPtR+Z5Kh974GzgSnvtJIkzZ2Wu6Q+w+AU1Dbg6dYN\nJ7kaOBNYlmQX8FFgKUBVrQMuB44GPpkEYG9VTQDHAhu7tmcDn6uqL7V+riSpHy2B8XhV3TDzsF9U\nVRfN0P8u4F1TtO8ETnnmGpKkUWoJjK8m+RiD00d/va+xqm7vrSpJ0rzTEhj7fkw3MdRWwN+b/XIk\nSfPVjIFRVWfNRSGSpPmt5XkYxyb5TJIbuuVVSd7Zf2mSpPmk5bbaq4AbgV/tlu8H3tdXQZKk+akl\nMJZV1V8APweoqr0cwO21kqTDQ0tg/DjJ0QwudJPkNODxXquSJM07LXdJ/T6wCfi1JN8AxoALeq1K\nkjTvtNwltS3JGcDLgQD39TW3lCRp/mq5S2o78IfAz6rqLsNCkhanlmsY5zF4/sVfJNmS5APdxIGS\npEVkxsCoqu9V1X+oqlcC/xg4mcHjVSVJi0jLRW+SvBR4S/d6msEpKknSIjJjYCT5FoNpyf8b8OZu\nNllJ0iLTcoTxtqq6r/dKJEnzWstF7x84l5QkybmkJElNeptLKsmGJHuSTPk87gxckWRHku1JXjHU\ntybJfV3fZY37IknqUZ9zSV0FrNlP/znAyu61FvhUt/0lwJVd/yrgoiSrGj5PktSjlove/5KDmEuq\nqm5JsmI/Q84HPltVBdyW5IVJXgKsAHbsuxsryTXd2HsaapUk9aTlh3u3A2cArwb+GXBSVW2fhc9e\nDjw4tLyra5uufUpJ1ibZmmTrQw89dNDFrLjs+infzzS2Zfxs1DQb4/ow25+90PelZRv7GzOK/R/l\n/+ZaWFpOSVFVe6vq7vk4l1RVra+qiaqaGBsbG3U5knTYavqld092A8cPLR/XtS2dpl2SNELTHmEk\neU3353N6+uxNwNu6u6VOAx6vqu8DW4CVSU5IcgRwYTdWkjRC+zvCuAJ4JfBN4BX7GTelJFcDZwLL\nkuwCPsrg6IGqWgdsBs4FdgA/Ad7e9e1NcimD334sATZU1d0H+vmSpNm1v8B4Ksl6YHmSKyZ3VtW/\n2N+Gq+qiGfoLeO80fZsZBIokaZ7YX2C8Hngt8DvAtrkpR5I0X00bGFX1MHBNku9U1bfnsCZJ0jzU\nclvtI0k2dtN87Eny+STH9V6ZJGleaQmMP2Nwl9Kvdq8vdm2SpEWkJTCOqao/6368t7eqrmIwPYgk\naRFpCYyHk7w1yZLu9Vbgkb4LkyTNLy2B8Q7gHwH/B/g+g4kH395nUZKk+WfGqUGq6nvAG+agFknS\nPNY0+aAkSQaGJKmJgSFJajJjYCT5yND7vmaulSTNc/ub3vyDSX6LX3wc6zf7L0mSNB/t7y6pe4E3\nAycm+Vq3fHSSl1fVfXNSnSRp3tjfKakfAB9m8LyKM4GPd+2XJbm157okSfPM/o4wfge4HPg14E+A\n7cCPq8of7UnSIjTtEUZVfbiqfht4APjPDJ5+N5bk60m+2LLxJGuS3JdkR5LLpuj/gyR3dK+7kjyd\n5EVd3wNJ7uz6th7U3kmSZs2Mv/QGbqyqrcDWJO+pqtOTLJtppSRLgCuB1wG7gC1JNlXVPfvGVNXH\ngI91488D3l9Vjw5t5qzuuRySpBGb8bbaqvrDocV/0rW1/CW+GthRVTur6kngGuD8/Yy/CLi6YbuS\npBE4oB/uHeCT95YDDw4t7+raniHJ84A1wOeHPw64Ocm2JGsPpE5J0uxrOSU1F84DvjHpdNTpVbU7\nyTHAl5PcW1W3TF6xC5O1AOPj43NTrSQtQn1ODbIbOH5o+biubSoXMul0VFXt7v7cA2xkcIrrGapq\nfVVNVNXE2JjPdZKkvvQZGFuAlUlOSHIEg1DYNHlQkhcAZwBfGGo7MslR+94DZwN39VirJGkGvZ2S\nqqq9SS4FbmRwS+6Gqro7ybu7/nXd0DcBN1XVj4dWPxbYmGRfjZ+rqi/1VaskaWa9XsOoqs3A5klt\n6yYtXwVcNaltJ3BKn7VJkg6M05tLkpoYGJKkJgaGJKmJgSFJamJgSJKaGBiSpCYGhiSpiYEhSWpi\nYEiSmhgYkqQmBoYkqYmBIUlqYmBIkpoYGJKkJgaGJKmJgSFJamJgSJKa9BoYSdYkuS/JjiSXTdF/\nZpLHk9zRvS5vXVeSNLd6e0RrkiXAlcDrgF3AliSbquqeSUO/VlWvP8h1JUlzpM8jjNXAjqraWVVP\nAtcA58/BupKkHvQZGMuBB4eWd3Vtk706yfYkNyQ56QDXlSTNkd5OSTW6HRivqieSnAtcB6w8kA0k\nWQusBRgfH5/9CiVJQL9HGLuB44eWj+va/kZV/bCqnujebwaWJlnWsu7QNtZX1URVTYyNjc1m/ZKk\nIX0GxhZgZZITkhwBXAhsGh6Q5MVJ0r1f3dXzSMu6kqS51dspqaram+RS4EZgCbChqu5O8u6ufx1w\nAfCeJHuBnwIXVlUBU67bV62SpJn1eg2jO820eVLbuqH3nwA+0bquJGl0/KW3JKmJgSFJamJgSJKa\nGBiSpCYGhiSpiYEhSWpiYEiSmhgYkqQmBoYkqYmBIUlqYmBIkpoYGJKkJgaGJKmJgSFJamJgSJKa\nGBiSpCYGhiSpSa+BkWRNkvuS7Ehy2RT9FyfZnuTOJLcmOWWo74Gu/Y4kW/usU5I0s94e0ZpkCXAl\n8DpgF7Alyaaqumdo2HeBM6rqsSTnAOuBVw31n1VVD/dVoySpXZ9HGKuBHVW1s6qeBK4Bzh8eUFW3\nVtVj3eJtwHE91iNJOgR9BsZy4MGh5V1d23TeCdwwtFzAzUm2JVnbQ32SpAPQ2ympA5HkLAaBcfpQ\n8+lVtTvJMcCXk9xbVbdMse5aYC3A+Pj4nNQrSYtRn0cYu4Hjh5aP69p+QZKTgU8D51fVI/vaq2p3\n9+ceYCODU1zPUFXrq2qiqibGxsZmsXxJ0rA+A2MLsDLJCUmOAC4ENg0PSDIOXAtcUlX3D7UfmeSo\nfe+Bs4G7eqxVkjSD3k5JVdXeJJcCNwJLgA1VdXeSd3f964DLgaOBTyYB2FtVE8CxwMau7dnA56rq\nS33VKkmaWa/XMKpqM7B5Utu6offvAt41xXo7gVMmt0uSRsdfekuSmhgYkqQmBoYkqYmBIUlqYmBI\nkpoYGJKkJgaGJKmJgSFJamJgSJKaGBiSpCYGhiSpiYEhSWpiYEiSmhgYkqQmBoYkqYmBIUlqYmBI\nkpr0GhhJ1iS5L8mOJJdN0Z8kV3T925O8onVdSdLc6i0wkiwBrgTOAVYBFyVZNWnYOcDK7rUW+NQB\nrCtJmkN9HmGsBnZU1c6qehK4Bjh/0pjzgc/WwG3AC5O8pHFdSdIc6jMwlgMPDi3v6tpaxrSsK0ma\nQ6mqfjacXACsqap3dcuXAK+qqkuHxvwl8O+q6uvd8leADwIrZlp3aBtrGZzOAng5cN9BlrwMePgg\n112o3OfD32LbX3CfD9RLq2qsZeCzD/IDWuwGjh9aPq5raxmztGFdAKpqPbD+UItNsrWqJg51OwuJ\n+3z4W2z7C+5zn/o8JbUFWJnkhCRHABcCmyaN2QS8rbtb6jTg8ar6fuO6kqQ51NsRRlXtTXIpcCOw\nBNhQVXcneXfXvw7YDJwL7AB+Arx9f+v2VaskaWZ9npKiqjYzCIXhtnVD7wt4b+u6PTvk01oLkPt8\n+Fts+wvuc296u+gtSTq8ODWIJKnJogqMQ5mqZKFq2OdfT/LNJH+d5AOjqHG2Nezzxd33e2eSW5Oc\nMoo6Z1PDPp/f7fMdSbYmOX0Udc6m1umDkvxmkr3drf4LWsP3fGaSx7vv+Y4kl89qAVW1KF4MLp7/\nL+BE4Ajg28CqSWPOBW4AApwGfGvUdc/BPh8D/Cbwb4EPjLrmOdrnVwO/0r0/Z5F8z8/n/5+CPhm4\nd9R1973PQ+P+B4ProReMuu45+J7PBP6yrxoW0xHGoUxVslDNuM9VtaeqtgBPjaLAHrTs861V9Vi3\neBuD3/ksZC37/ER1f6MARwIL/eJl6/RB/xz4PLBnLovrycinTFpMgXEoU5UsVIfb/rQ40H1+J4Oj\nyoWsaZ+TvCnJvcD1wDvmqLa+zLjPSZYDb6Kb1PQw0Prf9qu70483JDlpNgtYTIEh/YIkZzEIjA+O\nupa5UFUbq+rXgTcCfzzqeubAfwQ+WFU/H3Uhc+h2YLyqTgb+E3DdbG58MQXGoUxVslAdbvvTommf\nk5wMfBo4v6oemaPa+nJA33NV3QKcmGRZ34X1qGWfJ4BrkjwAXAB8Mskb56a8Xsy4z1X1w6p6onu/\nGVg6m9/zYgqMQ5mqZKFajFOszLjPScaBa4FLqur+EdQ421r2+WVJ0r1/BfAcYCEH5Yz7XFUnVNWK\nqloB/Hfgd6tqVv/FPcdavucXD33Pqxn8HT9r33Ovv/SeT+oQpipZqFr2OcmLga3ALwM/T/I+Bnde\n/HBkhR+Cxu/5cuBoBv/iBNhbC3iyusZ9/ocM/jH0FPBT4C1DF8EXnMZ9Pqw07vMFwHuS7GXwPV84\nm9+zv/SWJDVZTKekJEmHwMCQJDUxMCRJTQwMSVITA0OS1MTAkCQ1MTAkSU0WzQ/3pNnQTS3x9xn8\n0PEzVXXTiEuS5oxHGNI0kixJ8vEkd3cPWzqxqq6rqn8KvBt4yzTr/VKSv0qyZIq+PzrYB1UlOSLJ\nLUn8h55GwsCQpvchYGdVnQRcAfzuUN9HgCunWe8dwLVV9fRsFtM9A+ErTBNUUt8MDGkKSY4E3lRV\nH++avgu8rJuY8t8DN1TV7dOsfjHwhaFt/ask9yf5OvDySZ/z1iT/s3uc5p/uOypJ8q+7R3F+PcnV\nQ0cl13Xbl+ach7bS1F4LHJ/kjm75RcDNDJ7g9lrgBUleNnmSu24W0ROr6oFu+ZUMZhU9lcH/324H\ntnV9v8HgaOE1VfVUkk8CFyf5DoPJAk8Blg6vA9zF4JG60pwzMKSpnQpcvi8Qknwa2F5VVzA4PTWd\nZcAPhpb/LrCxqn7SbWd4OurfBl4JbOlmzf0lBo8SfRHwhar6GfCzJF/ct0JVPZ3kySRHVdWPDnUn\npQPhKSlpar/CYIp7uovMZwNf3O8aAz8Fntv4GQH+vKpO7V4vr6o/aljvOcDPGj9DmjUGhjS1+4HT\nuvfvB66vqu/OtFJVPQYsSbIvNG4B3tjdOXUUcN7Q8K8AFyQ5BiDJi5K8FPgGcF6S5yZ5PvD6fSsk\nORp4uKqeOsT9kw6Yp6SkqV0N3JBkB/BNYO0BrHsTcDpwc1XdnuS/At9mcLppy75BVXVPko8ANyV5\nFvAU8N6quq07dbUd+L/AncDj3WpnAdcf2q5JB8cHKEmzrHsE6vur6pJD2Mbzq+qJJM9jcJSytguf\na4HLDpNHy2qB8QhDmmXdX+xfTbLkEH6LsT7JKgbXQ/682+YRwHWGhUbFIwxJUhMvekuSmhgYkqQm\nBoYkqYmBIUlqYmBIkpoYGJKkJgaGJKmJgSFJamJgSJKa/D8uWuK1AhBpuQAAAABJRU5ErkJggg==\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEQCAYAAABfiGi4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvFvnyVgAAE99JREFUeJzt3X2QZXV95/H3xwHUBKIJMyZkmGFC\nRCvEjQizLJTZjQaTICpkS0xwRYNhM7WsbtToZsG46Jp/NKk1JasJmQ0EtBLEKJJhgUJj2EVcIcxM\nYOQhULOIYZBdHtRRgihDvvvHPfOrtjPdfbqnz7109/tVdWvOw++e+/31w3z6PP1OqgpJkgCeMekC\nJElPH4aCJKkxFCRJjaEgSWoMBUlSYyhIkhpDQZLUGAqSpMZQkCQ1hoIkqTlg0gXM1+rVq2vDhg2T\nLkOSlpRt27Y9UlVr5mq35EJhw4YNbN26ddJlSNKSkuSrfdp5+EiS1BgKkqTGUJAkNYaCJKkZLBSS\nPCvJ3yS5LckdSf7LPto8M8nlSXYmuTnJhqHqkSTNbcg9he8CP19VLwaOAU5OcsK0NmcD36iq5wN/\nAHxwwHokSXMYLBRq5LFu9sDuNf3Zn6cBl3bTnwJOSpKhapIkzW7QcwpJViW5FXgI+FxV3TytyVrg\nfoCq2gPsBg4dsiZJ0swGvXmtqp4CjknyXOAzSV5UVbfPdztJNgGbANavX7/gejace/W833PfB161\n4M+TpKVmLFcfVdU3geuBk6etegBYB5DkAOA5wKP7eP/mqtpYVRvXrJnzLm1J0gINefXRmm4PgSTP\nBn4B+LtpzbYAv9ZNnw78dVVNP+8gSRqTIQ8fHQZcmmQVo/D5ZFX9jyTvB7ZW1RbgIuDjSXYCXwfO\nGLAeSdIcBguFqtoBvGQfy8+fMv0E8LqhapAkzY93NEuSGkNBktQYCpKkxlCQJDWGgiSpMRQkSY2h\nIElqDAVJUmMoSJIaQ0GS1BgKkqTGUJAkNYaCJKkxFCRJjaEgSWoMBUlSYyhIkhpDQZLUGAqSpMZQ\nkCQ1hoIkqTEUJEmNoSBJagwFSVJjKEiSGkNBktQMFgpJ1iW5PsmdSe5I8rZ9tHlZkt1Jbu1e5w9V\njyRpbgcMuO09wDuranuSQ4BtST5XVXdOa/eFqnr1gHVIknoabE+hqh6squ3d9LeBu4C1Q32eJGn/\njeWcQpINwEuAm/ex+sQktyW5NslPz/D+TUm2Jtn68MMPD1ipJK1sg4dCkoOBTwNvr6pvTVu9HTii\nql4M/Dfgyn1to6o2V9XGqtq4Zs2aYQuWpBVs0FBIciCjQPizqrpi+vqq+lZVPdZNXwMcmGT1kDVJ\nkmY25NVHAS4C7qqqD83Q5se6diQ5vqvn0aFqkiTNbsirj14KvBH4cpJbu2XvBtYDVNWFwOnAOUn2\nAN8BzqiqGrAmSdIsBguFqroRyBxtPgJ8ZKgaJEnz4x3NkqTGUJAkNYaCJKkxFCRJjaEgSWoMBUlS\nYyhIkhpDQZLUGAqSpMZQkCQ1hoIkqTEUJEmNoSBJagwFSVJjKEiSGkNBktQYCpKkxlCQJDWGgiSp\nMRQkSY2hIElqDAVJUmMoSJIaQ0GS1BgKkqRmsFBIsi7J9UnuTHJHkrfto02SXJBkZ5IdSY4dqh5J\n0twOGHDbe4B3VtX2JIcA25J8rqrunNLmlcBR3etfAH/U/StJmoDB9hSq6sGq2t5Nfxu4C1g7rdlp\nwMdq5CbguUkOG6omSdLsxnJOIckG4CXAzdNWrQXunzK/i38aHJKkMRny8BEASQ4GPg28vaq+tcBt\nbAI2Aaxfv34Rq1u4Dedevc/l933gVWOuRJIWz5x7Ckl+MMkzuukXJDk1yYF9Nt61+zTwZ1V1xT6a\nPACsmzJ/eLfs+1TV5qraWFUb16xZ0+ejJUkL0Ofw0Q3As5KsBT4LvBG4ZK43JQlwEXBXVX1ohmZb\ngDd1VyGdAOyuqgd7VS5JWnR9Dh+lqh5Pcjbwh1X1e0lu7fG+lzIKkC9Paf9uYD1AVV0IXAOcAuwE\nHgfePN8OSJIWT69QSHIi8Abg7G7ZqrneVFU3ApmjTQFv6VGDJGkM+hw+ehtwHvCZqrojyZHA9cOW\nJUmahD57Cj9aVafunamqe5N8YcCaJEkT0mdP4byeyyRJS9yMewpJXsnoJPDaJBdMWfVDjIawkCQt\nM7MdPvoasBU4Fdg2Zfm3gXcMWZQkaTJmDIWqug24LcmfV9WTY6xJkjQhfU40H5/kfcARXfswupr0\nyCELkySNX59QuIjR4aJtwFPDliNJmqQ+obC7qq4dvBJJ0sT1CYXrk/w+cAXw3b0L9z4rQZK0fPQJ\nhb1PQts4ZVkBP7/45UiSJmnOUKiql4+jEEnS5PV5nsKPJrkoybXd/NHdiKmSpGWmzzAXlwDXAT/e\nzd8DvH2ogiRJk9MnFFZX1SeBfwSoqj14aaokLUt9QuEfkhzK6OQye5+QNmhVkqSJ6HP10TsZPTbz\nJ5N8EVgDnD5oVZKkiehz9dG2JD8HvJDREBd3OxaSJC1Pfa4+2gH8NvBEVd1uIEjS8tXnnMJrGD0/\n4ZNJbknyriTrB65LkjQBc4ZCVX21qn6vqo4D/g3wM8BXBq9MkjR2fU40k+QI4Fe711OMDidJkpaZ\nOUMhyc3AgcBfAK+rqnsHr0qSNBF99hTeVFV3D16JJGni+pxo/qZjH0nSyuDYR5KkZrCxj5JcnOSh\nJLfPsP5lSXYnubV7nT+vyiVJi67POYWFjn10CfAR4GOztPlCVb26x7YkSWPQJxR+iwWMfVRVNyTZ\nsF/VSZLGqs/YR9sHHPvoxCS3AV8D3lVVd+yrUZJNwCaA9eu9mVqShtLr5rXuPMI+/8PeD9uBI6rq\nsSSnAFcCR83w+ZuBzQAbN26sRa5DktTpc6J5EFX1rap6rJu+BjgwyepJ1SNJmiUUkry0+/eZQ3xw\nkh9Lkm76+K6WR4f4LElSP7MdProAOA74EnDsfDec5DLgZcDqJLuA9zIaLoOqupDRyepzkuwBvgOc\nUVUeGpKkCZotFJ5MshlYm+SC6Sur6jdn23BVvX6O9R9hdMmqJOlpYrZQeDXwCuCXgG3jKUeSNEkz\nhkJVPQJ8IsldVXXbGGuSJE1In6uPHk3ymW7IioeSfDrJ4YNXJkkauz6h8KeM7mj+8e51VbdMkrTM\n9AmF51XVn1bVnu51CaOhLiRJy0yfUHgkyZlJVnWvM/F+AklalvqEwq8DvwL8X+BBRvcXvHnIoiRJ\nk9FnQLyvAqeOoRZJ0oRNbOwjSdLTj6EgSWoMBUlSM2coJHnPlOlBRkyVJD09zDZ09n9KciLf/+jN\nLw1fkiRpUma7+ujvgNcBRyb5Qjd/aJIXVtXdY6lOkjRWsx0++ibwbmAno+cifLhbfm6S/z1wXZKk\nCZhtT+GXgPOBnwQ+BOwA/qGqvHFNkpapGfcUqurdVXUScB/wcWAVsCbJjUmuGlN9kqQxmvOOZuC6\nqtoKbE1yTlX9bJLVQxcmSRq/OS9JrarfnjJ7VrfskaEKkiRNzrxuXvMJbJK0vHlHsySpMRQkSY2h\nIElqDAVJUmMoSJKawUIhycVJHkpy+wzrk+SCJDuT7Ehy7FC1SJL6GXJP4RLg5FnWvxI4qnttAv5o\nwFokST0MFgpVdQPw9VmanAZ8rEZuAp6b5LCh6pEkzW2S5xTWAvdPmd/VLZMkTUifsY8mLskmRoeY\nWL9+/YSrWXwbzr16Xu3v+8CrBv+M+Zqtppk+eyH90PdbqV/bcfzOrFST3FN4AFg3Zf7wbtk/UVWb\nq2pjVW1cs2bNWIqTpJVokqGwBXhTdxXSCcDuqnpwgvVI0oo32OGjJJcxemLb6iS7gPcCBwJU1YXA\nNcApjJ7s9jjgw3skacIGC4Wqev0c6wt4y1CfL0maP+9oliQ1hoIkqTEUJEmNoSBJagwFSVJjKEiS\nGkNBktQYCpKkxlCQJDWGgiSpMRQkSY2hIElqDAVJUmMoSJIaQ0GS1BgKkqTGUJAkNYaCJKkxFCRJ\njaEgSWoMBUlSYyhIkhpDQZLUGAqSpMZQkCQ1hoIkqRk0FJKcnOTuJDuTnLuP9WcleTjJrd3r3w5Z\njyRpdgcMteEkq4CPAr8A7AJuSbKlqu6c1vTyqnrrUHVIkvobck/heGBnVd1bVd8DPgGcNuDnSZL2\n05ChsBa4f8r8rm7ZdK9NsiPJp5Ks29eGkmxKsjXJ1ocffniIWiVJTP5E81XAhqr6GeBzwKX7alRV\nm6tqY1VtXLNmzVgLlKSVZMhQeACY+pf/4d2ypqoerarvdrN/Ahw3YD2SpDkMGQq3AEcl+YkkBwFn\nAFumNkhy2JTZU4G7BqxHkjSHwa4+qqo9Sd4KXAesAi6uqjuSvB/YWlVbgN9MciqwB/g6cNZQ9UiS\n5jZYKABU1TXANdOWnT9l+jzgvCFrkCT1N+kTzZKkpxFDQZLUGAqSpMZQkCQ1hoIkqTEUJEmNoSBJ\nagwFSVJjKEiSGkNBktQYCpKkxlCQJDWGgiSpMRQkSY2hIElqDAVJUmMoSJIaQ0GS1BgKkqTGUJAk\nNYaCJKkxFCRJjaEgSWoMBUlSYyhIkppBQyHJyUnuTrIzybn7WP/MJJd3629OsmHIeiRJsxssFJKs\nAj4KvBI4Gnh9kqOnNTsb+EZVPR/4A+CDQ9UjSZrbkHsKxwM7q+reqvoe8AngtGltTgMu7aY/BZyU\nJAPWJEmaxZChsBa4f8r8rm7ZPttU1R5gN3DogDVJkmZxwKQL6CPJJmBTN/tYkrsXuKnVwCPz+ux5\nHtCab/uFmOdnzLvPC7GQfg/4tRpLn59mvq/P4/g5fBro/X1eRl+P/fnZPqJPoyFD4QFg3ZT5w7tl\n+2qzK8kBwHOAR6dvqKo2A5v3t6AkW6tq4/5uZymxzyuDfV4ZxtHnIQ8f3QIcleQnkhwEnAFsmdZm\nC/Br3fTpwF9XVQ1YkyRpFoPtKVTVniRvBa4DVgEXV9UdSd4PbK2qLcBFwMeT7AS+zig4JEkTMug5\nhaq6Brhm2rLzp0w/AbxuyBqm2e9DUEuQfV4Z7PPKMHif49EaSdJeDnMhSWqWZSisxOE1evT5t5Lc\nmWRHks8n6XV52tPZXH2e0u61SSrJkr9SpU+fk/xK972+I8mfj7vGxdbjZ3t9kuuT/G33833KJOpc\nLEkuTvJQkttnWJ8kF3Rfjx1Jjl3UAqpqWb0YndT+P8CRwEHAbcDR09r8e+DCbvoM4PJJ1z2GPr8c\n+IFu+pyV0Oeu3SHADcBNwMZJ1z2G7/NRwN8CP9zNP2/SdY+hz5uBc7rpo4H7Jl33fvb5XwHHArfP\nsP4U4FogwAnAzYv5+ctxT2ElDq8xZ5+r6vqqerybvYnRfSNLWZ/vM8DvMhpT64lxFjeQPn3+DeCj\nVfUNgKp6aMw1LrY+fS7gh7rp5wBfG2N9i66qbmB0NeZMTgM+ViM3Ac9Ncthiff5yDIWVOLxGnz5P\ndTajvzSWsjn73O1Wr6uqq8dZ2ID6fJ9fALwgyReT3JTk5LFVN4w+fX4fcGaSXYyudvwP4yltYub7\n+z4vS2KYCy2eJGcCG4Gfm3QtQ0ryDOBDwFkTLmXcDmB0COlljPYGb0jyz6rqmxOtalivBy6pqv+a\n5ERG9z69qKr+cdKFLUXLcU9hPsNrMNvwGktInz6T5BXA7wCnVtV3x1TbUObq8yHAi4D/meQ+Rsde\ntyzxk819vs+7gC1V9WRVfQW4h1FILFV9+nw28EmAqvoS8CxGYwQtV71+3xdqOYbCShxeY84+J3kJ\n8MeMAmGpH2eGOfpcVburanVVbaiqDYzOo5xaVVsnU+6i6POzfSWjvQSSrGZ0OOnecRa5yPr0+e+B\nkwCS/BSjUHh4rFWO1xbgTd1VSCcAu6vqwcXa+LI7fFQrcHiNnn3+feBg4C+6c+p/X1WnTqzo/dSz\nz8tKzz5fB/xikjuBp4D/WFVLdi+4Z5/fCfz3JO9gdNL5rKX8R16SyxgF++ruPMl7gQMBqupCRudN\nTgF2Ao8Db17Uz1/CXztJ0iJbjoePJEkLZChIkhpDQZLUGAqSpMZQkCQ1hoIkqTEUJEnNsrt5Tdpf\nSX4ZeBWjkTcvqqrPTrgkaWzcU9CKlmRVkg93D6T5cpIjq+rKqvoN4N8BvzrD+56d5H8lWbWPde9L\n8q4F1nNQkhu6MbmksTMUtNKdB9xbVT8NXMDoAUx7vQf46Azv+3Xgiqp6ajGL6Z4Z8HlmCCNpaIaC\nVqwkPwj866r6cLfoK8Dzu4HGPghcW1XbZ3j7G4C/nLKt30lyT5IbgRdO+5wzk/xNkluT/PHevYsk\n/7l7zOSNSS6bsndxZbd9aezcRdVK9gpgXZJbu/kfAf6K0UNaXgE8J8nzu0HImm60ziOr6r5u/jhG\ngyoew+h3ajuwrVv3U4z+6n9pVT2Z5A+BNyS5C3gt8GJGg5219wC3A/98kB5LczAUtJIdA5y/9z/9\nJH8C7KiqCxgdSprJamDqQ2v+JfCZvY87TTJ1hNaTgOOAW7rRaZ8NPMQogP6yqp4Ankhy1d43VNVT\nSb6X5JCq+vb+dlKaDw8faSX7YUZDD+992NIvAlfN+o6R7zAas7+PAJdW1THd64VV9b4e73smy+O5\n0lpiDAWtZPcweiIbwDuAq7unlc2qqr4BrEqyNxhuAH65uyLpEOA1U5p/Hjg9yfMAkvxIkiOALwKv\nSfKsJAcDr977hiSHAo9U1ZP72T9p3jx8pJXsMuDa7mFLXwI2zeO9nwV+Fvirqtqe5HLgNkaHhm7Z\n26iq7kzyHuCz3XOjnwTeUlU3dYeZdgD/D/gysLt728uBq/eva9LC+JAdaQGSHAu8o6reuB/bOLiq\nHkvyA4z2NjZ1AXMFcG5V3bNY9Up9uacgLUD3n/f1SVbtx70Km5Mczej8xKXdNg8CrjQQNCnuKUiS\nGk80S5IaQ0GS1BgKkqTGUJAkNYaCJKkxFCRJjaEgSWoMBUlSYyhIkpr/D3oc7amjIOUyAAAAAElF\nTkSuQmCC\n", "text/plain": [ - "" + "
" ] }, "metadata": {}, @@ -205,7 +195,7 @@ ], "source": [ "# plot \n", - "plt.hist(thetasq, bins=np.linspace(0, 0.5, 300))\n", + "plt.hist(thetasquare, bins=np.linspace(0, 1, 50))\n", "plt.xlabel(r'$\\theta^2$ (deg)')\n", "plt.ylabel(\"# of events\")\n", "plt.show()" @@ -237,7 +227,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.2" + "version": "3.6.5" }, "toc": { "nav_menu": { diff --git a/examples/plot_theta_square.py b/examples/plot_theta_square.py index 1ecd98e857f..a3216a651f0 100644 --- a/examples/plot_theta_square.py +++ b/examples/plot_theta_square.py @@ -1,52 +1,52 @@ -""" Create a theta-square plot . - """ +Create a theta-square plot . +""" + +import sys import matplotlib.pyplot as plt import numpy as np from astropy import units as u +from astropy.coordinates.angle_utilities import angular_separation from ctapipe.calib import CameraCalibrator from ctapipe.image import hillas_parameters from ctapipe.image import tailcuts_clean from ctapipe.io import event_source from ctapipe.reco import HillasReconstructor -from ctapipe.utils import datasets, linalg +from ctapipe.utils import datasets + + +if len(sys.argv) >= 2: + filename = sys.argv[1] +else: + # importing data from avaiable datasets in ctapipe + filename = datasets.get_dataset_path("gamma_test_large.simtel.gz") -# importing data from avaiable datasets in ctapipe -filename = datasets.get_dataset_path("gamma_test_large.simtel.gz") # reading the Monte Carlo file for LST source = event_source(filename, allowed_tels={1, 2, 3, 4}) -# pointing direction of the telescopes -point_azimuth = {} -point_altitude = {} - reco = HillasReconstructor() calib = CameraCalibrator(r1_product="HESSIOR1Calibrator") off_angles = [] for event in source: - # The direction the incident particle. Converting Monte Carlo Shower - # parameter theta and phi to corresponding to 3 components (x,y,z) of a - # vector - shower_azimuth = event.mc.az # same as in Monte Carlo file i.e. phi - shower_altitude = np.pi * u.rad / 2 - event.mc.alt # altitude = 90 - theta - shower_direction = linalg.set_phi_theta(shower_azimuth, shower_altitude) # calibrating the event calib.calibrate(event) hillas_params = {} subarray = event.inst.subarray + # pointing direction of the telescopes + point_azimuth = {} + point_altitude = {} + for tel_id in event.dl0.tels_with_data: # telescope pointing direction point_azimuth[tel_id] = event.mc.tel[tel_id].azimuth_raw * u.rad - point_altitude[tel_id] = ( - np.pi / 2 - event.mc.tel[tel_id].altitude_raw - ) * u.rad + point_altitude[tel_id] = event.mc.tel[tel_id].altitude_raw * u.rad # print(point_azimuth,point_altitude) # Camera Geometry required for hillas parametrization @@ -74,36 +74,30 @@ if len(hillas_params) < 2: continue - reco.get_great_circles( - hillas_params, event.inst.subarray, point_azimuth, point_altitude - ) - - # fit the gamma's direction of origin - # return reconstructed direction (3 components) with errors on the values - reco_direction, reco_dir_err = reco.fit_origin_crosses() - - # In case fit fails to get any real value - if np.isnan(reco_direction).any(): - continue + reco_result = reco.predict(hillas_params, event.inst, point_altitude, point_azimuth) # get angular offset between reconstructed shower direction and MC # generated shower direction - off_angle = linalg.angle(reco_direction, shower_direction) + off_angle = angular_separation( + event.mc.az, + event.mc.alt, + reco_result.az, + reco_result.alt + ) # Appending all estimated off angles off_angles.append(off_angle.to(u.deg).value) -# calculate theta square -thetasq = [] -for i in off_angles: - thetasq.append(i**2) +# calculate theta square for angles which are not nan +off_angles = np.array(off_angles) +thetasquare = off_angles[np.isfinite(off_angles)]**2 # To plot thetasquare The number of events in th data files for LSTCam is not # significantly high to give a nice thetasquare plot for gammas One can use # deedicated MC file for LST get nice plot plt.figure(figsize=(10, 8)) -plt.hist(thetasq, bins=np.linspace(0, 10, 50)) +plt.hist(thetasquare, bins=np.linspace(0, 1, 50)) plt.title(r'$\theta^2$ plot') -plt.xlabel(r'$\theta^2$') -plt.ylabel("# of events") +plt.xlabel(r'$\theta^2$ (deg)') +plt.ylabel('# of events') plt.show()