Explore Calibrated Data

import ctapipe
from ctapipe.utils.datasets import get_dataset_path
from ctapipe.io import EventSource, EventSeeker
from ctapipe.visualization import CameraDisplay
from ctapipe.instrument import CameraGeometry
from matplotlib import pyplot as plt
from astropy import units as u
import numpy as np
%matplotlib inline
plt.style.use("ggplot")

print(ctapipe.__version__)
print(ctapipe.__file__)

0.12.1.dev1+g71bebb8f
/home/docs/checkouts/readthedocs.org/user_builds/ctapipe/envs/v0.12.0-rtd/lib/python3.8/site-packages/ctapipe/__init__.py

Let’s first open a raw event file and get an event out of it:

filename = get_dataset_path("gamma_test.simtel.gz")
source = EventSource(filename, max_events=2)

for event in source:
    print(event.index.event_id)

/home/docs/checkouts/readthedocs.org/user_builds/ctapipe/envs/v0.12.0-rtd/lib/python3.8/site-packages/ctapipe/io/simteleventsource.py:76: UnknownPixelShapeWarning: Unkown pixel_shape -1.0 for camera_type LSTCam
  warnings.warn(
/home/docs/checkouts/readthedocs.org/user_builds/ctapipe/envs/v0.12.0-rtd/lib/python3.8/site-packages/ctapipe/io/simteleventsource.py:76: UnknownPixelShapeWarning: Unkown pixel_shape -1.0 for camera_type NectarCam
  warnings.warn(
/home/docs/checkouts/readthedocs.org/user_builds/ctapipe/envs/v0.12.0-rtd/lib/python3.8/site-packages/ctapipe/io/simteleventsource.py:76: UnknownPixelShapeWarning: Unkown pixel_shape -1.0 for camera_type CHEC
  warnings.warn(
/home/docs/checkouts/readthedocs.org/user_builds/ctapipe/envs/v0.12.0-rtd/lib/python3.8/site-packages/ctapipe/io/simteleventsource.py:76: UnknownPixelShapeWarning: Unkown pixel_shape -1.0 for camera_type SCTCam
  warnings.warn(

408
409

filename

PosixPath('/home/docs/.cache/ctapipe/cccta-dataserver.in2p3.fr/data/ctapipe-extra/v0.3.3/gamma_test.simtel.gz')

source

SimTelEventSource

Read events from a SimTelArray data file (in EventIO format).

allowed_tels	None	list of allowed tel_ids, others will be ignored. If None, all telescopes in the input stream will be included (default: None)
back_seekable	False	Require the event source to be backwards seekable. This will reduce in slower read speed for gzipped files and is not possible for zstd compressed files (default: False)
calib_scale	1.0	Factor to transform ADC counts into number of photoelectrons. Corrects the DC_to_PHE factor. (default: 1.0)
calib_shift	0.0	Factor to shift the R1 photoelectron samples. Can be used to simulate mis-calibration. (default: 0.0)
focal_length_choice	nominal	if both nominal and effective focal lengths are available in the SimTelArray file, which one to use when constructing the SubarrayDescription (which will be used in CameraFrame to TelescopeFrame coordinate transforms. The 'nominal' focal length is the one used during the simulation, the 'effective' focal length is computed using specialized ray-tracing from a point light source (default: nominal)
gain_selector_type	ThresholdGainSelector	GainSelector to use. (default: ThresholdGainSelector)
input_url	/home/docs/.cache/ctapipe/cccta-dataserver.in2p3.fr/data/ctapipe-extra/v0.3.3/gamma_test.simtel.gz	Path to the input file containing events. (default: traitlets.Undefined)
max_events	2	Maximum number of events that will be read from the file (default: None)
skip_calibration_events	True	Skip calibration events (default: True)

event

ctapipe.containers.ArrayEventContainer:
                       index.*: event indexing information
                          r0.*: Raw Data
                          r1.*: R1 Calibrated Data
                         dl0.*: DL0 Data Volume Reduced Data
                         dl1.*: DL1 Calibrated image
                         dl2.*: DL2 reconstruction info
                  simulation.*: Simulated Event Information
                     trigger.*: central trigger information
                         count: number of events processed
                    pointing.*: Array and telescope pointing positions
                 calibration.*: Container for calibration coefficients for the
                                current event
                         mon.*: container for event-wise monitoring data (MON)

print(event.r1)

{'tel': {11: {'selected_gain_channel': array([0, 0, 0, ..., 0, 0, 0], dtype=int8),
              'waveform': array([[-0.08697608,  0.18923593,  0.18923593, ...,  0.58820885,
         0.03578481, -0.17904675],
       [ 0.65603745,  0.65603745,  0.16448696, ..., -0.22875345,
        -0.5236837 , -0.13044335],
       [-0.31448284, -0.15485168, -0.12292545, ..., -0.21870415,
        -0.15485168, -0.28255662],
       ...,
       [-0.2288092 ,  0.1161234 , -0.04066414, ...,  0.02205088,
        -0.07202165,  0.02205088],
       [ 0.0596489 , -0.06216849, -0.09262283, ..., -0.18398587,
        -0.36671194, -0.12307718],
       [-0.38044545, -0.24650918,  0.05484743, ...,  0.52362436,
         0.55710846,  0.52362436]], dtype=float32)},
         21: {'selected_gain_channel': array([0, 0, 0, ..., 0, 0, 0], dtype=int8),
              'waveform': array([[-0.23289204, -0.2635209 , -0.07974766, ..., -0.04911879,
         0.01213896, -0.17163429],
       [ 1.1236377 ,  0.5958728 ,  0.25437793, ..., -0.02502699,
        -0.05607199, -0.18025196],
       [-0.17667109, -0.04912165,  0.3973014 , ...,  0.58862555,
         0.6205129 ,  0.716175  ],
       ...,
       [ 0.40350705,  0.8413429 ,  0.77879494, ..., -0.31579474,
        -0.00305482,  0.528603  ],
       [-0.09964581, -0.0381696 ,  0.08478284, ...,  0.3306877 ,
         0.20773527,  0.23847337],
       [ 0.42958432,  0.15098573,  0.27480733, ..., -0.03474665,
        -0.00379125, -0.06570205]], dtype=float32)},
         24: {'selected_gain_channel': array([0, 0, 0, ..., 0, 0, 0], dtype=int8),
              'waveform': array([[ 0.20489624,  0.33959898, -0.13186066, ..., -0.40126616,
        -0.16553634, -0.5359689 ],
       [-0.22010921, -0.12690239,  0.18378705, ..., -0.344385  ,
        -0.2822471 ,  0.02844233],
       [-0.09333542, -0.31169948, -0.18692002, ..., -0.31169948,
        -0.4052841 , -0.46767384],
       ...,
       [-0.21383263, -0.43987572, -0.21383263, ..., -0.31070825,
        -0.11695702, -0.3430001 ],
       [ 0.16586055,  0.07120538,  0.13430882, ...,  0.57603294,
         0.85999846,  0.8915502 ],
       [-0.05169091, -0.11653713,  0.04557843, ..., -0.18138336,
        -0.18138336, -0.37592202]], dtype=float32)},
         26: {'selected_gain_channel': array([0, 0, 0, ..., 0, 0, 0], dtype=int8),
              'waveform': array([[ 0.10818457, -0.17343406, -0.07956119, ..., -0.17343406,
        -0.2985979 , -0.26730695],
       [ 0.28453484,  0.03107848, -0.41247013, ...,  0.12612462,
         0.09444257,  0.60135525],
       [-0.01938689, -0.05021886,  0.07310902, ...,  0.01144508,
         0.35059676,  0.13477296],
       ...,
       [-0.48684907, -0.29694626, -0.01209205, ...,  0.2411117 ,
         0.01955842, -0.1703444 ],
       [-0.15805551, -0.15805551,  0.2170638 , ..., -0.31435522,
        -0.47065493, -0.12679556],
       [-0.17201822, -0.23635136, -0.20418479, ...,  0.21398063,
         0.27831376,  0.37481347]], dtype=float32)},
         61: {'selected_gain_channel': array([0, 0, 0, ..., 0, 0, 0], dtype=int8),
              'waveform': array([[-0.09033689,  0.14919089, -0.21010077, ...,  0.029427  ,
        -0.21010077, -0.09033689],
       [-0.20168626,  0.16380924,  0.07243536, ..., -0.11031239,
         0.16380924,  0.25518313],
       [-0.02433773, -0.02433773, -0.14491166, ...,  0.33738407,
         0.09623621,  0.09623621],
       ...,
       [ 0.07372501, -0.3746946 , -0.10564283, ..., -0.0159589 ,
        -0.0159589 ,  0.16340894],
       [-0.09563459, -0.00609556, -0.00609556, ..., -0.09563459,
        -0.09563459, -0.18517363],
       [ 0.0242704 ,  0.12670045, -0.07815965, ...,  0.53642064,
         0.33156055, -0.28301972]], dtype=float32)},
         63: {'selected_gain_channel': array([0, 0, 0, ..., 0, 0, 0], dtype=int8),
              'waveform': array([[ 8.3050177e-02, -8.0907416e-01, -2.1432461e-01, ...,
         8.2648712e-01,  2.3173757e-01,  8.3050177e-02],
       [ 8.3567634e-02, -7.2033242e-03, -3.7028715e-01, ...,
         2.6510957e-01, -2.7951619e-01, -4.6105811e-01],
       [-9.6725412e-02, -6.9741560e-03, -9.6725412e-02, ...,
        -4.5573044e-01, -9.6725412e-02,  1.7252836e-01],
       ...,
       [ 1.3219908e-01, -6.0026228e-02,  3.2442439e-01, ...,
         3.6086425e-02, -1.5613888e-01,  3.6086425e-02],
       [-8.2004257e-02, -4.6903652e-04, -4.6903652e-04, ...,
         1.6260140e-01, -4.6903652e-04,  8.1066184e-02],
       [-1.8674533e-01,  1.0895659e-01, -1.8674533e-01, ...,
         1.0895659e-01, -2.8531262e-01, -8.8178024e-02]], dtype=float32)},
         118: {'selected_gain_channel': array([0, 0, 0, ..., 0, 0, 0], dtype=int8),
               'waveform': array([[-0.00312013, -0.00312013, -0.01841514, ..., -0.01841514,
        -0.00312013, -0.00312013],
       [-0.01518433, -0.01518433, -0.01518433, ..., -0.01518433,
        -0.01518433, -0.01518433],
       [-0.00777916, -0.00777916, -0.00777916, ..., -0.00777916,
        -0.00777916, -0.00777916],
       ...,
       [-0.01964806, -0.01964806, -0.01964806, ..., -0.01964806,
        -0.00275011, -0.01964806],
       [-0.01970303, -0.01970303, -0.01970303, ..., -0.01970303,
        -0.01970303, -0.01970303],
       [-0.01387683,  0.00185392, -0.01387683, ...,  0.12769988,
         0.08050765,  0.04904616]], dtype=float32)},
         119: {'selected_gain_channel': array([0, 0, 0, ..., 0, 0, 0], dtype=int8),
               'waveform': array([[-0.005349  , -0.005349  , -0.005349  , ..., -0.005349  ,
        -0.005349  , -0.005349  ],
       [-0.0152637 , -0.0152637 , -0.0152637 , ...,  0.17803247,
         0.17803247,  0.12970842],
       [-0.00617304, -0.00617304, -0.00617304, ..., -0.00617304,
        -0.00617304, -0.00617304],
       ...,
       [-0.0155689 , -0.0155689 , -0.0155689 , ..., -0.0155689 ,
        -0.0155689 , -0.0155689 ],
       [-0.01582186, -0.01582186, -0.01582186, ..., -0.01582186,
        -0.01582186, -0.01582186],
       [-0.01675882, -0.01675882, -0.01675882, ..., -0.01675882,
        -0.01675882, -0.01675882]], dtype=float32)}}}

Perform basic calibration:

Here we will use a CameraCalibrator which is just a simple wrapper that runs the three calibraraton and trace-integration phases of the pipeline, taking the data from levels:

R0 → R1 → DL0 → DL1

You could of course do these each separately, by using the classes R1Calibrator, DL0Reducer, and DL1Calibrator. Note that we have not specified any configuration to the CameraCalibrator, so it will be using the default algorithms and thresholds, other than specifying that the product is a “HESSIOR1Calibrator” (hopefully in the near future that will be automatic).

from ctapipe.calib import CameraCalibrator

calib = CameraCalibrator(subarray=source.subarray)
calib(event)

Now the r1, dl0 and dl1 containers are filled in the event

• r1.tel[x]: contains the “r1-calibrated” waveforms, after gain-selection, pedestal subtraciton, and gain-correction

• dl0.tel[x]: is the same but with optional data volume reduction (some pixels not filled), in this case this is not performed by default, so it is the same as r1

• dl1.tel[x]: contains the (possibly re-calibrated) waveforms as dl0, but also the time-integrated image that has been calculated using a ImageExtractor (a NeighborPeakWindowSum by default)

for tel_id in event.dl1.tel:
    print("TEL{:03}: {}".format(tel_id, source.subarray.tel[tel_id]))
    print("  - r0  wave shape  : {}".format(event.r0.tel[tel_id].waveform.shape))
    print("  - r1  wave shape  : {}".format(event.r1.tel[tel_id].waveform.shape))
    print("  - dl1 image shape : {}".format(event.dl1.tel[tel_id].image.shape))

TEL011: MST_MST_NectarCam
  - r0  wave shape  : (1, 1855, 30)
  - r1  wave shape  : (1855, 30)
  - dl1 image shape : (1855,)
TEL021: MST_MST_NectarCam
  - r0  wave shape  : (1, 1855, 30)
  - r1  wave shape  : (1855, 30)
  - dl1 image shape : (1855,)
TEL024: MST_MST_NectarCam
  - r0  wave shape  : (1, 1855, 30)
  - r1  wave shape  : (1855, 30)
  - dl1 image shape : (1855,)
TEL026: MST_MST_NectarCam
  - r0  wave shape  : (1, 1855, 30)
  - r1  wave shape  : (1855, 30)
  - dl1 image shape : (1855,)
TEL061: SST_ASTRI_CHEC
  - r0  wave shape  : (1, 2048, 25)
  - r1  wave shape  : (2048, 25)
  - dl1 image shape : (2048,)
TEL063: SST_ASTRI_CHEC
  - r0  wave shape  : (1, 2048, 25)
  - r1  wave shape  : (2048, 25)
  - dl1 image shape : (2048,)
TEL118: MST_SCT_SCTCam
  - r0  wave shape  : (1, 11328, 32)
  - r1  wave shape  : (11328, 32)
  - dl1 image shape : (11328,)
TEL119: MST_SCT_SCTCam
  - r0  wave shape  : (1, 11328, 32)
  - r1  wave shape  : (11328, 32)
  - dl1 image shape : (11328,)

Some image processing:

Let’s look at the image

from ctapipe.visualization import CameraDisplay
tel_id = sorted(event.r1.tel.keys())[1]
sub = source.subarray
geometry = sub.tel[tel_id].camera.geometry
image = event.dl1.tel[tel_id].image

disp = CameraDisplay(geometry, image=image)

from ctapipe.image import tailcuts_clean, hillas_parameters

mask = tailcuts_clean(geometry, image, picture_thresh=10, boundary_thresh=5, min_number_picture_neighbors=2)
cleaned = image.copy()
cleaned[~mask] = 0
disp = CameraDisplay(geometry, image=cleaned)

params = hillas_parameters(geometry, cleaned)
print(params)
params

{'intensity': 40.55668926239014,
 'kurtosis': 1.7355258945159615,
 'length': <Quantity 0.0329767 m>,
 'length_uncertainty': <Quantity 0.00222047 m>,
 'phi': <Angle 2.85056754 rad>,
 'psi': <Angle -0.32029021 rad>,
 'r': <Quantity 0.52409634 m>,
 'skewness': 0.15350147214718432,
 'width': <Quantity 0.01579441 m>,
 'width_uncertainty': <Quantity 0.00150044 m>,
 'x': <Quantity -0.50205821 m>,
 'y': <Quantity 0.15038126 m>}

ctapipe.containers.CameraHillasParametersContainer:
                     intensity: total intensity (size)
                      skewness: measure of the asymmetry
                      kurtosis: measure of the tailedness
                             x: centroid x coordinate [m]
                             y: centroid x coordinate [m]
                             r: radial coordinate of centroid [m]
                           phi: polar coordinate of centroid [deg]
                        length: standard deviation along the major-axis [m]
            length_uncertainty: uncertainty of length [m]
                         width: standard spread along the minor-axis [m]
             width_uncertainty: uncertainty of width [m]
                           psi: rotation angle of ellipse [deg]

params = hillas_parameters(geometry, cleaned)

plt.figure(figsize=(10,10))
disp = CameraDisplay(geometry, image=image)
disp.add_colorbar()
disp.overlay_moments(params, color='red', lw=3)
disp.highlight_pixels(mask, color='white', alpha=0.3, linewidth=2)

plt.xlim(params.x.to_value(u.m) - 0.5, params.x.to_value(u.m) + 0.5)
plt.ylim(params.y.to_value(u.m) - 0.5, params.y.to_value(u.m) + 0.5)

(-0.34961873844871744, 0.6503812615512825)

source.metadata

{'is_simulation': False}

More complex image processing:

Let’s now explore how stereo reconstruction works.

first, look at a summed image from multiple telescopes

For this, we want to use a CameraDisplay again, but since we can’t sum and display images with different cameras, we’ll just sub-select images from a particular camera type

These are the telescopes that are in this event:

tels_in_event = set(event.dl1.tel.keys())  # use a set here, so we can intersect it later
tels_in_event

{11, 21, 24, 26, 61, 63, 118, 119}

cam_ids = set(sub.get_tel_ids_for_type("MST_MST_NectarCam"))
cam_ids

{5,
 6,
 7,
 8,
 9,
 10,
 11,
 12,
 13,
 14,
 15,
 16,
 17,
 18,
 19,
 20,
 21,
 22,
 23,
 24,
 25,
 26,
 27,
 28}

cams_in_event = tels_in_event.intersection(cam_ids)
first_tel_id = list(cams_in_event)[0]
tel = sub.tel[first_tel_id]
print("{}s in event: {}".format(tel, cams_in_event))

MST_MST_NectarCams in event: {24, 26, 11, 21}

Now, let’s sum those images:

image_sum = np.zeros_like(tel.camera.geometry.pix_x.value)  # just make an array of 0's in the same shape as the camera

for tel_id in cams_in_event:
    image_sum += event.dl1.tel[tel_id].image

And finally display the sum of those images

plt.figure(figsize=(8,8))

disp = CameraDisplay(tel.camera.geometry, image=image_sum)
disp.overlay_moments(params, with_label=False)
plt.title("Sum of {}x {}".format(len(cams_in_event), tel))

Text(0.5, 1.0, 'Sum of 4x MST_MST_NectarCam')

let’s also show which telescopes those were. Note that currently ArrayDisplay’s value field is a vector by tel_index, not tel_id, so we have to convert to a tel_index. (this may change in a future version to be more user-friendly)

from ctapipe.visualization import ArrayDisplay

nectarcam_subarray = sub.select_subarray(cam_ids, name="NectarCam")

hit_pattern = np.zeros(shape=nectarcam_subarray.num_tels)
hit_pattern[[nectarcam_subarray.tel_indices[x] for x in cams_in_event ]] = 100

plt.set_cmap(plt.cm.Accent)
plt.figure(figsize=(8,8))

ad = ArrayDisplay(nectarcam_subarray)
ad.values = hit_pattern
ad.add_labels()

<Figure size 432x288 with 0 Axes>