2018 LST Bootcamp walkthrough

ctapipe

LST Analysis Bootcamp

Padova, 26.11.2018

Maximilian Nöthe (@maxnoe) & Kai A. Brügge (@mackaiver)

[1]:
import matplotlib.pyplot as plt
import numpy as np

%matplotlib inline
[2]:
plt.rcParams['figure.figsize'] = (12, 8)
plt.rcParams['font.size'] = 14
plt.rcParams['figure.figsize']
[2]:
[12.0, 8.0]

Table of Contents

General Information

Design

  • DL0 → DL3 analysis

  • Currently some R0 → DL0 code to be able to analyze simtel files

  • ctapipe is built upon the Scientific Python Stack, core dependencies are

    • numpy

    • scipy

    • astropy

Developement

  • ctapipe is developed as Open Source Software (Currently under MIT License) at https://github.com/cta-observatory/ctapipe

  • We use the “Github-Workflow”:

    • Few people (e.g. @kosack, @mackaiver) have write access to the main repository

    • Contributors fork the main repository and work on branches

    • Pull Requests are merged after Code Review and automatic execution of the test suite

  • Early developement stage ⇒ backwards-incompatible API changes might and will happen

  • Many open design questions ⇒ Core Developer Meeting in the second week of December in Dortmund

What’s there?

  • Reading simtel simulation files

  • Simple calibration, cleaning and feature extraction functions

  • Camera and Array plotting

  • Coordinate frames and transformations

  • Stereo-reconstruction using line intersections

What’s still missing?

  • Easy to use IO of analysis results to standard data formats (e.g. FITS, hdf5)

  • Easy to use “analysis builder”

  • A “Standard Analysis”

  • Good integration with machine learning techniques

  • IRF calculation

  • Defining APIs for IO, instrument description access etc.

  • Most code only tested on HESSIO simulations

  • Documentation, e.g. formal definitions of coordinate frames

What can you do?

  • Report issues

    • Hard to get started? Tell us where you are stuck

    • Tell user stories

    • Missing features

  • Start contributing

    • ctapipe needs more workpower

    • Implement new reconstruction features

A simple hillas analysis

Reading in simtel files

[3]:
from ctapipe.io import EventSource
from ctapipe.utils.datasets import get_dataset_path

input_url = get_dataset_path('gamma_test_large.simtel.gz')

# EventSource() automatically detects what kind of file we are giving it,
# if already supported by ctapipe
source = EventSource(input_url, max_events=49)

print(type(source))
/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))
<class 'ctapipe.io.simteleventsource.SimTelEventSource'>
[4]:
for event in source:
    print('Id: {}, E = {:1.3f}, Telescopes: {}'.format(event.count, event.simulation.shower.energy, len(event.r0.tel)))
Id: 0, E = 0.571 TeV, Telescopes: 4
Id: 1, E = 1.864 TeV, Telescopes: 9
Id: 2, E = 1.864 TeV, Telescopes: 4
Id: 3, E = 1.864 TeV, Telescopes: 17
Id: 4, E = 1.864 TeV, Telescopes: 2
Id: 5, E = 0.464 TeV, Telescopes: 7
Id: 6, E = 0.017 TeV, Telescopes: 2
Id: 7, E = 76.426 TeV, Telescopes: 4
Id: 8, E = 76.426 TeV, Telescopes: 16
Id: 9, E = 76.426 TeV, Telescopes: 3
Id: 10, E = 0.267 TeV, Telescopes: 5
Id: 11, E = 0.010 TeV, Telescopes: 2
Id: 12, E = 1.407 TeV, Telescopes: 3
Id: 13, E = 1.407 TeV, Telescopes: 3
Id: 14, E = 0.121 TeV, Telescopes: 10
Id: 15, E = 0.032 TeV, Telescopes: 4
Id: 16, E = 0.073 TeV, Telescopes: 2
Id: 17, E = 0.129 TeV, Telescopes: 3
Id: 18, E = 0.129 TeV, Telescopes: 2
Id: 19, E = 10.220 TeV, Telescopes: 2
Id: 20, E = 10.220 TeV, Telescopes: 4
Id: 21, E = 10.220 TeV, Telescopes: 4
Id: 22, E = 0.055 TeV, Telescopes: 2
Id: 23, E = 0.100 TeV, Telescopes: 3
Id: 24, E = 0.079 TeV, Telescopes: 5
Id: 25, E = 0.127 TeV, Telescopes: 7
Id: 26, E = 3.923 TeV, Telescopes: 4
Id: 27, E = 3.923 TeV, Telescopes: 2
Id: 28, E = 3.923 TeV, Telescopes: 4
Id: 29, E = 0.147 TeV, Telescopes: 6
Id: 30, E = 0.974 TeV, Telescopes: 6
Id: 31, E = 0.974 TeV, Telescopes: 3
Id: 32, E = 0.474 TeV, Telescopes: 4
Id: 33, E = 0.474 TeV, Telescopes: 4
Id: 34, E = 0.110 TeV, Telescopes: 2
Id: 35, E = 0.103 TeV, Telescopes: 2
Id: 36, E = 0.120 TeV, Telescopes: 4
Id: 37, E = 0.078 TeV, Telescopes: 4
Id: 38, E = 0.079 TeV, Telescopes: 2
Id: 39, E = 0.592 TeV, Telescopes: 4
Id: 40, E = 0.066 TeV, Telescopes: 6
Id: 41, E = 0.126 TeV, Telescopes: 3
Id: 42, E = 0.126 TeV, Telescopes: 2
Id: 43, E = 0.129 TeV, Telescopes: 2
Id: 44, E = 0.017 TeV, Telescopes: 3
Id: 45, E = 0.053 TeV, Telescopes: 2
Id: 46, E = 0.009 TeV, Telescopes: 2
Id: 47, E = 0.248 TeV, Telescopes: 4
Id: 48, E = 0.248 TeV, Telescopes: 11

