Getting Started with ctapipe

This hands-on was presented at the Paris CTA Consoritum meeting (K. Kosack)

Part 1: load and loop over data

[1]:
from ctapipe.io import EventSource
from ctapipe import utils
from matplotlib import pyplot as plt
import numpy as np
%matplotlib inline
[2]:
path = utils.get_dataset_path("gamma_prod5.simtel.zst")
[3]:
source = EventSource(path, max_events=5)

for event in source:
    print(event.count, event.index.event_id, event.simulation.shower.energy)
0 4009 0.07287972420454025 TeV
1 5101 1.7454713582992554 TeV
2 5103 1.7454713582992554 TeV
3 5104 1.7454713582992554 TeV
4 5105 1.7454713582992554 TeV
[4]:
event
[4]:
ctapipe.containers.ArrayEventContainer:
                       index.*: event indexing information with default None
                          r0.*: Raw Data with default None
                          r1.*: R1 Calibrated Data with default None
                         dl0.*: DL0 Data Volume Reduced Data with default None
                         dl1.*: DL1 Calibrated image with default None
                         dl2.*: DL2 reconstruction info with default None
                  simulation.*: Simulated Event Information with default None
                                with type <class
                                'ctapipe.containers.SimulatedEventContainer'>
                     trigger.*: central trigger information with default None
                         count: number of events processed with default 0
                    pointing.*: Array and telescope pointing positions with
                                default None
                 calibration.*: Container for calibration coefficients for the
                                current event with default None
                         mon.*: container for event-wise monitoring data (MON)
                                with default None
[5]:
event.r1
[5]:
ctapipe.containers.R1Container:
                        tel[*]: map of tel_id to R1CameraContainer with default
                                None
[6]:
for event in EventSource(path, max_events=5):
    print(event.count, event.r1.tel.keys())
0 dict_keys([9, 14, 104])
1 dict_keys([26, 27, 29, 47, 53, 59, 69, 73, 121, 122, 124, 148, 162, 166])
2 dict_keys([82, 92, 175, 179])
3 dict_keys([17, 26, 27, 29, 43, 47, 53, 57, 69, 112, 121, 122, 124, 125, 128, 140, 144, 148, 162, 166])
4 dict_keys([3, 4, 8, 9, 14, 15, 23, 24, 25, 29, 32, 36, 42, 46, 52, 56, 103, 104, 109, 110, 118, 119, 120, 126, 127, 129, 130, 137, 139, 143, 161])
[7]:
event.r0.tel[3]
[7]:
ctapipe.containers.R0CameraContainer:
                      waveform: numpy array containing ADC samples(n_channels,
                                n_pixels, n_samples) with default None
[8]:
r0tel = event.r0.tel[3]
[9]:
r0tel.waveform
[9]:
array([[[432, 413, 395, ..., 382, 387, 382],
        [419, 420, 428, ..., 392, 396, 405],
        [381, 392, 413, ..., 375, 387, 387],
        ...,
        [391, 384, 385, ..., 379, 395, 404],
        [370, 376, 388, ..., 391, 383, 383],
        [393, 411, 413, ..., 406, 392, 396]],

       [[413, 396, 398, ..., 397, 401, 395],
        [392, 407, 403, ..., 403, 388, 399],
        [405, 403, 393, ..., 393, 400, 409],
        ...,
        [399, 406, 407, ..., 391, 405, 397],
        [401, 399, 396, ..., 395, 400, 388],
        [399, 393, 403, ..., 395, 407, 401]]], dtype=uint16)
[10]:
r0tel.waveform.shape
[10]:
(2, 1855, 40)

note that this is (\(N_{channels}\), \(N_{pixels}\), \(N_{samples}\))

[11]:
plt.pcolormesh(r0tel.waveform[0])
[11]:
<matplotlib.collections.QuadMesh at 0x7f9dfa67d610>
../_images/tutorials_ctapipe_handson_13_1.png
[12]:
brightest_pixel = np.argmax(r0tel.waveform[0].sum(axis=1))
print(f"pixel {brightest_pixel} has sum {r0tel.waveform[0,1535].sum()}")
pixel 1831 has sum 15824
[13]:
plt.plot(r0tel.waveform[0,brightest_pixel], label="channel 0 (high-gain)")
plt.plot(r0tel.waveform[1,brightest_pixel], label="channel 1 (low-gain)")
plt.legend()
[13]:
<matplotlib.legend.Legend at 0x7f9df8572df0>
../_images/tutorials_ctapipe_handson_15_1.png
[14]:
from ipywidgets import interact

