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
/usr/local/lib/python3.8/site-packages/setuptools_scm/git.py:68: UserWarning: "/github/workspace" is shallow and may cause errors
  warnings.warn('"{}" is shallow and may cause errors'.format(wd.path))
[2]:
path = utils.get_dataset_path("gamma_test_large.simtel.gz")
[3]:
source = EventSource(path, max_events=4)

for event in source:
    print(event.count, event.index.event_id, event.simulation.shower.energy)
0 23703 0.5707105398178101 TeV
1 31007 1.8637498617172241 TeV
2 31010 1.8637498617172241 TeV
3 31012 1.8637498617172241 TeV
[4]:
event
[4]:
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.*: Reconstructed Shower Information
                  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)
[5]:
event.r0
[5]:
ctapipe.containers.R0Container:
                        tel[*]: map of tel_id to R0CameraContainer
[6]:
for event in EventSource(path, max_events=4):
    print(event.count, event.r0.tel.keys())
0 dict_keys([13, 21, 25, 34])
1 dict_keys([7, 13, 16, 19, 25, 28, 34, 36, 42])
2 dict_keys([28, 46, 58, 68])
3 dict_keys([2, 3, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 16, 18, 19, 20, 21])
[7]:
event.r0.tel[2]
[7]:
ctapipe.containers.R0CameraContainer:
                      waveform: numpy array containing ADC samples(n_channels,
                                n_pixels, n_samples)
[8]:
r0tel = event.r0.tel[3]
[9]:
r0tel.waveform
[9]:
array([[[267, 288, 286, ..., 297, 298, 312],
        [293, 285, 258, ..., 344, 373, 365],
        [333, 319, 289, ..., 306, 272, 261],
        ...,
        [255, 261, 270, ..., 258, 269, 271],
        [277, 312, 277, ..., 294, 323, 376],
        [319, 300, 347, ..., 283, 302, 265]],

       [[295, 304, 299, ..., 299, 298, 299],
        [298, 297, 300, ..., 300, 307, 304],
        [301, 305, 300, ..., 303, 300, 297],
        ...,
        [297, 296, 301, ..., 295, 297, 298],
        [295, 301, 296, ..., 300, 305, 302],
        [298, 300, 302, ..., 302, 297, 301]]], dtype=uint16)
[10]:
r0tel.waveform.shape
[10]:
(2, 1855, 30)

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