Each event is a DataContainer holding several Fields of data, which can be containers or just numbers. Let’s look a one event:

[5]:
event
[5]:
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)
[6]:
source.subarray.camera_types
[6]:
[CameraDescription(camera_name=ASTRICam, geometry=ASTRICam, readout=ASTRICam),
 CameraDescription(camera_name=FlashCam, geometry=FlashCam, readout=FlashCam),
 CameraDescription(camera_name=LSTCam, geometry=LSTCam, readout=LSTCam)]
[7]:
len(event.r0.tel), len(event.r1.tel)
[7]:
(11, 11)

Data calibration

The CameraCalibrator calibrates the event (obtaining the dl1 images).

[8]:
from ctapipe.calib import CameraCalibrator

calibrator = CameraCalibrator(subarray=source.subarray)
[9]:
calibrator(event)

Event displays

Let’s use ctapipe’s plotting facilities to plot the telescope images

[10]:
event.dl1.tel.keys()
[10]:
dict_keys([1, 2, 3, 4, 5, 7, 9, 11, 13, 16, 19])
[11]:
tel_id = 4
[12]:
geometry = source.subarray.tel[tel_id].camera.geometry
dl1 = event.dl1.tel[tel_id]

geometry, dl1
[12]:
(CameraGeometry(camera_name='LSTCam', pix_type=<PixelShape.HEXAGON: 'hexagon'>, npix=1855, cam_rot=0.0 deg, pix_rot=100.89299992867878 deg),
 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)
[13]:
dl1.image
[13]:
array([ 1.8977832 , -1.7927673 ,  2.2859168 , ..., -1.0049672 ,
        1.1161382 ,  0.36395848], dtype=float32)
[14]:
from ctapipe.visualization import CameraDisplay

display = CameraDisplay(geometry)

# right now, there might be one image per gain channel.
# This will change as soon as
display.image = dl1.image
display.add_colorbar()
../_images/tutorials_lst_analysis_bootcamp_2018_27_0.png

Image Cleaning

[15]:
from ctapipe.image.cleaning import tailcuts_clean
[16]:
# unoptimized cleaning levels, copied from
# https://github.com/tudo-astroparticlephysics/cta_preprocessing
cleaning_level = {
    'ASTRICam': (5, 7, 2),  # (5, 10)?
    'LSTCam': (3.5, 7.5, 2),  # ?? (3, 6) for Abelardo...
    'FlashCam': (4, 8, 2),  # there is some scaling missing?
}
[17]:
boundary, picture, min_neighbors = cleaning_level[geometry.camera_name]

clean = tailcuts_clean(
    geometry,
    dl1.image,
    boundary_thresh=boundary,
    picture_thresh=picture,
    min_number_picture_neighbors=min_neighbors
)
[18]:
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 5))

d1 = CameraDisplay(geometry, ax=ax1)
d2 = CameraDisplay(geometry, ax=ax2)

ax1.set_title('Image')
d1.image = dl1.image
d1.add_colorbar(ax=ax1)

ax2.set_title('Pulse Time')
d2.image = dl1.peak_time - np.average(dl1.peak_time, weights=dl1.image)
d2.cmap = 'RdBu_r'
d2.add_colorbar(ax=ax2)
d2.set_limits_minmax(-20,20)

d1.highlight_pixels(clean, color='red', linewidth=1)
../_images/tutorials_lst_analysis_bootcamp_2018_32_0.png

Image Parameters

[19]:
from ctapipe.image import hillas_parameters, leakage_parameters, concentration_parameters
from ctapipe.image import timing_parameters
from ctapipe.image import number_of_islands
from ctapipe.image import camera_to_shower_coordinates
[20]:
hillas = hillas_parameters(geometry[clean], dl1.image[clean])

print(hillas)
{'intensity': 293.30156207084656,
 'kurtosis': 2.1156918832117624,
 'length': <Quantity 0.09458967 m>,
 'length_uncertainty': <Quantity 0.00291695 m>,
 'phi': <Angle 0.71440377 rad>,
 'psi': <Angle -0.71608044 rad>,
 'r': <Quantity 0.67893267 m>,
 'skewness': -0.15071029797454388,
 'width': <Quantity 0.03181493 m>,
 'width_uncertainty': <Quantity 0.00096868 m>,
 'x': <Quantity 0.51292278 m>,
 'y': <Quantity 0.44481435 m>}
[21]:
display = CameraDisplay(geometry)

# set "unclean" pixels to 0
cleaned = dl1.image.copy()
cleaned[~clean] = 0.0

display.image = cleaned
display.add_colorbar()

display.overlay_moments(hillas, color='xkcd:red')
../_images/tutorials_lst_analysis_bootcamp_2018_36_0.png
[22]:
timing = timing_parameters(
    geometry,
    dl1.image,
    dl1.peak_time,
    hillas,
    clean
)

print(timing)
{'deviation': 0.5507465160185104,
 'intercept': 30.871464798739247,
 'slope': <Quantity -5.81771647 1 / m>}
[23]:
long, trans = camera_to_shower_coordinates(
    geometry.pix_x, geometry.pix_y,hillas.x, hillas.y, hillas.psi
)

plt.plot(long[clean], dl1.peak_time[clean], 'o')
plt.plot(long[clean], timing.slope * long[clean] + timing.intercept)
[23]:
[<matplotlib.lines.Line2D at 0x7fb1df4ee850>]
../_images/tutorials_lst_analysis_bootcamp_2018_38_1.png
[24]:
l = leakage_parameters(geometry, dl1.image, clean)
print(l)
{'intensity_width_1': 0.0,
 'intensity_width_2': 0.0,
 'pixels_width_1': 0.0,
 'pixels_width_2': 0.0}