@interact
def view_waveform(chan=0, pix_id=brightest_pixel):
    plt.plot(r0tel.waveform[chan, pix_id])

try making this compare 2 waveforms

Part 2: Explore the instrument description

This is all well and good, but we don’t really know what camera or telescope this is… how do we get instrumental description info?

Currently this is returned inside the event (it will soon change to be separate in next version or so)

[15]:
subarray = source.subarray
[16]:
subarray
[16]:
SubarrayDescription(name='MonteCarloArray', n_tels=180)
[17]:
subarray.peek()
../_images/tutorials_ctapipe_handson_21_0.png
[18]:
subarray.to_table()
[18]:
Table length=180
tel_idnametypepos_xpos_ypos_zcamera_nameoptics_namecamera_indexoptics_indextel_description
mmm
int16str5str3float32float32float32str9str5int64int64str17
1LSTLST-20.643-64.81734.3LSTCamLST21LST_LST_LSTCam
2LSTLST79.993996-0.76829.4LSTCamLST21LST_LST_LSTCam
3LSTLST-19.39665.231.0LSTCamLST21LST_LST_LSTCam
4LSTLST-120.0331.15133.1LSTCamLST21LST_LST_LSTCam
5MSTMST-0.017-0.00124.35FlashCamMST00MST_MST_FlashCam
6MSTMST-1.468-151.22131.0FlashCamMST00MST_MST_FlashCam
7MSTMST-3.1379998-325.24539.0FlashCamMST00MST_MST_FlashCam
8MSTMST1.4339999151.2225.0FlashCamMST00MST_MST_FlashCam
9MSTMST3.1039999325.24323.5FlashCamMST00MST_MST_FlashCam
.................................
171ASTRISST260.0920.045.0CHECASTRI12SST_ASTRI_CHEC
172ASTRISST260.0-920.065.0CHECASTRI12SST_ASTRI_CHEC
173ASTRISST-500.0815.015.0CHECASTRI12SST_ASTRI_CHEC
174ASTRISST-500.0-815.075.0CHECASTRI12SST_ASTRI_CHEC
175ASTRISST500.0815.045.0CHECASTRI12SST_ASTRI_CHEC
176ASTRISST500.0-815.053.0CHECASTRI12SST_ASTRI_CHEC
177ASTRISST-810.0655.012.0CHECASTRI12SST_ASTRI_CHEC
178ASTRISST-810.0-655.068.0CHECASTRI12SST_ASTRI_CHEC
179ASTRISST810.0655.020.0CHECASTRI12SST_ASTRI_CHEC
180ASTRISST810.0-655.041.0CHECASTRI12SST_ASTRI_CHEC
[19]:
subarray.tel[2]
[19]:
TelescopeDescription(type='LST', optics_name='LST', camera_name='LSTCam')
[20]:
subarray.tel[2].camera
[20]:
CameraDescription(name=LSTCam, geometry=LSTCam, readout=LSTCam)
[21]:
subarray.tel[2].optics
[21]:
OpticsDescription(name=LST, size_type=LST, reflector_shape=PARABOLIC, equivalent_focal_length=28.00 m, effective_focal_length=29.31 m, n_mirrors=1, mirror_area=386.73 m2)
[22]:
tel = subarray.tel[2]
[23]:
tel.camera
[23]:
CameraDescription(name=LSTCam, geometry=LSTCam, readout=LSTCam)
[24]:
tel.optics
[24]:
OpticsDescription(name=LST, size_type=LST, reflector_shape=PARABOLIC, equivalent_focal_length=28.00 m, effective_focal_length=29.31 m, n_mirrors=1, mirror_area=386.73 m2)
[25]:
tel.camera.geometry.pix_x
[25]:
$[0,~-0.037796701,~-0.047245475,~\dots,~0.67088033,~-0.45356484,~-0.50081024] \; \mathrm{m}$
[26]:
tel.camera.geometry.to_table()
[26]:
Table length=1855
pix_idpix_xpix_ypix_area
mmm2
int64float64float64float64
00.00.00.002079326892271638
1-0.03779670120038493-0.032732429036768330.002079326892271638
2-0.0472454749632934550.016366662083976240.002079326892271638
3-0.0094487737629085250.049099091120744570.002079326892271638
40.037796701200384930.032732429036768330.002079326892271638
50.047245474963293455-0.016366662083976240.002079326892271638
60.009448773762908525-0.049099091120744570.002079326892271638
7-0.07559330728842002-0.065464839769835820.002079326892271638
8-0.08504208105132854-0.0163657486490912420.002079326892271638
............
1845-0.70867688432060720.93288632847675070.002079326892271638
1846-0.67088033276595170.96561882545319940.002079326892271638
18470.453564836008831851.08017546211512070.002079326892271638
18480.50081024449666991.06380875119878530.002079326892271638
18491.16224173440920890.147289060475024020.002079326892271638
18501.17169050746812880.098189973012446750.002079326892271638
18510.7086768843206072-0.93288632847675070.002079326892271638
18520.6708803327659517-0.96561882545319940.002079326892271638
1853-0.45356483600883185-1.08017546211512070.002079326892271638
1854-0.5008102444966699-1.06380875119878530.002079326892271638
[27]:
tel.optics.mirror_area
[27]:
$386.73325 \; \mathrm{m^{2}}$
[28]:
from ctapipe.visualization import CameraDisplay
[29]:
disp = CameraDisplay(tel.camera.geometry)
../_images/tutorials_ctapipe_handson_33_0.png
[30]:
disp = CameraDisplay(tel.camera.geometry)
disp.image = r0tel.waveform[0,:,10]  # display channel 0, sample 0 (try others like 10)
../_images/tutorials_ctapipe_handson_34_0.png

