Make a theta-square plot

This is a basic example to analyze some events and make a \(\Theta^2\) plot

%matplotlib inline
from astropy import units as u
from astropy.coordinates.angle_utilities import angular_separation
from astropy.coordinates import SkyCoord, AltAz

import matplotlib.pyplot as plt
import numpy as np

from import EventSource
from ctapipe.visualization import CameraDisplay
from ctapipe.instrument import CameraGeometry
from ctapipe.calib import CameraCalibrator
from ctapipe.reco import HillasReconstructor
from ctapipe.image import hillas_parameters, tailcuts_clean
from ctapipe.utils import datasets

Get source events in MC dataset. Here we stop at 10 events, just to make this example run fast, but for real use, one would need more statistics.

filename = datasets.get_dataset_path("gamma_test_large.simtel.gz")
source = EventSource(filename, allowed_tels={1, 2, 3, 4}, max_events=10)
reco = HillasReconstructor()
calib = CameraCalibrator(subarray=source.subarray)
horizon_frame = AltAz()
off_angles = []
for event in source:

    # calibrating the event

    hillas_params = {}

    # pointing direction of the telescopes
    telescope_pointings = {}

    subarray = source.subarray

    # get hillas params for each event in different telescopes
    for tel_id in

        # telescope pointing direction
        telescope_pointings[tel_id] = SkyCoord(

        # Camera Geometry required for hillas parametrization
        camgeom =[tel_id].camera.geometry

        # note the [0] is for channel 0 which is high-gain channel
        image =[tel_id].image

        # Cleaning  of the image
        cleaned_image = image
        # create a clean mask of pixels above the threshold
        cleanmask = tailcuts_clean(
            camgeom, image, picture_thresh=10, boundary_thresh=5
        # set all rejected pixels to zero
        cleaned_image[~cleanmask] = 0

        # Calulate hillas parameters
        # It fails for empty pixels
            hillas = hillas_parameters(camgeom, cleaned_image)
        except Exception as e:

        if hillas.width.value > 0:
            hillas_params[tel_id] = hillas

    if len(hillas_params) < 2:

    array_pointing = SkyCoord(

    reco_result = reco.predict(
        hillas_params, source.subarray,
        array_pointing, telescope_pointings,

    # get angular offset between reconstructed shower direction and MC
    # generated shower direction
    mcshower = event.simulation.shower
    off_angle = angular_separation(, mcshower.alt,, reco_result.alt)

    # Appending all estimated off angles

calculate theta square for angles which are not nan

off_angles = np.array(off_angles)
thetasquare = off_angles[np.isfinite(off_angles)]**2

Plot the results

plt.hist(thetasquare, bins=10, range=[0,0.4])
plt.xlabel(r'$\theta^2$ (deg)')
plt.ylabel("# of events")

again, this plot is not beautiful since we have such low stats