[25]:
disp = CameraDisplay(geometry)
disp.image = dl1.image
disp.highlight_pixels(geometry.get_border_pixel_mask(1), linewidth=2, color='xkcd:red')
../_images/tutorials_lst_analysis_bootcamp_2018_40_0.png
[26]:
n_islands, island_id = number_of_islands(geometry, clean)

print(n_islands)
1
[27]:
conc = concentration_parameters(geometry, dl1.image, hillas)
print(conc)
{'cog': 0.21367737999682757,
 'core': 0.34470898401576955,
 'pixel': 0.10290955973463788}

Putting it all together / Stereo reconstruction

[28]:
import astropy.units as u
from astropy.coordinates import SkyCoord, AltAz

from ctapipe.io import EventSource
from ctapipe.utils.datasets import get_dataset_path

from ctapipe.calib import CameraCalibrator

from ctapipe.image import tailcuts_clean, number_of_islands
from ctapipe.image import hillas_parameters, leakage_parameters, concentration_parameters
from ctapipe.image import timing_parameters

from ctapipe.reco import HillasReconstructor

from ctapipe.io import HDF5TableWriter

from copy import deepcopy
import tempfile

# unoptimized cleaning levels, copied from
# https://github.com/tudo-astroparticlephysics/cta_preprocessing
cleaning_level = {
    'ASTRICam': (5, 7, 2),  # (5, 10)?
    'LSTCam': (3.5, 7.5, 2),  # ?? (3, 6) for Abelardo...
    'FlashCam': (4, 8, 2),  # there is some scaling missing?
}


input_url = get_dataset_path('gamma_test_large.simtel.gz')
source = EventSource(input_url)

calibrator = CameraCalibrator(subarray=source.subarray)

horizon_frame = AltAz()

reco = HillasReconstructor()

f = tempfile.NamedTemporaryFile(suffix='.hdf5')

with HDF5TableWriter(f.name, mode='w', group_name='events') as writer:

    for event in source:
        print('Id: {}, E = {:1.3f}, Telescopes: {}'.format(event.count, event.simulation.shower.energy, len(event.r0.tel)))

        calibrator(event)

        # mapping of telescope_id to parameters for stereo reconstruction
        hillas_containers = {}
        telescope_pointings = {}
        time_gradients = {}

        for telescope_id, dl1 in event.dl1.tel.items():
            geometry = source.subarray.tels[telescope_id].camera.geometry
            image = dl1.image
            peak_time = dl1.peak_time

            boundary, picture, min_neighbors = cleaning_level[geometry.camera_name]

            clean = tailcuts_clean(
                geometry,
                image,
                boundary_thresh=boundary,
                picture_thresh=picture,
                min_number_picture_neighbors=min_neighbors
            )

            # require more than five pixels after cleaning in each telescope
            if clean.sum() < 5:
                continue

            hillas_c = hillas_parameters(geometry[clean], image[clean])
            leakage_c = leakage_parameters(geometry, image, clean)
            n_islands, island_ids = number_of_islands(geometry, clean)

            # remove events with high leakage
            if leakage_c.intensity_width_2 > 0.2:
                continue

            timing_c = timing_parameters(geometry, image, peak_time, hillas_c, clean)

            hillas_containers[telescope_id] = hillas_c

            # ssts have no timing in prod3b, so we'll use the skewness
            time_gradients[telescope_id] = timing_c.slope.value if geometry.camera_name != 'ASTRICam' else hillas_c.skewness

            # this makes sure, that we get an arrow in the array plow for each telescope
            # might have the wrong direction though
            if abs(time_gradients[telescope_id]) < 0.2:
                time_gradients[telescope_id] = 1.0

            telescope_pointings[telescope_id] = SkyCoord(
                alt=event.pointing.tel[telescope_id].altitude,
                az=event.pointing.tel[telescope_id].azimuth,
                frame=horizon_frame
            )

        # the array pointing is needed for the creation of the TiltedFrame to perform the
        # impact point reconstruction
        array_pointing = SkyCoord(
            az=event.pointing.array_azimuth,
            alt=event.pointing.array_altitude,
            frame=horizon_frame
        )

        if len(hillas_containers) > 1:
            stereo = reco.predict(
                hillas_containers, source.subarray,
                array_pointing, telescope_pointings
            )

            writer.write('reconstructed', stereo)
            writer.write('true', event.simulation.shower)

            print('  Alt: {:.2f}°'.format(stereo.alt.deg))
            print('  Az: {:.2f}°'.format(stereo.az.deg))
            print('  Hmax: {:.0f}'.format(stereo.h_max))
            print('  CoreX: {:.1f}'.format(stereo.core_x))
            print('  CoreY: {:.1f}'.format(stereo.core_y))

        # save a nice event for plotting later
        if event.count == 3:
            plotting_event = deepcopy(event)
            plotting_hillas = hillas_containers
            plotting_timing = time_gradients
            plotting_stereo = stereo
Id: 0, E = 0.571 TeV, Telescopes: 4
Column tel_ids of container ReconstructedShowerContainer in table reconstructed not writable, skipping
  Alt: 66.37°
  Az: -7.85°
  Hmax: 8963 m
  CoreX: 114.8 m
  CoreY: -414.8 m
Id: 1, E = 1.864 TeV, Telescopes: 9
  Alt: 68.25°
  Az: -353.46°
  Hmax: 6798 m
  CoreX: -48.5 m
  CoreY: -394.6 m
Id: 2, E = 1.864 TeV, Telescopes: 4
  Alt: 68.49°
  Az: -354.35°
  Hmax: 7050 m
  CoreX: -419.5 m
  CoreY: -574.8 m