** aside: ** show demo using a CameraDisplay in interactive mode in ipython rather than notebook

Part 3: Apply some calibration and trace integration

[31]:
from ctapipe.calib import CameraCalibrator
[32]:
calib = CameraCalibrator(subarray=subarray)
[33]:
for event in EventSource(path, max_events=5):
    calib(event) # fills in r1, dl0, and dl1
    print(event.dl1.tel.keys())
dict_keys([9, 14, 104])
dict_keys([26, 27, 29, 47, 53, 59, 69, 73, 121, 122, 124, 148, 162, 166])
dict_keys([82, 92, 175, 179])
dict_keys([17, 26, 27, 29, 43, 47, 53, 57, 69, 112, 121, 122, 124, 125, 128, 140, 144, 148, 162, 166])
dict_keys([3, 4, 8, 9, 14, 15, 23, 24, 25, 29, 32, 36, 42, 46, 52, 56, 103, 104, 109, 110, 118, 119, 120, 126, 127, 129, 130, 137, 139, 143, 161])
[34]:
event.dl1.tel[3]
[34]:
ctapipe.containers.DL1CameraContainer:
                         image: Numpy array of camera image, after waveform
                                extraction.Shape: (n_pixel) with default None as
                                a 1-D array with dtype float32
                     peak_time: Numpy array containing position of the peak of
                                the pulse as determined by the extractor. Shape:
                                (n_pixel, ) with default None as a 1-D array
                                with dtype float32
                    image_mask: Boolean numpy array where True means the pixel
                                has passed cleaning. Shape: (n_pixel, ) with
                                default None as a 1-D array with dtype bool
                      is_valid: True if image extraction succeeded, False if
                                failed or in the case of TwoPass methods, that
                                the first pass only was returned. with default
                                False
                    parameters: Image parameters with default None with type
                                <class
                                'ctapipe.containers.ImageParametersContainer'>
[35]:
dl1tel = event.dl1.tel[3]
[36]:
dl1tel.image.shape # note this will be gain-selected in next version, so will be just 1D array of 1855
[36]:
(1855,)
[37]:
dl1tel.peak_time
[37]:
array([27.381046, 17.746258, 18.80205 , ..., 11.64868 , 31.961287,
       29.717525], dtype=float32)
[38]:
CameraDisplay(tel.camera.geometry, image=dl1tel.image)
[38]:
<ctapipe.visualization.mpl_camera.CameraDisplay at 0x7f9dfc968550>
../_images/tutorials_ctapipe_handson_44_1.png
[39]:
CameraDisplay(tel.camera.geometry, image=dl1tel.peak_time)
[39]:
<ctapipe.visualization.mpl_camera.CameraDisplay at 0x7f9df83ec580>
../_images/tutorials_ctapipe_handson_45_1.png

Now for Hillas Parameters

