<jupyter_start><jupyter_text>Tailcut interactive notebook [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/jeremiedecock/pywi-cta-notebooks/master?filepath=tuto_3a_interactive_tailcut_cleaning_with_bokeh.ipynb)  TODO- [ ] **improve the colormap and/or use the radius to distinguish background to signal**- [ ] add a colorbar- [ ] add the possibility to get images id to plot from a (file automatically made from others notebooks)  Import required modules and set some variables This notebook requires PyWI and PyWI-CTA for the I/O and for the signal processing. It also requires Bokeh to display images (as a much faster alternative to Matplotlib).<jupyter_code>import numpy as np

#from bokeh.plotting import figure, output_notebook, show    # For fix images

import bokeh

from bokeh.io import push_notebook, output_notebook, show    # For animations
from bokeh.plotting import figure, ColumnDataSource          # For animations

from bokeh.models import LogColorMapper, LogTicker, ContinuousColorMapper, ContinuousTicker, ColorBar
from bokeh.models import HoverTool

from bokeh.models.annotations import Title

import matplotlib as mpl
import math

import os

import ipywidgets
from ipywidgets import interact

import astropy.units as u

import pywicta
import pywicta.data
from pywicta.io import geometry_converter
from pywicta.io.images import image_generator
from pywicta.io.images import plot_ctapipe_image
from pywicta.io.images import plot_hillas_parameters_on_axes
from pywicta.io.images import get_mars_like_default_integrator_config
from pywicta.image.hillas_parameters import get_hillas_parameters
from pywicta.denoising.rejection_criteria import CTAMarsCriteria

from pywicta.denoising.tailcut import Tailcut
import warnings
warnings.filterwarnings('ignore')<jupyter_output><empty_output><jupyter_text>Make the dataset The next cells define the list of images to use in this notebook.Images can be fetched from Fits files or from Simtel files.Fits files are much lighter and much faster to process than Simtel files but they are specific to PyWI-CTA and thus you first have to generate them from Simtel files using [the following script](). Also, contrary to Simtel files, a Fits files contains only one "image" (i.e. an unique event viewed from one unique telescope).The others things to configure are:1. `CAM_ID` to define the camera to use:  ASTRICam, CHEC, DigiCam, FlashCam, NectarCam or LSTCam.2. `IMG_ID_LIST` to fetch specific images e.g. set `IMG_ID_LIST = [ "run104_2081200_1", "run105_1244901_4"]` if you only want the image of the event `2081200` from telescope `1` in run `104` and the image of the event `1244901` from telescope `4` in run `105`.Alternatively, one can set `TEL_FILTER_LIST` and/or `EVENT_FILTER_LIST` to set the list of desired telescopes (e.g. `TEL_FILTER_LIST = [1, 3, 4]`) and the list of desired events (e.g. `EVENT_FILTER_LIST = [2081200, 1244901]`).<jupyter_code>#CAM_ID = "ASTRICam"
#CAM_ID = "CHEC"
#CAM_ID = "DigiCam"
#CAM_ID = "FlashCam"
#CAM_ID = "NectarCam"
CAM_ID = "LSTCam"

TEL_FILTER_LIST = []
EVENT_FILTER_LIST = []

PATHS = pywicta.data.lst()
print(PATHS)<jupyter_output><empty_output><jupyter_text>Define preselection cuts `rejection_criteria` defines a function to apply a preselection cut.<jupyter_code>#rejection_criteria = lambda ref_image: not 200 < np.nansum(ref_image) < 250
#rejection_criteria = CTAMarsCriteria(cam_id=CAM_ID)
rejection_criteria = None<jupyter_output><empty_output><jupyter_text>Get images<jupyter_code>integrator_config_dict = get_mars_like_default_integrator_config(CAM_ID)

image_dict = {"run{}_{}_{}_{}GeV".format(int(image.meta['run_id']),
                                         int(image.meta['event_id']),
                                         int(image.meta['tel_id']),
                                         int(float(image.meta['mc_energy']) * 1000)): image
              for image
              in image_generator(PATHS,
                                 cam_filter_list=[CAM_ID],
                                 tel_filter_list=TEL_FILTER_LIST,
                                 ev_filter_list=EVENT_FILTER_LIST,
                                 ctapipe_format=False,
                                 mc_rejection_criteria=rejection_criteria,
                                 **integrator_config_dict)}<jupyter_output><empty_output><jupyter_text>Plot<jupyter_code>palette_1 = bokeh.palettes.Viridis256
palette_2 = bokeh.palettes.Viridis256

palette_2[0] = "#eeeeee"
color_mapper = bokeh.models.mappers.LinearColorMapper(palette=palette_1,
                                                      low=0,
                                                      high=0)