Id: 3, E = 1.864 TeV, Telescopes: 17
  Alt: 68.47°
  Az: -353.41°
  Hmax: 6655 m
  CoreX: 42.2 m
  CoreY: 67.4 m
Id: 4, E = 1.864 TeV, Telescopes: 2
  Alt: 68.64°
  Az: -353.31°
  Hmax: 6396 m
  CoreX: -1097.3 m
  CoreY: 474.8 m
Id: 5, E = 0.464 TeV, Telescopes: 7
  Alt: 71.92°
  Az: -1.07°
  Hmax: 8684 m
  CoreX: -66.1 m
  CoreY: 289.6 m
Id: 6, E = 0.017 TeV, Telescopes: 2
  Alt: 70.52°
  Az: -358.18°
  Hmax: 10667 m
  CoreX: -159.0 m
  CoreY: -232.3 m
Id: 7, E = 76.426 TeV, Telescopes: 4
  Alt: 73.99°
  Az: -354.24°
  Hmax: 9817 m
  CoreX: 11.9 m
  CoreY: -1121.4 m
Id: 8, E = 76.426 TeV, Telescopes: 16
  Alt: 75.26°
  Az: -353.94°
  Hmax: 5815 m
  CoreX: 137.1 m
  CoreY: 449.8 m
Id: 9, E = 76.426 TeV, Telescopes: 3
  Alt: 75.00°
  Az: -354.03°
  Hmax: 7485 m
  CoreX: 913.0 m
  CoreY: -995.4 m
Id: 10, E = 0.267 TeV, Telescopes: 5
Id: 11, E = 0.010 TeV, Telescopes: 2
Id: 12, E = 1.407 TeV, Telescopes: 3
Id: 13, E = 1.407 TeV, Telescopes: 3
Id: 14, E = 0.121 TeV, Telescopes: 10
  Alt: 67.41°
  Az: -0.73°
  Hmax: 13016 m
  CoreX: -163.9 m
  CoreY: -137.2 m
Id: 15, E = 0.032 TeV, Telescopes: 4
  Alt: 71.90°
  Az: -0.47°
  Hmax: 7851 m
  CoreX: 21.2 m
  CoreY: 65.7 m
Id: 16, E = 0.073 TeV, Telescopes: 2
Id: 17, E = 0.129 TeV, Telescopes: 3
  Alt: 70.31°
  Az: -353.32°
  Hmax: 4184 m
  CoreX: 186.6 m
  CoreY: -130.5 m
Id: 18, E = 0.129 TeV, Telescopes: 2
Id: 19, E = 10.220 TeV, Telescopes: 2
Id: 20, E = 10.220 TeV, Telescopes: 4
Id: 21, E = 10.220 TeV, Telescopes: 4
Id: 22, E = 0.055 TeV, Telescopes: 2
  Alt: 73.02°
  Az: -1.67°
  Hmax: 7389 m
  CoreX: -335.0 m
  CoreY: 8.7 m
Id: 23, E = 0.100 TeV, Telescopes: 3
  Alt: 68.86°
  Az: -350.50°
  Hmax: 7785 m
  CoreX: -123.7 m
  CoreY: -148.0 m
Id: 24, E = 0.079 TeV, Telescopes: 5
  Alt: 68.34°
  Az: -353.47°
  Hmax: 10671 m
  CoreX: -113.7 m
  CoreY: 91.9 m
Id: 25, E = 0.127 TeV, Telescopes: 7
  Alt: 70.11°
  Az: -4.55°
  Hmax: 8605 m
  CoreX: 23.1 m
  CoreY: -198.6 m
Id: 26, E = 3.923 TeV, Telescopes: 4
  Alt: 70.37°
  Az: -345.42°
  Hmax: 5987 m
  CoreX: 373.0 m
  CoreY: -194.6 m
Id: 27, E = 3.923 TeV, Telescopes: 2
Id: 28, E = 3.923 TeV, Telescopes: 4
Id: 29, E = 0.147 TeV, Telescopes: 6
  Alt: 70.10°
  Az: -3.79°
  Hmax: 5555 m
  CoreX: 131.2 m
  CoreY: 224.3 m
Id: 30, E = 0.974 TeV, Telescopes: 6
  Alt: 73.41°
  Az: -5.79°
  Hmax: 8489 m
  CoreX: 460.3 m
  CoreY: 328.2 m
Id: 31, E = 0.974 TeV, Telescopes: 3
  Alt: 73.11°
  Az: -6.73°
  Hmax: 7629 m
  CoreX: -84.8 m
  CoreY: 554.9 m
Id: 32, E = 0.474 TeV, Telescopes: 4
  Alt: 70.06°
  Az: -10.68°
  Hmax: 7234 m
  CoreX: -123.8 m
  CoreY: -478.5 m
Id: 33, E = 0.474 TeV, Telescopes: 4
  Alt: 69.83°
  Az: -13.61°
  Hmax: 7167 m
  CoreX: 448.8 m
  CoreY: -303.2 m
Id: 34, E = 0.110 TeV, Telescopes: 2
  Alt: 71.91°
  Az: -358.11°
  Hmax: 8842 m
  CoreX: -221.9 m
  CoreY: -7.7 m
Id: 35, E = 0.103 TeV, Telescopes: 2
Id: 36, E = 0.120 TeV, Telescopes: 4
  Alt: 68.60°
  Az: -8.66°
  Hmax: 10768 m
  CoreX: 118.9 m
  CoreY: -206.3 m
Id: 37, E = 0.078 TeV, Telescopes: 4
  Alt: 68.58°
  Az: -355.25°
  Hmax: 9472 m
  CoreX: -170.6 m
  CoreY: -115.5 m