[11]:
plt.pcolormesh(r0tel.waveform[0])
[11]:
<matplotlib.collections.QuadMesh at 0x7fe49aa27c70>
../_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 344 has sum 8642
[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 0x7fe49a98f790>
../_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', num_tels=98)
[17]:
subarray.peek()
../_images/tutorials_ctapipe_handson_21_0.png
[18]:
subarray.to_table()
[18]:
Table length=98
tel_idpos_xpos_ypos_znametypenum_mirrorscamera_typetel_description
mmm
int16float64float64float64str5str3int64str8str18
1-20.065.016.0LSTLST1LSTCamLST_LST_LSTCam
2-20.0-65.016.0LSTLST1LSTCamLST_LST_LSTCam
380.00.016.0LSTLST1LSTCamLST_LST_LSTCam
4-120.00.016.0LSTLST1LSTCamLST_LST_LSTCam
50.00.010.0MSTMST1FlashCamMST_MST_FlashCam
60.0151.199996948242210.0MSTMST1FlashCamMST_MST_FlashCam
70.0-151.199996948242210.0MSTMST1FlashCamMST_MST_FlashCam
8146.6559906005859475.599998474121110.0MSTMST1FlashCamMST_MST_FlashCam
9146.65599060058594-75.599998474121110.0MSTMST1FlashCamMST_MST_FlashCam
...........................
89956.7870483398438739.8229980468755.0ASTRISST2ASTRICamSST_ASTRI_ASTRICam
90956.7870483398438-739.8229980468755.0ASTRISST2ASTRICamSST_ASTRI_ASTRICam
91-239.196990966796881109.73498535156255.0ASTRISST2ASTRICamSST_ASTRI_ASTRICam
92-239.19699096679688-1109.73498535156255.0ASTRISST2ASTRICamSST_ASTRI_ASTRICam
93-956.7870483398438739.8229980468755.0ASTRISST2ASTRICamSST_ASTRI_ASTRICam
94-956.7870483398438-739.8229980468755.0ASTRISST2ASTRICamSST_ASTRI_ASTRICam
951195.9840087890625369.91198730468755.0ASTRISST2ASTRICamSST_ASTRI_ASTRICam
961195.9840087890625-369.91198730468755.0ASTRISST2ASTRICamSST_ASTRI_ASTRICam
97-1195.9840087890625369.91198730468755.0ASTRISST2ASTRICamSST_ASTRI_ASTRICam
98-1195.9840087890625-369.91198730468755.0ASTRISST2ASTRICamSST_ASTRI_ASTRICam
[19]:
subarray.tel[2]
[19]:
TelescopeDescription(type=LST, name=LST, optics=LST, camera=LSTCam)
[20]:
subarray.tel[2].camera
[20]:
CameraDescription(camera_name=LSTCam, geometry=LSTCam, readout=LSTCam)
[21]:
subarray.tel[2].optics
[21]:
OpticsDescription(name=LST, equivalent_focal_length=28.00 m, num_mirrors=1, mirror_area=386.73 m2)
[22]:
tel = subarray.tel[2]
[23]:
tel.camera
[23]:
CameraDescription(camera_name=LSTCam, geometry=LSTCam, readout=LSTCam)
[24]:
tel.optics
[24]:
OpticsDescription(name=LST, equivalent_focal_length=28.00 m, num_mirrors=1, mirror_area=386.73 m2)
[25]:
tel.camera.geometry.pix_x
[25]:
$[0,~-0.0094487737,~-0.047244198,~\dots,~-0.6519913,~-0.6141959,~-0.62364468] \; \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.009448773714197630.0490990911301186360.002079326892271638
2-0.047244198246667620.016366907822831360.002079326892271638
3-0.037795424532469986-0.032732183307287280.002079326892271638
40.00944877371419763-0.0490990911301186360.002079326892271638
50.04724419824666762-0.016366907822831360.002079326892271638
60.0377954245324699860.032732183307287280.002079326892271638
70.06614174532306863-0.114565088253984980.002079326892271638
80.056692972312859496-0.065466000782033660.002079326892271638
............
1845-0.71814285557966-0.85104567723254190.002079326892271638
1846-0.6803376189247268-0.81831163031019890.002079326892271638
1847-0.6897863947508899-0.76921252820557840.002079326892271638
1848-0.6614400785362163-0.91650977598876320.002079326892271638
1849-0.6708888430985637-0.86741073241481950.002079326892271638
1850-0.7086940797534969-0.90014477933716240.002079326892271638
1851-0.6992453039273338-0.9492438814417830.002079326892271638
1852-0.6519913027100532-0.96560887809338370.002079326892271638
1853-0.6141958992088292-0.93287672349210340.002079326892271638
1854-0.6236446750349923-0.88377762138748280.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=4):
    calib(event) # fills in r1, dl0, and dl1
    print(event.dl1.tel.keys())
dict_keys([13, 21, 25, 34])
dict_keys([7, 13, 16, 19, 25, 28, 34, 36, 42])
dict_keys([28, 46, 58, 68])
dict_keys([2, 3, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 16, 18, 19, 20, 21])
[34]:
event.dl1.tel[3]
[34]:
ctapipe.containers.DL1CameraContainer:
                         image: Numpy array of camera image, after waveform
                                extraction.Shape: (n_pixel) as a 1-D array with
                                type float32
                     peak_time: Numpy array containing position of the peak of
                                the pulse as determined by the extractor. Shape:
                                (n_pixel) as a 1-D array with type float32
                    image_mask: Boolean numpy array where True means the pixel
                                has passed cleaning. Shape: (n_pixel) as a 1-D
                                array with type bool
                  parameters.*: Parameters derived from images
[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([34.393623, 49.123714, 37.32338 , ..., 50.50625 , 45.477562,
       42.10054 ], dtype=float32)
[38]:
CameraDisplay(tel.camera.geometry, image=dl1tel.image)
[38]:
<ctapipe.visualization.mpl_camera.CameraDisplay at 0x7fe49aa80a60>
../_images/tutorials_ctapipe_handson_44_1.png
[39]:
CameraDisplay(tel.camera.geometry, image=dl1tel.peak_time)
[39]:
<ctapipe.visualization.mpl_camera.CameraDisplay at 0x7fe492fd3fa0>
../_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 0x7fe492eb5850>
../_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(-1.0,0)
plt.ylim(0,1.0)
[44]:
(0.0, 1.0)
../_images/tutorials_ctapipe_handson_51_1.png
[45]:
params = hillas_parameters(tel.camera.geometry, cleaned)
print(params)
{'intensity': 24052.02079820633,
 'kurtosis': 5.966535690818616,
 'length': <Quantity 0.20528489 m>,
 'length_uncertainty': <Quantity 0.00147495 m>,
 'phi': <Angle 2.23935873 rad>,
 'psi': <Angle -1.32078288 rad>,
 'r': <Quantity 0.9278434 m>,
 'skewness': 1.5762077646356902,
 'width': <Quantity 0.07930658 m>,
 'width_uncertainty': <Quantity 0.00058711 m>,
 'x': <Quantity -0.57513164 m>,
 'y': <Quantity 0.72809132 m>}
[46]:
disp = CameraDisplay(tel.camera.geometry, image=cleaned)
disp.cmap = plt.cm.coolwarm
disp.add_colorbar()
plt.xlim(-1.0,0)
plt.ylim(0,1.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, name=LST, optics=LST, camera=LSTCam),
 TelescopeDescription(type=SST, name=ASTRI, optics=ASTRI, camera=ASTRICam),
 TelescopeDescription(type=MST, name=MST, optics=MST, camera=FlashCam)]