def plot_hillas_parameters(geom,
                           image=None,
                           hillas_params=None,
                           hillas_implementation=2):
    """Plot the shower ellipse and direction on bokeh."""
    
    centroid = (0, 0)
    angle = 0.
    length = 0.
    width = 0.
    lines_xs = [[0, 0], [0, 0]]
    lines_ys = [[0, 0], [0, 0]]

    try:
        if hillas_params is None and image is not None:
            hillas_params = get_hillas_parameters(geom, image, implementation=hillas_implementation)

            centroid = (hillas_params.x.value, hillas_params.y.value)
            length = hillas_params.length.value
            width = hillas_params.width.value
            angle = hillas_params.psi.to(u.rad).value

            #print("centroid:", centroid)
            #print("length:",   length)
            #print("width:",    width)
            #print("angle:",    angle)

            p0_x = centroid[0]
            p0_y = centroid[1]

            p1_x = p0_x + math.cos(angle)
            p1_y = p0_y + math.sin(angle)

            p2_x = p0_x + math.cos(angle + math.pi)
            p2_y = p0_y + math.sin(angle + math.pi)

            lines_xs = [[p1_x, p2_x], [0, p0_x]]
            lines_ys = [[p1_y, p2_y], [0, p0_y]]

    except HillasParameterizationError as err:
        print(err)
    
    return centroid, angle, length, width, lines_xs, lines_ys
tailcut = Tailcut()
geom1d = geometry_converter.get_geom1d(CAM_ID)

if geom1d.pix_type == 'hexagonal':
    radius = math.sqrt(geom1d.pix_area.value[0]/(2. * math.sqrt(3.))) # assuming an hexagon (see:https://fr.wikipedia.org/wiki/Hexagone#Calcul_de_l'aire)
elif geom1d.pix_type == 'rectangular':
    radius = math.sqrt(geom1d.pix_area.value[0]) / 2.
else:
    raise NotImplementedError("Unknown camera type {}".format(geom1d.pix_type))
    

hover = HoverTool(
            tooltips=[
                ("PE", "@pe"),
            ]
        )
TOOLS = "crosshair,pan,wheel_zoom,box_zoom,reset,tap,previewsave,box_select,poly_select,lasso_select"

CAM_SIZE = max(geom1d.pix_x.value.max(), geom1d.pix_y.value.max()) * 1.05

FIGURE_SIZE = 600
#fig = figure(plot_width=FIGURE_SIZE, plot_height=FIGURE_SIZE, tools=TOOLS)
fig = figure(plot_width=FIGURE_SIZE,
             plot_height=FIGURE_SIZE,
             x_range=(-CAM_SIZE, CAM_SIZE),
             y_range=(-CAM_SIZE, CAM_SIZE))

output_notebook()

#colors = data_to_colors(img)

title = Title()
title.text = "-"
fig.title = title

## add a circle renderer with a size, color, and alpha
#circles = fig.circle(geom1d.pix_x.value,
#                     geom1d.pix_y.value,
#                     #size=5,              # The size (diameter) values for the markers **in screen space units** (i.e. aspect changes with figure size or zoom).
#                     radius=radius,        # The radius values for circle markers (**in "data space" units**, by default).
#                     fill_color=["#ffffff" for pix in geom1d.pix_x],
#                     line_color=["#000000" for pix in geom1d.pix_x],
#                     alpha=1.)

source = ColumnDataSource(
             data=dict(
                 x=geom1d.pix_x.value,
                 y=geom1d.pix_y.value,
                 fill_color=["#ffffff" for pix in geom1d.pix_x],
                 line_color=["#000000" for pix in geom1d.pix_x],
                 pe=[0. for pix in geom1d.pix_x],
             )
         )

circles = fig.circle("x",
                     "y",
                     radius=radius,        # The radius values for circle markers (**in "data space" units**, by default).
                     fill_color=bokeh.transform.transform('pe', color_mapper), # "fill_color",
                     line_color=bokeh.transform.transform('pe', color_mapper), # "line_color",
                     source=source)

shower_lines = fig.multi_line(xs=[[0, 0],[0, 0]],
                              ys=[[0, 0],[0, 0]],
                              color=["red", "green"],
                              line_width=2,
                              alpha=0.75)

shower_ellipse = fig.ellipse(x=0.,
                             y=0.,
                             width=0.,
                             height=0.,
                             angle=0.,
                             color="red",
                             alpha=0.5)

fig.add_tools(hover)

# show the results
handle = show(fig, notebook_handle=True)