Id: 38, E = 0.079 TeV, Telescopes: 2
Id: 39, E = 0.592 TeV, Telescopes: 4
Id: 40, E = 0.066 TeV, Telescopes: 6
  Alt: 68.82°
  Az: -6.06°
  Hmax: 10677 m
  CoreX: 67.5 m
  CoreY: -14.0 m
Id: 41, E = 0.126 TeV, Telescopes: 3
  Alt: 71.57°
  Az: -355.17°
  Hmax: 7758 m
  CoreX: -210.1 m
  CoreY: -308.9 m
Id: 42, E = 0.126 TeV, Telescopes: 2
Id: 43, E = 0.129 TeV, Telescopes: 2
Id: 44, E = 0.017 TeV, Telescopes: 3
Id: 45, E = 0.053 TeV, Telescopes: 2
  Alt: 70.77°
  Az: -355.33°
  Hmax: 10482 m
  CoreX: 461.5 m
  CoreY: -97.6 m
Id: 46, E = 0.009 TeV, Telescopes: 2
Id: 47, E = 0.248 TeV, Telescopes: 4
  Alt: 71.94°
  Az: -0.06°
  Hmax: 10458 m
  CoreX: 467.2 m
  CoreY: 96.6 m
Id: 48, E = 0.248 TeV, Telescopes: 11
  Alt: 71.99°
  Az: -359.96°
  Hmax: 10187 m
  CoreX: 8.3 m
  CoreY: -147.2 m
Id: 49, E = 0.248 TeV, Telescopes: 4
  Alt: 71.83°
  Az: -359.98°
  Hmax: 8856 m
  CoreX: -450.0 m
  CoreY: 149.1 m
Id: 50, E = 5.327 TeV, Telescopes: 2
  Alt: 66.71°
  Az: -349.98°
  Hmax: 6980 m
  CoreX: -1059.0 m
  CoreY: 871.1 m
Id: 51, E = 5.327 TeV, Telescopes: 18
  Alt: 67.08°
  Az: -351.11°
  Hmax: 6522 m
  CoreX: -340.3 m
  CoreY: 145.6 m
Id: 52, E = 5.327 TeV, Telescopes: 2
Id: 53, E = 0.095 TeV, Telescopes: 2
Id: 54, E = 0.937 TeV, Telescopes: 2
Id: 55, E = 0.937 TeV, Telescopes: 2
  Alt: 69.82°
  Az: -6.77°
  Hmax: 9255 m
  CoreX: -650.0 m
  CoreY: -704.9 m
Id: 56, E = 0.937 TeV, Telescopes: 4
  Alt: 69.79°
  Az: -7.33°
  Hmax: 9587 m
  CoreX: -601.1 m
  CoreY: 133.3 m
Id: 57, E = 0.937 TeV, Telescopes: 2
  Alt: 70.15°
  Az: -6.80°
  Hmax: 9136 m
  CoreX: 414.7 m
  CoreY: 866.7 m
Id: 58, E = 0.161 TeV, Telescopes: 2
Id: 59, E = 0.407 TeV, Telescopes: 14
  Alt: 69.64°
  Az: -8.06°
  Hmax: 9542 m
  CoreX: -76.2 m
  CoreY: -88.8 m
Id: 60, E = 0.037 TeV, Telescopes: 3
  Alt: 67.22°
  Az: -358.15°
  Hmax: 13977 m
  CoreX: -437.6 m
  CoreY: 51.7 m
Id: 61, E = 0.020 TeV, Telescopes: 4
  Alt: 70.55°
  Az: -0.70°
  Hmax: 9237 m
  CoreX: 1.8 m
  CoreY: 123.1 m
Id: 62, E = 2.364 TeV, Telescopes: 4
  Alt: 65.91°
  Az: -359.40°
  Hmax: 7748 m
  CoreX: 573.7 m
  CoreY: 730.2 m
Id: 63, E = 0.710 TeV, Telescopes: 8
  Alt: 71.35°
  Az: -349.97°
  Hmax: 10259 m
  CoreX: 302.3 m
  CoreY: 222.5 m
Id: 64, E = 0.710 TeV, Telescopes: 2
  Alt: 71.12°
  Az: -348.22°
  Hmax: 10447 m
  CoreX: 677.4 m
  CoreY: -573.6 m
Id: 65, E = 1.158 TeV, Telescopes: 2
Id: 66, E = 0.055 TeV, Telescopes: 5
  Alt: 69.85°
  Az: -356.16°
  Hmax: 8673 m
  CoreX: -13.5 m
  CoreY: 160.4 m
Id: 67, E = 0.032 TeV, Telescopes: 2
Id: 68, E = 0.232 TeV, Telescopes: 3
  Alt: 68.56°
  Az: -5.28°
  Hmax: 10128 m
  CoreX: 484.4 m
  CoreY: -442.1 m
Id: 69, E = 0.232 TeV, Telescopes: 2
Id: 70, E = 22.299 TeV, Telescopes: 5
  Alt: 66.80°
  Az: -358.30°
  Hmax: 6775 m
  CoreX: 139.1 m
  CoreY: -855.5 m
Id: 71, E = 22.299 TeV, Telescopes: 12
  Alt: 63.74°
  Az: -356.20°
  Hmax: 7888 m
  CoreX: -170.7 m
  CoreY: 271.0 m
Id: 72, E = 0.118 TeV, Telescopes: 3
  Alt: 67.25°
  Az: -4.93°
  Hmax: 9275 m
  CoreX: 186.7 m
  CoreY: 290.8 m
Id: 73, E = 0.118 TeV, Telescopes: 3
  Alt: 67.35°
  Az: -4.97°
  Hmax: 8865 m
  CoreX: 205.6 m
  CoreY: -187.7 m