[40]:
from ctapipe.image import hillas_parameters, tailcuts_clean
[41]:
image = dl1tel.image
mask = tailcuts_clean(tel.camera.geometry, image, picture_thresh=10, boundary_thresh=5)
mask
[41]:
array([False, False, False, ..., False, False, False])
[42]:
CameraDisplay(tel.camera.geometry, image=mask)
[42]:
<ctapipe.visualization.mpl_camera.CameraDisplay at 0x7f9df8539490>
../_images/tutorials_ctapipe_handson_49_1.png
[43]:
cleaned = image.copy()
cleaned[~mask] = 0
[44]:
disp = CameraDisplay(tel.camera.geometry, image=cleaned)
disp.cmap = plt.cm.coolwarm
disp.add_colorbar()
plt.xlim(0.5, 1.0)
plt.ylim(-1.0, 0.0)
[44]:
(-1.0, 0.0)
../_images/tutorials_ctapipe_handson_51_1.png
[45]:
params = hillas_parameters(tel.camera.geometry, cleaned)
print(params)
{'intensity': 32.78351974487305,
 'kurtosis': 1.2333006676347844,
 'length': <Quantity 0.02219065 m>,
 'length_uncertainty': <Quantity 0.00093599 m>,
 'phi': <Angle -0.88747705 rad>,
 'psi': <Angle -1.56935711 rad>,
 'r': <Quantity 1.11337613 m>,
 'skewness': -0.2716153815520628,
 'width': <Quantity 0.01614423 m>,
 'width_uncertainty': <Quantity 0.00226475 m>,
 'x': <Quantity 0.70295289 m>,
 'y': <Quantity -0.86340237 m>}
[46]:
disp = CameraDisplay(tel.camera.geometry, image=cleaned)
disp.cmap = plt.cm.coolwarm
disp.add_colorbar()
plt.xlim(0.5, 1.0)
plt.ylim(-1.0, 0.0)
disp.overlay_moments(params, color='white', lw=2)
../_images/tutorials_ctapipe_handson_53_0.png

Part 4: Let’s put it all together:

  • loop over events, selecting only telescopes of the same type (e.g. LST:LSTCam)

  • for each event, apply calibration/trace integration

  • calculate Hillas parameters

  • write out all hillas paremeters to a file that can be loaded with Pandas

first let’s select only those telescopes with LST:LSTCam

[47]:
subarray.telescope_types
[47]:
(TelescopeDescription(type='LST', optics_name='LST', camera_name='LSTCam'),
 TelescopeDescription(type='MST', optics_name='MST', camera_name='NectarCam'),
 TelescopeDescription(type='MST', optics_name='MST', camera_name='FlashCam'),
 TelescopeDescription(type='SST', optics_name='ASTRI', camera_name='CHEC'))
[48]:
subarray.get_tel_ids_for_type("LST_LST_LSTCam")
[48]:
(1, 2, 3, 4)

Now let’s write out program

[49]:
data = utils.get_dataset_path("gamma_prod5.simtel.zst")
source = EventSource(data) # remove the max_events limit to get more stats
[50]:
for event in source:
    calib(event)

    for tel_id, tel_data in event.dl1.tel.items():
        tel = source.subarray.tel[tel_id]
        mask = tailcuts_clean(tel.camera.geometry, tel_data.image)
        if np.count_nonzero(mask) > 0:
            params = hillas_parameters(tel.camera.geometry[mask], tel_data.image[mask])
[51]:
from ctapipe.io import HDF5TableWriter

[52]:
with HDF5TableWriter(filename='hillas.h5', group_name='dl1', overwrite=True) as writer:

    source = EventSource(data, allowed_tels=[1,2,3,4],  max_events=10)
    for event in source:
        calib(event)

        for tel_id, tel_data in event.dl1.tel.items():
            tel = source.subarray.tel[tel_id]
            mask = tailcuts_clean(tel.camera.geometry, tel_data.image)
            params = hillas_parameters(tel.camera.geometry[mask], tel_data.image[mask])
            writer.write("hillas", params)

We can now load in the file we created and plot it

[53]:
!ls *.h5
hillas.h5
[54]:
import pandas as pd

hillas = pd.read_hdf("hillas.h5", key='/dl1/hillas')
hillas
[54]:
intensity skewness kurtosis x y r phi length length_uncertainty width width_uncertainty psi
0 47.098688 -0.530399 1.881643 0.695206 -0.847657 1.096282 -50.643096 0.032743 0.002240 0.026627 0.002372 -51.309393
1 126.314263 0.020218 1.742869 0.464446 -0.999899 1.102500 -65.085472 0.055256 0.002119 0.025090 0.001144 -67.386572
[55]:
_ = hillas.hist(figsize=(8,8))
../_images/tutorials_ctapipe_handson_66_0.png

If you do this yourself, chose a larger file to loop over more events to get better statistics