@interact(plot_type=["raw", "ref.", "cleaned"], image_key=image_dict.keys(), picture_thresh=10, boundary_thresh=5, pixels_clusters=["off", "scipy", "mars"])
def compute_hillas_and_display(plot_type, image_key, picture_thresh=10, boundary_thresh=5, pixels_clusters="off"):
    image = image_dict[image_key]
    cam_id = image.meta['cam_id']
    image.meta['npe'] = np.nansum(image.reference_image)
    show_background = True
    
    if plot_type == "raw":

        img_1d = geometry_converter.image_2d_to_1d(image.input_image, cam_id)

        shower_centroid, shower_angle, shower_length, shower_width, shower_lines_xs, shower_lines_ys = plot_hillas_parameters(geom1d)
        show_background = True
            
    elif plot_type == "ref.":

        img_1d = geometry_converter.image_2d_to_1d(image.reference_image, cam_id)
        
        #reference_image_mask = np.full(reference_image_1d.shape, False)
        #reference_image_mask[reference_image_1d > 0] = True
        
        shower_centroid, shower_angle, shower_length, shower_width, shower_lines_xs, shower_lines_ys = plot_hillas_parameters(geom1d, image=img_1d)
        show_background = False

    elif plot_type == "cleaned":

        calibrated_image = image.input_image
        tailcut_cleaned_img = tailcut.clean_image(calibrated_image,
                                                  high_threshold=picture_thresh,
                                                  low_threshold=boundary_thresh,
                                                  pixels_clusters_filtering=pixels_clusters,
                                                  cam_id=cam_id)
        show_background = False

        if np.nanmax(tailcut_cleaned_img) == 0:

            img_1d = None
            shower_centroid, shower_angle, shower_length, shower_width, shower_lines_xs, shower_lines_ys = plot_hillas_parameters(geom1d)
            
        else:
            
            img_1d = geometry_converter.image_2d_to_1d(tailcut_cleaned_img, cam_id)
            
            shower_centroid, shower_angle, shower_length, shower_width, shower_lines_xs, shower_lines_ys = plot_hillas_parameters(geom1d, image=img_1d)

            #print("cen_x:...", hillas_params.cen_x)
            #print("cen_y:...", hillas_params.cen_y)
            #print("length:..", hillas_params.length)
            #print("width:...", hillas_params.width)
            #print("size:....", hillas_params.size, "PE")
            #print("NPE: ....", np.sum(tailcut_cleaned_img), "PE")
            #print("psi:.....", hillas_params.psi)
            #print("miss:....", hillas_params.miss)
            #print("phi:.....", hillas_params.phi)
            #print("r:.......", hillas_params.r)
            #print("kurtosis:", hillas_params.kurtosis)
            #print("skewness:", hillas_params.skewness)
    
    # Update the plot
    if img_1d is None:
        circles.data_source.data['fill_color'] = ["#ffffff" for color in geom1d.pix_x]
        circles.data_source.data['line_color'] = ["#aaaaaa" for color in geom1d.pix_x]
    else:
        circles.data_source.data['pe'] = img_1d
        color_mapper.low = img_1d.min()
        color_mapper.high = img_1d.max()
        
        if show_background:
            color_mapper.palette = palette_1
        else:
            color_mapper.palette = palette_2
        
    shower_lines.data_source.data['xs'] = shower_lines_xs
    shower_lines.data_source.data['ys'] = shower_lines_ys
    
    shower_ellipse.glyph.x = shower_centroid[0]
    shower_ellipse.glyph.y = shower_centroid[1]
    shower_ellipse.glyph.angle = shower_angle
    shower_ellipse.glyph.width = shower_length
    shower_ellipse.glyph.height =  shower_width
        
    #print(image.meta)
    
    if plot_type == "raw":
        plot_desc = "calibrated image"
    elif plot_type == "ref.":
        plot_desc = "MC image"
    elif plot_type == "cleaned":
        plot_desc = "Tailcut clean"
    else:
        plot_desc = plot_type
        
    title.text = "Run{} Ev{} Tel{} {:0.3f}TeV {}NPE {}".format(image.meta['run_id'],
                                                               image.meta['event_id'],
                                                               image.meta['tel_id'],
                                                               image.meta['mc_energy'],
                                                               image.meta['npe'],
                                                               plot_desc)
    fig.title = title
    
    push_notebook(handle=handle)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10239332503842494, 'fraction_numerical': 0.016321452096904047, 'mean_word_length': 2.2988894253983583, 'pattern_counts': {'":': 6, '<': 23, '<?xml version=': 0, '>': 22, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 34, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '228398', 'license': 'MIT', 'n_tokens_mistral': 4416, 'n_tokens_neox': 4233, 'n_words': 898, 'path': 'tuto_3a_interactive_tailcut_cleaning_with_bokeh.ipynb', 'repo_name': 'jeremiedecock/pywi-cta-notebooks', 'size': '13663'}