Id: 74, E = 0.118 TeV, Telescopes: 5
  Alt: 67.02°
  Az: -5.38°
  Hmax: 9031 m
  CoreX: -61.6 m
  CoreY: -336.3 m
Id: 75, E = 0.101 TeV, Telescopes: 2
Id: 76, E = 0.065 TeV, Telescopes: 2
Id: 77, E = 3.581 TeV, Telescopes: 2
  Alt: 71.05°
  Az: -350.78°
  Hmax: 9400 m
  CoreX: -382.5 m
  CoreY: -1149.8 m
Id: 78, E = 3.581 TeV, Telescopes: 5
  Alt: 71.11°
  Az: -350.72°
  Hmax: 8273 m
  CoreX: 188.4 m
  CoreY: -667.7 m
Id: 79, E = 0.768 TeV, Telescopes: 5
Id: 80, E = 0.290 TeV, Telescopes: 2
Id: 81, E = 0.060 TeV, Telescopes: 2
Id: 82, E = 2.967 TeV, Telescopes: 2
Id: 83, E = 2.967 TeV, Telescopes: 2
Id: 84, E = 0.142 TeV, Telescopes: 4
  Alt: 69.42°
  Az: -358.78°
  Hmax: 10535 m
  CoreX: -149.3 m
  CoreY: -367.3 m
Id: 85, E = 0.142 TeV, Telescopes: 2
Id: 86, E = 0.435 TeV, Telescopes: 2
Id: 87, E = 0.435 TeV, Telescopes: 2
Id: 88, E = 0.339 TeV, Telescopes: 4
Id: 89, E = 1.851 TeV, Telescopes: 4
  Alt: 69.58°
  Az: -2.82°
  Hmax: 7769 m
  CoreX: -879.3 m
  CoreY: 105.7 m
Id: 90, E = 1.851 TeV, Telescopes: 6
  Alt: 69.47°
  Az: -2.52°
  Hmax: 7770 m
  CoreX: -711.7 m
  CoreY: -250.7 m
Id: 91, E = 0.189 TeV, Telescopes: 11
  Alt: 70.61°
  Az: -3.62°
  Hmax: 8683 m
  CoreX: -51.5 m
  CoreY: -92.9 m
Id: 92, E = 0.037 TeV, Telescopes: 3
Id: 93, E = 1.620 TeV, Telescopes: 4
  Alt: 69.32°
  Az: -7.08°
  Hmax: 7291 m
  CoreX: -867.7 m
  CoreY: -195.8 m
Id: 94, E = 1.620 TeV, Telescopes: 2
  Alt: 69.31°
  Az: -7.58°
  Hmax: 7869 m
  CoreX: -346.0 m
  CoreY: -1083.3 m
Id: 95, E = 1.620 TeV, Telescopes: 3
  Alt: 69.07°
  Az: -7.32°
  Hmax: 8700 m
  CoreX: 771.3 m
  CoreY: -465.3 m
Id: 96, E = 1.620 TeV, Telescopes: 3
  Alt: 69.27°
  Az: -7.34°
  Hmax: 8186 m
  CoreX: 93.0 m
  CoreY: 864.1 m
Id: 97, E = 1.620 TeV, Telescopes: 7
  Alt: 69.19°
  Az: -6.66°
  Hmax: 7739 m
  CoreX: 317.6 m
  CoreY: -596.2 m
Id: 98, E = 0.112 TeV, Telescopes: 3
  Alt: 70.67°
  Az: -358.12°
  Hmax: 9389 m
  CoreX: -22.3 m
  CoreY: -344.7 m
Id: 99, E = 0.322 TeV, Telescopes: 2
  Alt: 69.40°
  Az: -359.91°
  Hmax: 6798 m
  CoreX: -538.2 m
  CoreY: 62.1 m
Id: 100, E = 0.322 TeV, Telescopes: 10
  Alt: 69.43°
  Az: -0.75°
  Hmax: 7703 m
  CoreX: -178.3 m
  CoreY: -223.8 m
Id: 101, E = 6.104 TeV, Telescopes: 4
Id: 102, E = 0.245 TeV, Telescopes: 3
Id: 103, E = 0.636 TeV, Telescopes: 2
Id: 104, E = 0.636 TeV, Telescopes: 4
Id: 105, E = 1.824 TeV, Telescopes: 6
  Alt: 70.25°
  Az: -0.72°
  Hmax: 7238 m
  CoreX: -612.5 m
  CoreY: -411.4 m
Id: 106, E = 1.824 TeV, Telescopes: 2
Id: 107, E = 1.824 TeV, Telescopes: 3
  Alt: 70.19°
  Az: -0.60°
  Hmax: 6885 m
  CoreX: 822.3 m
  CoreY: -124.4 m
Id: 108, E = 0.538 TeV, Telescopes: 5
  Alt: 69.24°
  Az: -348.65°
  Hmax: 8039 m
  CoreX: -368.7 m
  CoreY: 93.4 m
Id: 109, E = 0.538 TeV, Telescopes: 2
[29]:
from astropy.coordinates.angle_utilities import angular_separation
import pandas as pd


df_rec = pd.read_hdf(f.name, key='events/reconstructed')
df_true = pd.read_hdf(f.name, key='events/true')


theta = angular_separation(
    df_rec.az.values * u.deg, df_rec.alt.values * u.deg,
    df_true.az.values * u.deg, df_true.alt.values * u.deg,
)

plt.hist(theta.to(u.deg).value**2, bins=25, range=[0, 0.3])
plt.xlabel(r'$\theta² / deg²$')
None
../_images/tutorials_lst_analysis_bootcamp_2018_45_0.png