[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_test_large.simtel.gz")
source = EventSource(data, allowed_tels=[1,2,3,4],  max_events=10) # 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)
        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 x y r phi length length_uncertainty width width_uncertainty psi skewness kurtosis
0 487.140499 -0.972742 0.384292 1.045900 158.442970 0.242785 0.005490 0.046837 1.497981e-03 63.170863 0.324822 1.996282
1 24233.714660 -0.573750 0.723683 0.923529 128.408085 0.214088 0.001583 0.082672 6.635284e-04 -75.979331 1.663701 6.302296
2 80.082324 0.564543 0.712082 0.908719 51.592385 0.045534 0.003388 0.015214 9.861667e-04 64.285034 0.414317 2.773099
3 101.234468 0.340860 0.896085 0.958725 69.173821 0.051080 0.002661 0.020644 1.047590e-03 81.540893 0.206285 2.098711
4 76.273210 -0.457182 0.851789 0.966727 118.223829 0.044145 0.003359 0.020285 1.321683e-03 73.286384 -0.310947 2.766870
5 30.754364 -0.658798 0.686477 0.951455 133.821326 0.033088 0.003372 0.000001 3.615390e-08 -19.104678 -0.060794 2.277771
6 147.958880 -1.010481 0.264038 1.044408 165.356059 0.108207 0.008961 0.035885 1.917601e-03 39.983540 1.330492 5.058719
7 263.484926 -1.007649 -0.031596 1.008144 -178.204001 0.061533 0.002495 0.047032 2.086092e-03 -13.703852 0.437843 2.732760
8 92.449732 -0.656211 0.144985 0.672037 167.541102 0.413698 0.043318 0.029936 1.439745e-03 15.975587 1.970612 5.054535
9 108.459053 -1.109201 0.179422 1.123619 170.811523 0.100173 0.005348 0.019330 9.515764e-04 48.966843 0.802975 2.236590
10 192.653661 0.687058 -0.033318 0.687866 -2.776301 0.260307 0.031718 0.044900 2.247933e-03 -34.426418 -3.327180 12.441415
11 246.094206 0.766179 -0.546651 0.941200 -35.506985 0.108631 0.004146 0.032701 1.462323e-03 77.818626 0.106322 2.433676
12 189.636198 0.418447 -0.307282 0.519153 -36.291270 0.112610 0.005283 0.024952 1.102765e-03 27.429674 0.428147 2.669737
13 17.267195 -0.664943 0.627072 0.913985 136.678923 0.023731 0.001892 0.000000 NaN 40.893720 -0.662687 1.439154
14 54.082850 -0.200106 1.015972 1.035491 101.142350 0.035135 0.002116 0.012037 1.696142e-03 -69.567000 -0.347361 1.784459
15 107.240586 -0.537786 0.798414 0.962642 123.962933 0.130395 0.009234 0.028054 1.478532e-03 -40.199983 1.022891 3.151161
16 184.539509 -0.641409 0.844656 1.060589 127.212031 0.087639 0.006696 0.067054 2.713469e-03 6.058281 1.737274 5.308467
17 62.723886 0.043613 -0.431527 0.433725 -84.228892 0.092635 0.002665 0.011940 1.516465e-03 -75.427859 0.009396 1.207598
18 492.546910 -0.009988 -0.402209 0.402333 -91.422563 0.231939 0.004799 0.034155 1.319447e-03 -78.585111 -0.059526 1.843314
19 151.799579 0.223207 -0.504964 0.552095 -66.153343 0.129388 0.006049 0.024831 1.371782e-03 88.796737 0.677150 2.326918
20 105.300175 -0.184234 -0.821752 0.842151 -102.636589 0.478566 0.047833 0.097434 6.273209e-03 -50.775228 -2.041619 5.207802
[55]:
_ = hillas.hist(figsize=(8,8))
../_images/tutorials_ctapipe_handson_66_0.png

If you do this yourself, loop over more events to get better statistics

[ ]: