Notes:
The idea here is to define a benchmark to optimise cleaning independently of any reconstruction that would come after.
This to avoid optimising the cleaning as a function of the whole reconstruction as:
This benchmark uses the the ground thruth image in photo-electron from MC simulations by computing the distance between the cleaned image and the ground truth as a function of cleaning method/parameters and finding the minimum of this distance (average on many events).
This also allow to study the cleaning as a function of event info (such as energy, signal amplitude... )
Of course, this supposes that the calibration has been previously optimised.
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:90% !important; }</style>"))
import ctapipe
print(ctapipe.__version__)
0.10.5
from ctapipe.io import EventSource
from ctapipe.utils import datasets
from ctapipe.calib import CameraCalibrator
from ctapipe.image import tailcuts_clean, dilate
from ctapipe.visualization import CameraDisplay
from ctapipe.instrument import CameraGeometry
import matplotlib.pyplot as plt
import numpy as np
from tqdm import tqdm
from scipy.stats import poisson
import os
from ctapipe.utils import get_dataset_path
from ctapipe.io import read_table
from ctapipe.instrument import SubarrayDescription # for working with CTA instruments
from astropy.table import join
from ctapipe.utils.download import download_file_cached
import copy
import astropy.units as u
import tables
from astropy.table import Table, vstack
%matplotlib inline
ls ../../prepared_data/
ls: cannot access '../../prepared_data/': No such file or directory
remote_url = "http://cccta-dataserver.in2p3.fr/data/Prod5_Paranal_North_20deg_ctapipe_v0.10.5_DL1/"
filename = "gamma_20deg_0deg_run107___cta-prod5-paranal_desert-2147m-Paranal-dark_cone10_merged.DL1.h5"
filename = download_file_cached(filename, default_url=remote_url)
filename
PosixPath('/home/runner/.cache/ctapipe/cccta-dataserver.in2p3.fr/data/Prod5_Paranal_North_20deg_ctapipe_v0.10.5_DL1/gamma_20deg_0deg_run107___cta-prod5-paranal_desert-2147m-Paranal-dark_cone10_merged.DL1.h5')
subarray = SubarrayDescription.from_hdf(filename)
subarray.info()
subarray.peek()
Subarray : MonteCarloArray
Num Tels : 180
Footprint: 4.92 km2
Type Count Tel IDs
----------------- ----- ---------------
MST_MST_FlashCam 28 5-29,125-127
LST_LST_LSTCam 4 1-4
SST_ASTRI_CHEC 120 30-99,131-180
MST_MST_NectarCam 28 100-124,128-130
subarray.tel_ids
array([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,
14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26,
27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39,
40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52,
53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65,
66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78,
79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91,
92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104,
105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117,
118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130,
131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143,
144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156,
157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169,
170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180], dtype=int16)
telescope_types = subarray.telescope_types
def read_images_from_telescope_type(filename, telescope_type):
images_tables = []
for tel_id in subarray.get_tel_ids_for_type(telescope_type):
images = read_table(filename, f"/dl1/event/telescope/images/tel_{tel_id:03d}")
sim_images = read_table(filename, f"/simulation/event/telescope/images/tel_{tel_id:03d}")
images_tables.append(join(images, sim_images, keys=['event_id', 'tel_id', 'obs_id']))
return vstack(images_tables)
subarray.to_table()
| tel_id | pos_x | pos_y | pos_z | name | type | num_mirrors | camera_type | tel_description |
|---|---|---|---|---|---|---|---|---|
| m | m | m | ||||||
| int16 | float64 | float64 | float64 | str5 | str3 | int64 | str9 | str17 |
| 1 | -20.64299964904785 | -64.81700134277344 | 34.29999923706055 | LST | LST | 1 | LSTCam | LST_LST_LSTCam |
| 2 | 79.9939956665039 | -0.7680000066757202 | 29.399999618530273 | LST | LST | 1 | LSTCam | LST_LST_LSTCam |
| 3 | -19.395999908447266 | 65.19999694824219 | 31.0 | LST | LST | 1 | LSTCam | LST_LST_LSTCam |
| 4 | -120.03299713134766 | 1.1510000228881836 | 33.099998474121094 | LST | LST | 1 | LSTCam | LST_LST_LSTCam |
| 5 | -0.017000000923871994 | -0.0010000000474974513 | 24.350000381469727 | MST | MST | 1 | FlashCam | MST_MST_FlashCam |
| 6 | -1.468000054359436 | -151.2209930419922 | 31.0 | MST | MST | 1 | FlashCam | MST_MST_FlashCam |
| 7 | -3.1379997730255127 | -325.2449951171875 | 39.0 | MST | MST | 1 | FlashCam | MST_MST_FlashCam |
| 8 | 1.4339998960494995 | 151.22000122070312 | 25.0 | MST | MST | 1 | FlashCam | MST_MST_FlashCam |
| 9 | 3.1039998531341553 | 325.2430114746094 | 23.5 | MST | MST | 1 | FlashCam | MST_MST_FlashCam |
| ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 171 | 260.0 | 920.0 | 45.0 | ASTRI | SST | 2 | CHEC | SST_ASTRI_CHEC |
| 172 | 260.0 | -920.0 | 65.0 | ASTRI | SST | 2 | CHEC | SST_ASTRI_CHEC |
| 173 | -500.0 | 815.0 | 15.0 | ASTRI | SST | 2 | CHEC | SST_ASTRI_CHEC |
| 174 | -500.0 | -815.0 | 75.0 | ASTRI | SST | 2 | CHEC | SST_ASTRI_CHEC |
| 175 | 500.0 | 815.0 | 45.0 | ASTRI | SST | 2 | CHEC | SST_ASTRI_CHEC |
| 176 | 500.0 | -815.0 | 53.0 | ASTRI | SST | 2 | CHEC | SST_ASTRI_CHEC |
| 177 | -810.0 | 655.0 | 12.0 | ASTRI | SST | 2 | CHEC | SST_ASTRI_CHEC |
| 178 | -810.0 | -655.0 | 68.0 | ASTRI | SST | 2 | CHEC | SST_ASTRI_CHEC |
| 179 | 810.0 | 655.0 | 20.0 | ASTRI | SST | 2 | CHEC | SST_ASTRI_CHEC |
| 180 | 810.0 | -655.0 | 41.0 | ASTRI | SST | 2 | CHEC | SST_ASTRI_CHEC |
telescope_types
[TelescopeDescription(type=MST, name=MST, optics=MST, camera=FlashCam), TelescopeDescription(type=LST, name=LST, optics=LST, camera=LSTCam), TelescopeDescription(type=SST, name=ASTRI, optics=ASTRI, camera=CHEC), TelescopeDescription(type=MST, name=MST, optics=MST, camera=NectarCam)]
tel_type = telescope_types[0]
geometry = subarray.tels[subarray.get_tel_ids_for_type(tel_type)[0]].camera.geometry
image_table = read_images_from_telescope_type(filename, tel_type)
def residuals_after_cleaning(cleaned_image, true_image):
return (cleaned_image-true_image)
import copy
def add_residuals_to_table(image_table):
cleaned_images = copy.deepcopy(image_table['image'])
cleaned_images[~image_table['image_mask']]=0
image_table['residuals'] = residuals_after_cleaning(cleaned_images, image_table['true_image'])
image_table['accuracy'] = np.linalg.norm(image_table['residuals'], axis=1)
add_residuals_to_table(image_table)
image_table[:3]
| obs_id | event_id | tel_id | image [1764] | peak_time [1764] | image_mask [1764] | true_image [1764] | residuals [1764] | accuracy |
|---|---|---|---|---|---|---|---|---|
| int32 | int64 | int16 | float32 | float32 | bool | int32 | float64 | float64 |
| 860 | 21604 | 5 | -0.6 .. 0.6 | 85.24 .. 64.97 | False .. False | 0 .. 0 | 0.0 .. 0.0 | 57.00771538250507 |
| 860 | 21609 | 5 | 3.2 .. -0.7 | 53.77 .. 54.86 | False .. False | 0 .. 0 | 0.0 .. 0.0 | 50.943988823544494 |
| 860 | 28506 | 5 | 5.6 .. -5.7 | 65.41 .. 48.08 | False .. False | 0 .. 0 | 0.0 .. 0.0 | 12.671622882409407 |
plt.hist(image_table['residuals'].ravel(), log=True, bins=100, range=(-20, 20));
print("residuals mean: ", np.mean(np.abs(image_table['residuals'])))
residuals mean: 0.2250397098023214
plt.hist(image_table['accuracy'], bins=100, range=(0, 100));
from matplotlib.colors import Normalize
def display_row(geometry, image_table, row_index=0):
fig, axes = plt.subplots(1, 3, figsize=(20,5))
row = image_table[row_index]
display = CameraDisplay(geometry, row['image'], ax=axes[0])
display.add_colorbar()
display.highlight_pixels(row['image_mask'], color='red', alpha=0.3)
display.axes.set_title('image')
display = CameraDisplay(geometry, row['true_image'], ax=axes[1])
display.add_colorbar()
display.axes.set_title('true_image')
if 'residuals' in row.colnames:
display = CameraDisplay(geometry, row['residuals'], ax=axes[2], cmap='RdBu')
max_pe = np.max(np.abs(row['residuals']))
display.add_colorbar()
display.set_limits_minmax(-max_pe, max_pe)
display.axes.set_title('residuals')
return axes
image_table[4]
| obs_id | event_id | tel_id | image [1764] | peak_time [1764] | image_mask [1764] | true_image [1764] | residuals [1764] | accuracy |
|---|---|---|---|---|---|---|---|---|
| int32 | int64 | int16 | float32 | float32 | bool | int32 | float64 | float64 |
| 4447 | 77409 | 5 | -4.2 .. 1.4 | 20.07 .. 52.82 | False .. False | 0 .. 0 | 0.0 .. 0.0 | 15.574017021914528 |
display_row(geometry, image_table, 4);
def thresholds_grid(image_table, pt_array=np.linspace(3, 12, 10)):
acc = []
picture_threshold = []
boundary_threshold = []
for pt in pt_array:
for bt in np.linspace(0, pt, len(pt_array)):
picture_threshold.append(pt)
boundary_threshold.append(bt)
tailcut_opt = dict(picture_thresh=pt, boundary_thresh=bt)
image_mask = [tailcuts_clean(geometry, image, **tailcut_opt) for image in image_table['image']]
image_table['image_mask'] = image_mask
add_residuals_to_table(image_table)
# acc.append(np.mean(image_table['accuracy']))
acc.append((np.linalg.norm(image_table['residuals'].ravel(), ord=2))/image_table['residuals'].ravel().shape[0])
return np.array(picture_threshold), np.array(boundary_threshold), np.array(acc)
def best_thresholds(picture_threshold, boundary_threshold, accuracy):
"""
return picture_threshold, boundary_threshold
"""
return picture_threshold[np.argmin(acc)], boundary_threshold[np.argmin(acc)]
def plot_threshold_heatmap(picture_threshold, boundary_threshold, accuracy):
fig, ax = plt.subplots(figsize=(10, 6))
im = ax.tricontourf(picture_threshold,boundary_threshold,acc)
cbar = plt.colorbar(im)
ax.set_xlabel('picture threshold')
ax.set_ylabel('boundary threhsold')
cbar.set_label('accuracy')
ax.axis('equal')
return ax
for tel_type in telescope_types:
print(f"---- {tel_type} ----")
geometry = subarray.tels[subarray.get_tel_ids_for_type(tel_type)[0]].camera.geometry
image_table = read_images_from_telescope_type(filename, tel_type)[:1000]
add_residuals_to_table(image_table)
print("Example:")
display_row(geometry, image_table, 0)
plt.show()
pt, bt, acc = thresholds_grid(image_table, np.linspace(4, 20, 10))
print(f"best threshold for {tel_type}: {best_thresholds(pt, bt, acc)}")
plot_threshold_heatmap(pt, bt, acc)
plt.show()
---- MST_MST_FlashCam ---- Example:
best threshold for MST_MST_FlashCam: (12.88888888888889, 2.864197530864198)
---- LST_LST_LSTCam ---- Example:
best threshold for LST_LST_LSTCam: (5.777777777777778, 1.9259259259259258)
---- SST_ASTRI_CHEC ---- Example:
best threshold for SST_ASTRI_CHEC: (4.0, 0.0)
---- MST_MST_NectarCam ---- Example:
best threshold for MST_MST_NectarCam: (5.777777777777778, 0.0)