ArrayDisplay

[30]:
from ctapipe.visualization import ArrayDisplay


angle_offset = plotting_event.pointing.array_azimuth


disp = ArrayDisplay(source.subarray)

disp.set_vector_hillas(
    plotting_hillas,
    time_gradient=plotting_timing,
    angle_offset=angle_offset,
    length=500
)

plt.scatter(
    plotting_event.simulation.shower.core_x, plotting_event.simulation.shower.core_y,
    s=200, c='k', marker='x', label='True Impact',
)
plt.scatter(
    plotting_stereo.core_x, plotting_stereo.core_y,
    s=200, c='r', marker='x', label='Estimated Impact',
)

plt.legend()
plt.xlim(-400, 400)
plt.ylim(-400, 400)
[30]:
(-400.0, 400.0)
../_images/tutorials_lst_analysis_bootcamp_2018_47_1.png

LST Mono with output

  • Let’s use the HDF5TableWriter to save the dl1 Hillas parameter data to an hdf5 file

  • This is not ideal yet and one of the major points to be discussed in two weeks

[31]:
from ctapipe.io import HDF5TableWriter
from ctapipe.core.container import Container, Field

input_url = get_dataset_path('gamma_test_large.simtel.gz')

source = EventSource(
    input_url,
    allowed_tels=[1, 2, 3, 4], # only use the first LST
)

calibrator = CameraCalibrator(subarray=source.subarray)

class EventInfo(Container):
    event_id = Field('event_id')
    obs_id =  Field('obs_id')
    telescope_id = Field('telescope_id')


with HDF5TableWriter(filename='hillas.h5', group_name='dl1', mode='w') as writer:

    for event in source:
        print('Id: {}, E = {:1.3f}, Telescopes: {}'.format(event.count, event.simulation.shower.energy, len(event.r0.tel)))

        calibrator(event)

        for telescope_id, dl1 in event.dl1.tel.items():

            geometry = source.subarray.tels[telescope_id].camera.geometry
            image = dl1.image
            peak_time = dl1.peak_time

            boundary, picture, min_neighbors = cleaning_level[geometry.camera_name]

            clean = tailcuts_clean(
                geometry,
                image,
                boundary_thresh=boundary,
                picture_thresh=picture,
                min_number_picture_neighbors=min_neighbors
            )

            if clean.sum() < 5:
                continue

            event_info = EventInfo(event_id=event.index.event_id, obs_id=event.index.obs_id, telescope_id=telescope_id)
            hillas_c = hillas_parameters(geometry[clean], image[clean])
            leakage_c = leakage_parameters(geometry, image, clean)
            timing_c = timing_parameters(geometry, image, peak_time, hillas_c, clean)

            writer.write('events', [event_info, event.simulation.shower, hillas_c, leakage_c, timing_c])

Id: 0, E = 1.864 TeV, Telescopes: 2
Id: 1, E = 0.017 TeV, Telescopes: 2
Id: 2, E = 0.010 TeV, Telescopes: 2
Id: 3, E = 0.121 TeV, Telescopes: 4
Id: 4, E = 0.032 TeV, Telescopes: 3
Id: 5, E = 0.073 TeV, Telescopes: 1
Id: 6, E = 0.129 TeV, Telescopes: 1
Id: 7, E = 0.079 TeV, Telescopes: 2
Id: 8, E = 0.127 TeV, Telescopes: 3
Id: 9, E = 0.147 TeV, Telescopes: 1
Id: 10, E = 0.103 TeV, Telescopes: 1
Id: 11, E = 0.078 TeV, Telescopes: 1
Id: 12, E = 0.592 TeV, Telescopes: 1
Id: 13, E = 0.066 TeV, Telescopes: 2
Id: 14, E = 0.017 TeV, Telescopes: 3
Id: 15, E = 0.009 TeV, Telescopes: 2
Id: 16, E = 0.248 TeV, Telescopes: 4
Id: 17, E = 5.327 TeV, Telescopes: 3
Id: 18, E = 0.407 TeV, Telescopes: 4
Id: 19, E = 0.037 TeV, Telescopes: 3
Id: 20, E = 0.020 TeV, Telescopes: 4
Id: 21, E = 0.710 TeV, Telescopes: 3
Id: 22, E = 0.055 TeV, Telescopes: 2
Id: 23, E = 0.032 TeV, Telescopes: 2
Id: 24, E = 0.118 TeV, Telescopes: 1
Id: 25, E = 0.142 TeV, Telescopes: 1
Id: 26, E = 0.189 TeV, Telescopes: 4
Id: 27, E = 0.037 TeV, Telescopes: 3
Id: 28, E = 0.112 TeV, Telescopes: 1
Id: 29, E = 0.322 TeV, Telescopes: 4
[32]:
import pandas as pd

df = pd.read_hdf('hillas.h5', key='dl1/events')
df.set_index(['obs_id', 'event_id', 'telescope_id'], inplace=True)
df.head()
[32]:
energy alt az core_x core_y h_first_int x_max shower_primary_id intensity x ... psi skewness kurtosis pixels_width_1 pixels_width_2 intensity_width_1 intensity_width_2 slope intercept deviation
obs_id event_id telescope_id
7514 31012 2 1.863750 68.478978 6.384091 51.675117 69.654037 18171.726562 379.538452 0 474.993227 -0.963622 ... 61.248556 0.231937 2.031285 0.011321 0.017790 0.653862 0.862198 -5.434742 32.110888 0.729585
3 1.863750 68.478978 6.384091 51.675117 69.654037 18171.726562 379.538452 0 24090.585138 -0.574563 ... -75.514828 1.560784 5.846910 0.009164 0.018868 0.114064 0.272069 1.909613 33.687022 0.598446
90914 2 0.016544 70.676229 3.107417 -113.833435 -142.343872 15776.052734 184.545456 0 77.619497 0.564774 ... 53.199432 -0.120016 2.473936 0.000000 0.000000 0.000000 0.000000 6.165023 30.539039 0.464106
4 0.016544 70.676229 3.107417 -113.833435 -142.343872 15776.052734 184.545456 0 110.688694 0.341674 ... 82.347736 0.050001 2.086657 0.000000 0.000000 0.000000 0.000000 -2.437174 30.185566 0.492260
153614 1 0.010250 68.949967 2.666187 -51.531204 -106.752579 19824.572266 241.111115 0 76.273210 -0.457182 ... 73.286384 -0.310947 2.766870 0.000000 0.000000 0.000000 0.000000 -2.750752 30.057657 0.185710

