Note
Go to the end to download the full example code.
Event Display with Hillas Parameters#
This script displays events with Hillas parameter reconstruction and visualization of all parameters on the camera image.
8 import numpy as np
9 import matplotlib.pyplot as plt
10 from matplotlib.patches import Ellipse, Arc
11 from matplotlib.lines import Line2D
12
13 from ctapipe.instrument import SubarrayDescription
14
15 from ctapipe.coordinates import TelescopeFrame
16 from ctapipe.image.toymodel import Gaussian
17 from ctapipe.visualization import CameraDisplay
18 from ctapipe.image.cleaning import tailcuts_clean
19 from ctapipe.image.hillas import hillas_parameters, HillasParameterizationError
20
21 import astropy.units as u
Camera Display with Hillas Parameters annotated#
28 def display_event_with_annotated_hillas(
29 image, geom, picture_thresh=10, boundary_thresh=5
30 ):
31 """
32 Display an event with detailed annotations showing what each Hillas parameter represents.
33 """
34
35 # Clean the image
36 mask = tailcuts_clean(
37 geom, image, picture_thresh=picture_thresh, boundary_thresh=boundary_thresh
38 )
39
40 # Calculate Hillas parameters
41 try:
42 hillas = hillas_parameters(geom[mask], image[mask])
43 except HillasParameterizationError:
44 print("Could not parametrize event")
45 return None
46
47 fig, ax = plt.subplots(figsize=(12, 10))
48
49 # Display the camera image
50 display = CameraDisplay(geom, ax=ax, cmap="gray")
51 display.image = image * mask
52 # display.highlight_pixels(mask, color='red', linewidth=2, alpha=0.4)
53 display.add_colorbar(ax=ax, label="Charge [p.e.]")
54
55 # Define colors for different elements using a colorblind-friendly palette
56 cog_color = "red"
57 ellipse_color = "#D55E00" # vermilion
58 length_color = "#0072B2" # blue
59 width_color = "#009E73" # green
60 angle_color = "#CC79A7" # pink
61 radial_color = "#E69F00" # orange
62
63 fontsize = 12
64
65 x = hillas.fov_lon
66 y = hillas.fov_lat
67
68 # 1. Center of Gravity (x, y)
69 ax.plot(
70 x.value,
71 y.value,
72 "o",
73 color=cog_color,
74 markersize=15,
75 markeredgewidth=3,
76 markerfacecolor="none",
77 label="COG (x, y)",
78 )
79
80 ax.plot(
81 x.value,
82 y.value,
83 "o",
84 color=cog_color,
85 markersize=15,
86 markeredgewidth=3,
87 markerfacecolor="none",
88 label="COG (FOV lon, lat)",
89 )
90
91 # 2. Hillas Ellipse (shows length and width)
92 # Note: Ellipse angle is rotation of the width (horizontal) axis
93 # Since height=length (major axis) should be along psi, we add 90 degrees
94 ellipse = Ellipse(
95 xy=(x.value, y.value),
96 width=hillas.width.value * 2,
97 height=hillas.length.value * 2,
98 angle=np.degrees(hillas.psi.value) + 90,
99 fill=False,
100 color=ellipse_color,
101 linewidth=3,
102 linestyle="--",
103 label="Hillas Ellipse",
104 )
105 ax.add_patch(ellipse)
106
107 # 3. Length axis (major axis)
108 length_end_x = x.value + hillas.length.value * np.cos(hillas.psi.value)
109 length_end_y = y.value + hillas.length.value * np.sin(hillas.psi.value)
110 length_start_x = x.value
111 length_start_y = y.value
112
113 ax.plot(
114 [length_start_x, length_end_x],
115 [length_start_y, length_end_y],
116 color=length_color,
117 linewidth=3,
118 label="Length",
119 zorder=10,
120 )
121
122 long_axis_end_x = x.value + 2 * hillas.length.value * np.cos(hillas.psi.value)
123 long_axis_end_y = y.value + 2 * hillas.length.value * np.sin(hillas.psi.value)
124 long_axis_start_x = x.value - 2 * hillas.length.value * np.cos(hillas.psi.value)
125 long_axis_start_y = y.value - 2 * hillas.length.value * np.sin(hillas.psi.value)
126
127 ax.plot(
128 [long_axis_start_x, long_axis_end_x],
129 [long_axis_start_y, long_axis_end_y],
130 color=length_color,
131 linewidth=3,
132 label="long-axis",
133 zorder=10,
134 alpha=0.5,
135 ls="--",
136 )
137
138 # Annotate length
139 mid_length_x = x.value + 0.7 * hillas.length.value * np.cos(hillas.psi.value)
140 mid_length_y = y.value + 0.7 * hillas.length.value * np.sin(hillas.psi.value)
141 ax.annotate(
142 f"Length\n{hillas.length:.3f}",
143 xy=(mid_length_x, mid_length_y),
144 xytext=(mid_length_x + 0.3, mid_length_y - 0.2),
145 color=length_color,
146 fontsize=fontsize,
147 fontweight="bold",
148 arrowprops=dict(arrowstyle="->", color=length_color, lw=2),
149 bbox=dict(boxstyle="round,pad=0.5", facecolor="black", alpha=1),
150 )
151
152 # 4. Width axis (minor axis)
153 width_end_x = x.value + hillas.width.value * np.sin(hillas.psi.value)
154 width_end_y = y.value - hillas.width.value * np.cos(hillas.psi.value)
155 width_start_x = x.value
156 width_start_y = y.value
157
158 ax.plot(
159 [width_start_x, width_end_x],
160 [width_start_y, width_end_y],
161 color=width_color,
162 linewidth=3,
163 label="Width",
164 zorder=10,
165 )
166
167 # Annotate width
168 mid_width_x = x.value + 0.7 * hillas.width.value * np.sin(hillas.psi.value)
169 mid_width_y = y.value - 0.7 * hillas.width.value * np.cos(hillas.psi.value)
170 ax.annotate(
171 f"Width\n{hillas.width:.3f}",
172 xy=(mid_width_x, mid_width_y),
173 xytext=(mid_width_x - 0.7, mid_width_y + 0.5),
174 color=width_color,
175 fontsize=fontsize,
176 fontweight="bold",
177 arrowprops=dict(arrowstyle="->", color=width_color, lw=2),
178 bbox=dict(boxstyle="round,pad=0.5", facecolor="black", alpha=1),
179 )
180
181 # 5. Angle psi (orientation of major axis)
182 # Draw the angle between the length axis and the x-axis (horizontal)
183 # We'll draw this at the end of the length axis
184
185 # Draw a horizontal reference line at the end of the length axis
186 psi_ref_length = 0.25
187 ref_line_x_start = length_end_x
188 ref_line_x_end = length_end_x + psi_ref_length
189 ref_line_y = length_end_y
190
191 ax.plot(
192 [ref_line_x_start, ref_line_x_end],
193 [ref_line_y, ref_line_y],
194 color=angle_color,
195 linewidth=3,
196 linestyle="--",
197 alpha=1,
198 zorder=10,
199 )
200
201 # Draw arc showing the psi angle at the end of the length axis
202 arc_radius = 0.25
203 # The arc should go from the horizontal (0°) to the length axis direction (psi)
204 angle_to_cog = np.degrees(hillas.psi.value)
205
206 # Normalize to 0-360
207 while angle_to_cog < 0:
208 angle_to_cog += 360
209 while angle_to_cog >= 360:
210 angle_to_cog -= 360
211
212 # Draw arc from horizontal (0°) to the direction pointing back to COG
213 arc = Arc(
214 (length_end_x, length_end_y),
215 2 * arc_radius,
216 2 * arc_radius,
217 angle=0,
218 theta1=0,
219 theta2=angle_to_cog,
220 color=angle_color,
221 linewidth=2.5,
222 zorder=10,
223 )
224 ax.add_patch(arc)
225
226 # Annotate psi angle - place it along the bisector of the arc
227 psi_bisector = angle_to_cog / 2
228 psi_label_x = length_end_x + arc_radius * 1.4 * np.cos(np.radians(psi_bisector))
229 psi_label_y = length_end_y + arc_radius * 1.4 * np.sin(np.radians(psi_bisector))
230 ax.annotate(
231 f"ψ = {np.degrees(hillas.psi.value):.1f}°",
232 xy=(psi_label_x, psi_label_y),
233 color=angle_color,
234 fontsize=fontsize,
235 fontweight="bold",
236 bbox=dict(boxstyle="round,pad=0.5", facecolor="black", alpha=0.7),
237 )
238
239 # 6. Radial distance r and angle phi from camera center
240 camera_center_x = 0
241 camera_center_y = 0
242
243 # Draw line from camera center to COG
244 ax.plot(
245 [camera_center_x, x.value],
246 [camera_center_y, y.value],
247 color=radial_color,
248 linewidth=2.5,
249 linestyle=":",
250 label="r (radial)",
251 zorder=5,
252 )
253
254 # Mark camera center
255 ax.plot(
256 camera_center_x,
257 camera_center_y,
258 "x",
259 color=radial_color,
260 markersize=20,
261 markeredgewidth=3,
262 )
263
264 # Annotate camera center
265 ax.annotate(
266 "Camera\nCenter",
267 xy=(camera_center_x, camera_center_y),
268 xytext=(-0.25, -0.15),
269 color=radial_color,
270 fontsize=10,
271 fontweight="bold",
272 bbox=dict(boxstyle="round,pad=0.5", facecolor="black", alpha=0.7),
273 )
274
275 # Annotate r
276 mid_r_x = x.value / 2
277 mid_r_y = y.value / 2
278 ax.annotate(
279 f"r = {hillas.r:.3f}",
280 xy=(mid_r_x, mid_r_y),
281 xytext=(mid_r_x + 0.2, mid_r_y - 0),
282 color=radial_color,
283 fontsize=fontsize,
284 fontweight="bold",
285 arrowprops=dict(arrowstyle="->", color=radial_color, lw=2),
286 bbox=dict(boxstyle="round,pad=0.5", facecolor="black", alpha=0.7),
287 )
288
289 # Draw arc for phi angle
290 phi_arc_radius = 0.25
291 phi_arc = Arc(
292 (camera_center_x, camera_center_y),
293 2 * phi_arc_radius,
294 2 * phi_arc_radius,
295 angle=0,
296 theta1=0,
297 theta2=np.degrees(hillas.phi.value),
298 color=radial_color,
299 linewidth=2,
300 linestyle="--",
301 zorder=5,
302 )
303 ax.add_patch(phi_arc)
304
305 # Annotate phi
306 phi_label_angle = hillas.phi.value / 2
307 phi_label_x = phi_arc_radius * 1.8 * np.cos(phi_label_angle)
308 phi_label_y = phi_arc_radius * 1.8 * np.sin(phi_label_angle)
309 ax.annotate(
310 f"φ = {np.degrees(hillas.phi.value):.1f}°",
311 xy=(phi_label_x, phi_label_y),
312 color=radial_color,
313 fontsize=fontsize,
314 fontweight="bold",
315 bbox=dict(boxstyle="round,pad=0.5", facecolor="black", alpha=0.7),
316 )
317
318 # Add reference x-axis line for angle measurement
319 ax.plot(
320 [-2.2, 2.2],
321 [camera_center_y, camera_center_y],
322 "w--",
323 linewidth=2,
324 alpha=0.8,
325 zorder=1,
326 )
327
328 # Add legend with all Hillas parameters
329 param_text = (
330 f"Hillas Parameters\n"
331 f"━━━━━━━━━━━━━━━━━━━━\n"
332 f"Intensity: {hillas.intensity:.1f} p.e.\n"
333 f"\n"
334 f"fov_long: {x:.4f}\n"
335 f"fov_lat: {y:.4f}\n"
336 f"r: {hillas.r:.4f}\n"
337 f"φ: {np.degrees(hillas.phi.value):.2f}°\n"
338 f"\n"
339 f"Length: {hillas.length:.4f}\n"
340 f"Width: {hillas.width:.4f}\n"
341 f"ψ: {np.degrees(hillas.psi.value):.2f}°\n"
342 f"\n"
343 f"Skewness: {hillas.skewness:.3f}\n"
344 f"Kurtosis: {hillas.kurtosis:.3f}"
345 )
346
347 ax.text(
348 0.02,
349 0.98,
350 param_text,
351 transform=ax.transAxes,
352 fontsize=10,
353 verticalalignment="top",
354 bbox=dict(boxstyle="round", facecolor="white", alpha=0.9),
355 family="monospace",
356 )
357
358 ax.set_title("Annotated Hillas Parameters", fontsize=14, pad=20, fontweight="bold")
359
360 # Custom legend
361 legend_elements = [
362 Line2D(
363 [0],
364 [0],
365 marker="o",
366 color="w",
367 markerfacecolor="none",
368 markeredgecolor=cog_color,
369 markersize=10,
370 markeredgewidth=2,
371 label="COG (fov_long, fov_lat)",
372 linestyle="None",
373 ),
374 Line2D(
375 [0],
376 [0],
377 color=ellipse_color,
378 linewidth=2,
379 linestyle="--",
380 label="Hillas Ellipse",
381 ),
382 Line2D([0], [0], color=length_color, linewidth=2, label="Length"),
383 Line2D([0], [0], color=width_color, linewidth=2, label="Width"),
384 Line2D([0], [0], color=angle_color, linewidth=2, label="Angle ψ"),
385 Line2D(
386 [0],
387 [0],
388 color=radial_color,
389 linewidth=2,
390 linestyle=":",
391 label="r, φ (radial)",
392 ),
393 ]
394 ax.legend(handles=legend_elements, loc="lower left", fontsize=10, framealpha=0.9)
395
396 plt.tight_layout()
397
398 return fig, hillas
Simulate an event#
405 subarray = SubarrayDescription.read("dataset://gamma_prod5.simtel.zst")
406 geom = subarray.tel[1].camera.geometry.transform_to(TelescopeFrame())
407
408 # Define Gaussian model parameters
409 x0 = -0.6 * u.deg
410 y0 = 1.2 * u.deg
411 sigma_length = 0.8 * u.deg
412 sigma_width = 0.2 * u.deg
413 psi = 65.0 * u.deg
414
415 model = Gaussian(x0, y0, sigma_length, sigma_width, psi)
416 image = model.generate_image(geom, intensity=10000)[0]
Hillas Parameters#
Hillas parameters and their meaning. Note that they are calculated after image cleaning.
Attribute |
Description |
|---|---|
intensity |
total intensity (size) |
skewness |
measure of the asymmetry |
kurtosis |
measure of the tailedness |
fov_lon |
longitude angle in a spherical system centered on the pointing position (deg) |
fov_lat |
latitude angle in a spherical system centered on the pointing position (deg) |
r |
radial coordinate of centroid (deg) |
phi |
polar coordinate of centroid (deg) |
length |
standard deviation along the major-axis (deg) |
length_uncertainty |
uncertainty of length (deg) |
width |
standard spread along the minor-axis (deg) |
width_uncertainty |
uncertainty of width (deg) |
psi |
rotation angle of ellipse (deg) |
psi_uncertainty |
uncertainty of psi (deg) |
Total running time of the script: (0 minutes 1.376 seconds)