5 rows × 27 columns

[33]:
plt.scatter(np.log10(df.energy), np.log10(df.intensity))
plt.xlabel('log10(E / TeV)')
plt.ylabel('log10(intensity)')
None
../_images/tutorials_lst_analysis_bootcamp_2018_51_0.png

Isn’t python slow?

  • Many of you might have heard: “Python is slow”.

  • That’s trueish.

  • All python objects are classes living on the heap, event integers.

  • Looping over lots of “primitives” is quite slow compared to other languages.

⇒ Vectorize as much as possible using numpy
⇒ Use existing interfaces to fast C / C++ / Fortran code
⇒ Optimize using cython or numba

But: “Premature Optimization is the root of all evil” — Donald Knuth

So profile to find exactly what is slow.

Why use python then?

  • Python works very well as glue for libraries of all kinds of languages

  • Python has a rich ecosystem for data science, physics, algorithms, astronomy

Example: Number of Islands

Find all groups of pixels, that survived the cleaning

[34]:
from ctapipe.image import toymodel
from ctapipe.instrument import CameraGeometry


geometry = CameraGeometry.from_name('LSTCam')

Let’s create a toy images with several islands;

[35]:
np.random.seed(42)

image = np.zeros(geometry.n_pixels)


for i in range(9):

    model = toymodel.Gaussian(
        x=np.random.uniform(-0.8, 0.8) * u.m,
        y=np.random.uniform(-0.8, 0.8) * u.m,
        width=np.random.uniform(0.05, 0.075) * u.m,
        length=np.random.uniform(0.1, 0.15) * u.m,
        psi=np.random.uniform(0, 2 * np.pi) * u.rad,
    )

    new_image, sig, bg = model.generate_image(
        geometry,
        intensity=np.random.uniform(1000, 3000),
        nsb_level_pe=5
    )
    image += new_image
[36]:
clean = tailcuts_clean(geometry, image, picture_thresh=10, boundary_thresh=5, min_number_picture_neighbors=2)
[37]:
disp = CameraDisplay(geometry)
disp.image = image
disp.highlight_pixels(clean, color='xkcd:red', linewidth=1.5)
disp.add_colorbar()
../_images/tutorials_lst_analysis_bootcamp_2018_57_0.png
[38]:
def num_islands_python(camera, clean):
    ''' A breadth first search to find connected islands of neighboring pixels in the cleaning set'''

    # the camera geometry has a [n_pixel, n_pixel] boolean array
    # that is True where two pixels are neighbors
    neighbors = camera.neighbor_matrix

    island_ids = np.zeros(camera.n_pixels)
    current_island = 0

    # a set to remember which pixels we already visited
    visited = set()

    # go only through the pixels, that survived cleaning
    for pix_id in np.where(clean)[0]:
        if pix_id not in visited:
            # remember that we already checked this pixel
            visited.add(pix_id)

            # if we land in the outer loop again, we found a new island
            current_island += 1
            island_ids[pix_id] = current_island

            # now check all neighbors of the current pixel recursively
            to_check = set(np.where(neighbors[pix_id] & clean)[0])
            while to_check:
                pix_id = to_check.pop()

                if pix_id not in visited:
                    visited.add(pix_id)
                    island_ids[pix_id] = current_island

                    to_check.update(np.where(neighbors[pix_id] & clean)[0])

    n_islands = current_island
    return n_islands, island_ids
[39]:
n_islands, island_ids = num_islands_python(geometry, clean)
[40]:
from matplotlib.colors import ListedColormap

cmap = plt.get_cmap('Paired')
cmap = ListedColormap(cmap.colors[:n_islands])
cmap.set_under('k')

disp = CameraDisplay(geometry)
disp.image = island_ids
disp.cmap = cmap
disp.set_limits_minmax(0.5, n_islands + 0.5)
disp.add_colorbar()
../_images/tutorials_lst_analysis_bootcamp_2018_60_0.png
[41]:
%timeit num_islands_python(geometry, clean)
3.16 ms ± 80.5 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
[42]:
from scipy.sparse.csgraph import connected_components

def num_islands_scipy(geometry, clean):
    neighbors = geometry.neighbor_matrix_sparse

    clean_neighbors = neighbors[clean][:, clean]
    num_islands, labels = connected_components(clean_neighbors, directed=False)

    island_ids = np.zeros(geometry.n_pixels)
    island_ids[clean] = labels + 1

    return num_islands, island_ids
[43]:
n_islands_s, island_ids_s = num_islands_scipy(geometry, clean)
[44]:
disp = CameraDisplay(geometry)
disp.image = island_ids_s
disp.cmap = cmap
disp.set_limits_minmax(0.5, n_islands_s + 0.5)
disp.add_colorbar()
../_images/tutorials_lst_analysis_bootcamp_2018_64_0.png
[45]:
%timeit num_islands_scipy(geometry, clean)
693 µs ± 13.9 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

A lot less code, and a factor 3 speed improvement