<jupyter_start><jupyter_text>Direction recontruction (DL2) **Recommended datasample(s):**Datasets of fully-analyzed showers used to obtain Instrument Response Functions, which in the default pipeline workflow are called ``gamma3``, ``proton2`` and ``electron``.**Data level(s):** DL2 (shower geometry + estimated energy + estimated particle classification)**Description:**This notebook contains benchmarks for the _protopipe_ pipeline regarding the shower geometry of events which have been completely analyzed.**Requirements and steps to reproduce:**- get a TRAINING file generated using ``protopipe-DL2`` or the equivalent command from the DIRAC Grid interface- execute the notebook with ``protopipe-BENCHMARK``,``protopipe-BENCHMARK launch --config_file configs/benchmarks.yaml -n DL2/benchmarks_DL2_direction-reconstruction``To obtain the list of all available parameters add ``--help-notebook``.**Development and testing:**    As with any other part of _protopipe_ and being part of the official repository, this notebook can be further developed by any interested contributor.   The execution of this notebook is not currently automatic, it must be done locally by the user _before_ pushing a pull-request.  Please, strip the output before pushing.**TODO:*** ...  Table of contents   - [Energy-dependent offset distribution](Energy-dependent-offset-distribution)   - [Angular-resolution-as-a-function-of-telescope-multiplicity](Angular-resolution-as-a-function-of-telescope-multiplicity)   - [Angular resolution for different containment radii and fixed signal efficiency](Angular-resolution-for-different-containment-radii-and-fixed-signal-efficiency)   - [PSF asymmetry](PSF-asymmetry)   - [True energy distributions](True-energy-distributions)  Imports<jupyter_code>import os
from pathlib import Path

%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from matplotlib.pyplot import rc
import matplotlib.style as style
from cycler import cycler

import numpy as np
import pandas as pd

import astropy.units as u
from astropy.coordinates import SkyCoord
from ctapipe.coordinates import NominalFrame
from pyirf.binning import (
    add_overflow_bins,
    create_bins_per_decade
)

from protopipe.benchmarks.utils import get_fig_size, string_to_boolean, get_fig_size
from protopipe.benchmarks.operations import compute_psf
from protopipe.benchmarks.plot import plot_psf<jupyter_output><empty_output><jupyter_text>Load data<jupyter_code># Parametrized cell
analyses_directory = None
output_directory = Path.cwd() # default output directory for plots
analysis_name = None
analysis_name_2 = None
gammas_infile_name = "DL2_tail_gamma_merged.h5"
protons_infile_name = "DL2_tail_proton_merged.h5"
electrons_infile_name = "DL2_tail_electron_merged.h5"
efficiency_cut = 0.9
export_data = False
superimpose_analysis_2 = False
use_seaborn = True
plots_scale = None
# Handle boolean variables (papermill reads them as strings)
[use_seaborn,
 export_data,
 superimpose_analysis_2] = string_to_boolean([use_seaborn,
                                              export_data,
                                              superimpose_analysis_2])
input_directory = Path(analyses_directory) / analysis_name / Path("data/DL2")
gammas = pd.read_hdf(os.path.join(input_directory, "gamma", gammas_infile_name), "/reco_events")
protons = pd.read_hdf(os.path.join(input_directory, "proton", protons_infile_name), "/reco_events")
electrons = pd.read_hdf(os.path.join(input_directory, "electron", electrons_infile_name), "/reco_events")
basic_selection_cut = (gammas["is_valid"]==True) & (gammas["NTels_reco"] >= 2)
selected_gammaness = gammas[basic_selection_cut]["gammaness"]
gammaness_cut = np.quantile(selected_gammaness, efficiency_cut)
selected_gammas = gammas[basic_selection_cut & (gammas["gammaness"] >= gammaness_cut)]
#selected_gammas = gammas[(gammas["is_valid"]==True) & (gammas["NTels_reco"] >= 2) & (gammas["gammaness"] >= 0.90)]
# First we check if a _plots_ folder exists already.  
# If not, we create it.
plots_folder = Path(output_directory) / "plots"
plots_folder.mkdir(parents=True, exist_ok=True)

# Next we check if a _data_ folder exists already.  
# If not, we create it.
data_folder = Path(output_directory) / "data"
data_folder.mkdir(parents=True, exist_ok=True)

input_directory_data_2 = Path(analyses_directory) / analysis_name_2/ "benchmarks_results/TRAINING"
# Plot aesthetics settings

scale = matplotlib_settings["scale"] if plots_scale is None else float(plots_scale)
style.use(matplotlib_settings["style"])
cmap = matplotlib_settings["cmap"]

if matplotlib_settings["style"] == "seaborn-colorblind":
    
    colors_order = ['#0072B2', '#D55E00', '#F0E442', '#009E73', '#CC79A7', '#56B4E9']
    rc('axes', prop_cycle=cycler(color=colors_order))

if use_seaborn:
    import seaborn as sns

    sns.set_theme(context=seaborn_settings["theme"]["context"] if "context" in seaborn_settings["theme"] else "talk",
                  style=seaborn_settings["theme"]["style"] if "style" in seaborn_settings["theme"] else "whitegrid",
                  palette=seaborn_settings["theme"]["palette"] if "palette" in seaborn_settings["theme"] else None,
                  font=seaborn_settings["theme"]["font"] if "font" in seaborn_settings["theme"] else "Fira Sans",
                  font_scale=seaborn_settings["theme"]["font_scale"] if "font_scale" in seaborn_settings["theme"] else 1.0,
                  color_codes=seaborn_settings["theme"]["color_codes"] if "color_codes" in seaborn_settings["theme"] else True
                  )
    
    sns.set_style(seaborn_settings["theme"]["style"], rc=seaborn_settings["rc_style"])
    sns.set_context(seaborn_settings["theme"]["context"],
                    font_scale=seaborn_settings["theme"]["font_scale"] if "font_scale" in seaborn_settings["theme"] else 1.0)<jupyter_output><empty_output><jupyter_text>Benchmarks Here we use events with the following cuts:- valid reconstructed events- at least 2 reconstructed images, regardless of the camera (on top of any other hardware trigger)- gammaness > 0.75 (mostly a conservative choice)<jupyter_code>min_true_energy = 0.006
max_true_energy = 660<jupyter_output><empty_output><jupyter_text>Energy-dependent offset distribution[back to top](Table-of-contents)<jupyter_code>n_bins = 4
true_energy_bin_edges = np.logspace(np.log10(min_true_energy),
                                    np.log10(max_true_energy), n_bins + 1)

plt.figure(figsize=get_fig_size(ratio=4./3., scale=scale))

plt.xlabel("Offset [deg]")
plt.ylabel("Number of events")


for i in range(len(true_energy_bin_edges)-1):
    
    low_E = true_energy_bin_edges[i]
    high_E = true_energy_bin_edges[i+1]
    
    selected_events = selected_gammas[(selected_gammas["true_energy"]>low_E) & (selected_gammas["true_energy"]<high_E)]
 
    plt.hist(selected_events["offset"], 
             bins=100,
             #range = [0,10],
             label=f"{low_E:.2f} < E_true [TeV] < {high_E:.2f}",
             histtype="step",
             linewidth=2)

plt.yscale("log")
plt.legend(loc="best")
plt.grid(which="both")

plt.savefig(plots_folder / f"DL2_offsets_{analysis_name}.png")

plt.show()
min_true_energy = [0.02, 0.2, 2, 20]
max_true_energy = [0.2, 2, 20, 200]

plt.figure(figsize=(10,5))
plt.xlabel("Offset [deg]")
plt.ylabel("Number of events")

for low_E, high_E in zip(min_true_energy, max_true_energy):
    
    selected_events = selected_gammas[(selected_gammas["true_energy"]>low_E) & (selected_gammas["true_energy"]<high_E)]
    
    plt.hist(selected_events["offset"], 
             bins=100,
             range = [0,10],
             label=f"{low_E} < E_true [TeV] < {high_E}",
             histtype="step",
             linewidth=2)

plt.yscale("log")
plt.legend(loc="best")
plt.grid(which="both")

plt.savefig(plots_folder / f"DL2_offsets_{analysis_name}.png")

plt.show()
min_true_energy = [0.02, 0.2, 2, 20]
max_true_energy = [0.2, 2, 20, 200]

plt.figure(figsize=(10,5))
plt.xlabel("Offset [deg]")
plt.ylabel("Number of events")

for low_E, high_E in zip(min_true_energy, max_true_energy):
    
    selected_events = selected_gammas[(selected_gammas["true_energy"]>low_E) & (selected_gammas["true_energy"]<high_E)]
    
    plt.hist(selected_events["offset"], 
             bins=100,
             range = [0,10],
             label=f"{low_E} < E_true [TeV] < {high_E}",
             histtype="step",
             linewidth=2)

plt.yscale("log")
plt.legend(loc="best")
plt.grid(which="both")

plt.savefig(plots_folder / f"DL2_offsets_{analysis_name}.png")

plt.show()<jupyter_output><empty_output><jupyter_text>Angular resolution as a function of telescope multiplicity[back to top](Table-of-contents) Here we compare how the multiplicity influences the performance of reconstructed events with a 90% gamma efficiency within a 68% containment radius.<jupyter_code>r_containment = 68

min_true_energy = 0.003
max_true_energy = 330
n_true_energy_bins = 21
true_energy_bin_edges = np.logspace(np.log10(min_true_energy),
                                    np.log10(max_true_energy),
                                    n_true_energy_bins)
true_energy_bin_centers = 0.5 * (true_energy_bin_edges[:-1]+true_energy_bin_edges[1:])

multiplicity_cuts = ['NTels_reco == 2',
                     'NTels_reco == 3',
                     'NTels_reco == 4',
                     'NTels_reco >= 2']

fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(16, 8))
axes = axes.flatten()

for cut_idx, cut in enumerate(multiplicity_cuts):
    data_mult = selected_gammas.query(cut)
    psf, err_psf = compute_psf(data_mult, true_energy_bin_edges, 68)
    plot_psf(axes[0], true_energy_bin_centers, psf, err_psf, label=multiplicity_cuts[cut_idx])
    
    y, tmp = np.histogram(data_mult['true_energy'], bins=true_energy_bin_edges)
    weights = np.ones_like(y)
    #weights = weights / float(np.sum(y))
    yerr = np.sqrt(y) * weights
    width = np.diff(true_energy_bin_edges)
    axes[1].bar(true_energy_bin_centers, y * weights, width=width, yerr=yerr, **{'label': multiplicity_cuts[cut_idx], 'lw': 2, 'fill': False})
    axes[1].set_ylabel('Number of events')
    
for ax in axes:
    #ax.set_xlim(limit)
    ax.set_xscale('log')
    ax.legend(loc='best')
    ax.grid(which='both', visible=True)
    ax.set_xlabel('True energy [TeV]')

plt.tight_layout()

fig.savefig(plots_folder / f"DL2_PSF_{analysis_name}.png")<jupyter_output><empty_output><jupyter_text>Angular resolution for different containment radii and fixed signal efficiency[back to top](Table-of-contents) Apply fixed signal efficiency cut (requires well defined ML separator and ML train-ing)Calculate angular resolution for 68%, 80%, and 95% containment radius.<jupyter_code>scale=0.75
plt.figure(figsize=(16*scale,9*scale))

true_energy_bins = create_bins_per_decade(10**-1.9 * u.TeV, 10**2.31 * u.TeV, 5).value

gamma_efficiency = 0.9
reconstructed_gammas = gammas.query("is_valid == True")
gammaness = reconstructed_gammas["gammaness"]
gammaness_cut = np.quantile(gammaness, gamma_efficiency)
selected_events = reconstructed_gammas.query(f"gammaness > {gammaness_cut}")

ax = plt.gca()

def angular_resolution_vs_true_energy(ax, events, true_energy_bins, containment):
    
    ang_res = []
    for i in range(len(true_energy_bins)-1):
        true_energy_mask = f"true_energy > {true_energy_bins[i]} & true_energy < {true_energy_bins[i+1]}"
        selected_offsets = events.query(true_energy_mask)["offset"]
        if len(selected_offsets)==0:
            ang_res.append(np.nan)
        else:
            ang_res.append(np.quantile(selected_offsets, containment/100.))

    ax.errorbar(
        0.5 * (true_energy_bins[:-1] + true_energy_bins[1:]),
        ang_res,
        xerr=0.5 * (true_energy_bins[:-1] - true_energy_bins[1:]),
        label=f'{containment}% containment radius',
        fmt='o',
    )
    
    return ax

angular_resolution_vs_true_energy(ax, selected_events, true_energy_bins, 68)
angular_resolution_vs_true_energy(ax, selected_events, true_energy_bins, 80)
angular_resolution_vs_true_energy(ax, selected_events, true_energy_bins, 95)

plt.xlabel("True energy [TeV]")
plt.xscale("log")
plt.ylabel("Angular resolution [deg]")
plt.legend()
plt.title(f"Reconstructed gammas with {gamma_efficiency*100}% signal efficiency")
plt.grid()<jupyter_output><empty_output><jupyter_text>H_max as a function of energy for gammas and protons Fixed gamma efficiency at 90%<jupyter_code>reconstructed_gammas = gammas.query("is_valid == True")
reconstructed_protons = protons.query("is_valid == True")

plt.figure(figsize=(12,6))

mask_gammaness = f"gammaness > 0.9"

plt.subplot(1, 2, 1)

hist_opt = {"bins":[100,100],
            "range": [[0.003, 300],[1,8]],
            "norm": LogNorm(vmin=1,vmax=1.e6),
            "cmap": cmap}

plt.hist2d(reconstructed_gammas.query(mask_gammaness)["reco_energy"],
           np.log10(reconstructed_gammas.query(mask_gammaness)["h_max"]),
           **hist_opt
          )
plt.xlabel("Reconstructed energy [TeV]")
plt.ylabel("log10(H max)")
plt.colorbar()
plt.title("DL2 gammas")

plt.subplot(1, 2, 2)

plt.hist2d(reconstructed_protons.query(mask_gammaness)["reco_energy"],
           np.log10(reconstructed_protons.query(mask_gammaness)["h_max"]),
           **hist_opt
          )
plt.xlabel("Reconstructed energy [TeV]")
plt.ylabel("log10(H max)")
plt.colorbar()
plt.title("DL2 protons")
None<jupyter_output><empty_output><jupyter_text>PSF asymmetry[back to top](Table-of-contents)<jupyter_code>reco_alt = selected_gammas.reco_alt
reco_az = selected_gammas.reco_az

# right now all reco_az for a 180° deg simualtion turn out to be all around -180°
#if ~np.count_nonzero(np.sign(reco_az) + 1):
reco_az = np.abs(reco_az)

# this is needed for projecting the angle onto the sky
reco_az_corr = reco_az * np.cos(np.deg2rad(selected_gammas.reco_alt))

true_alt = selected_gammas.iloc[0].true_alt
true_az = selected_gammas.iloc[0].true_az

daz = reco_az - true_az
daz_corr = daz * np.cos(np.deg2rad(reco_alt))
dalt = reco_alt - true_alt

plt.figure(figsize=(5, 5))

plt.xlabel("Mis-recontruction [deg]")
plt.ylabel("Number of events")

plt.hist(daz_corr, bins=100, alpha=0.5, label = "azimuth")
plt.hist(dalt, bins=100, alpha=0.5, label = "altitude")

plt.legend()
plt.yscale("log")
plt.grid()

print("Mean and STDs of sky-projected mis-reconstruction axes")
print('daz = {:.4f} +/- {:.4f} deg'.format(daz_corr.mean(), daz_corr.std()))
print('dalt = {:.4f} +/- {:.4f} deg'.format(dalt.mean(), dalt.std()))

plt.show()<jupyter_output><empty_output><jupyter_text>2D representation with **orange** events being those with **offset  20 TeV**<jupyter_code>angcut = (selected_gammas['offset'] < 1) & (selected_gammas['true_energy'] > 20)

plt.figure(figsize=(5,5))
ax = plt.gca()
FOV_size = 2.5 # deg

ax.scatter(daz_corr, dalt, alpha=0.1, s=1, label='no angular cut')
ax.scatter(daz_corr[angcut], dalt[angcut], alpha=0.05, s=1, label='offset < 1 deg & E_true > 20 TeV')

ax.set_aspect('equal')
ax.set_xlabel('cent. Az [deg]')
ax.set_ylabel('cent. Alt [deg]')
ax.set_xlim(-FOV_size,FOV_size)
ax.set_ylim(-FOV_size,FOV_size)
plt.tight_layout()
plt.grid(which="both")

fig.savefig(plots_folder / f"PSFasymmetry_2D_altaz_{analysis_name}.png")<jupyter_output><empty_output><jupyter_text>True energy distributions[back to top](Table-of-contents)<jupyter_code>#min_true_energy = 0.003
#max_true_energy = 330
true_energy_bins_edges = np.logspace(np.log10(min_true_energy), np.log10(max_true_energy), 6 + 1)

if len(np.unique(gammas["true_az"]))==1:
    true_az = np.unique(gammas["true_az"]) * u.deg
    true_alt = np.unique(gammas["true_alt"]) * u.deg
else:
    print("WARNING: diffuse simulations not yet supported.")

print(f"true AZ = {true_az}")
print(f"true ALT = {true_alt}")
get_fig_size(ratio=(9/16), scale=None)
plt.subplots_adjust(wspace=0.5, hspace=0.3)

center = SkyCoord(az=true_az, alt=true_alt, frame="altaz")
nominal_frame = NominalFrame(origin=center)

for i in range(len(true_energy_bins_edges)-1):
    
    plt.subplot(3,2,i+1)
    
    ax = plt.gca()
    ax.set_aspect("equal")

    reconstruction_mask = "is_valid == True and "
    true_energy_mask = f"true_energy > {true_energy_bins_edges[i]} and true_energy < {true_energy_bins_edges[i+1]}"
    selected_gammas = gammas.query(reconstruction_mask + true_energy_mask)

    reconstructed_coordinates = SkyCoord(az=selected_gammas.reco_az.values * u.degree,
                                         alt=selected_gammas.reco_alt.values * u.degree,
                                         frame="altaz")

    reconstructed_coordinates_nominal_frame = reconstructed_coordinates.transform_to(nominal_frame)

    hist_opt = {"bins":[100,100],
                "range":[[-10, 10], [-10, 10]],
                "norm":LogNorm(),
                "cmap":cmap}
    plt.hist2d(reconstructed_coordinates_nominal_frame.fov_lon.value,
               reconstructed_coordinates_nominal_frame.fov_lat.value,
               **hist_opt)
    plt.plot(0, 0, "*", markersize=15, color='#D55E00')
    plt.colorbar()
    plt.xlabel("FOV Longitude [deg]")
    plt.ylabel("FOV Latitude [deg]")
    plt.title(f"{true_energy_bins_edges[i]:.2f} TeV < True energy < {true_energy_bins_edges[i+1]:.2f} TeV")
None<jupyter_output><empty_output><jupyter_text>Same, but with a fixed gamma efficiency cut of 90%<jupyter_code>plt.figure(figsize=(15,15))
plt.subplots_adjust(wspace=0.3, hspace=0.3)

center = SkyCoord(az=true_az, alt=true_alt, frame="altaz")
nominal_frame = NominalFrame(origin=center)

for i in range(len(true_energy_bins_edges)-1):
    
    plt.subplot(3,2,i+1)
    
    ax = plt.gca()
    ax.set_aspect("equal")

    reconstruction_mask = "is_valid == True and "
    true_energy_mask = f"true_energy > {true_energy_bins_edges[i]} and true_energy < {true_energy_bins_edges[i+1]}"
    reconstructed_gammas_per_true_energy = gammas.query(reconstruction_mask + true_energy_mask)
    gammaness = reconstructed_gammas_per_true_energy["gammaness"]
    gammaness_cut = np.quantile(gammaness, gamma_efficiency)
    selected_gammas = reconstructed_gammas_per_true_energy.query(f"gammaness > {gammaness_cut}")
    
    selected_gammas

    reconstructed_coordinates = SkyCoord(az=selected_gammas.reco_az.values * u.degree,
                                         alt=selected_gammas.reco_alt.values * u.degree,
                                         frame="altaz")

    reconstructed_coordinates_nominal_frame = reconstructed_coordinates.transform_to(nominal_frame)

    hist_opt = {"bins":[100,100],
                "range":[[-10, 10], [-10, 10]],
                "norm":LogNorm(),
                "cmap":cmap}
    plt.hist2d(reconstructed_coordinates_nominal_frame.fov_lon.value,
               reconstructed_coordinates_nominal_frame.fov_lat.value,
               **hist_opt)
    plt.plot(0, 0, "*", markersize=20, color='#D55E00')
    plt.colorbar()
    plt.xlabel("FOV Longitude [deg]")
    plt.ylabel("FOV Latitude [deg]")
    plt.title(f"{true_energy_bins_edges[i]:.2f} TeV < True energy < {true_energy_bins_edges[i+1]:.2f} TeV")
None<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12943883688091626, 'fraction_numerical': 0.021494691700224885, 'mean_word_length': 4.536189924724956, 'pattern_counts': {'":': 13, '<': 63, '<?xml version=': 0, '>': 62, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '667614', 'license': 'CECILL-B', 'n_tokens_mistral': 7169, 'n_tokens_neox': 6693, 'n_words': 1235, 'path': 'protopipe/benchmarks/notebooks/DL2/benchmarks_DL2_direction-reconstruction.ipynb', 'repo_name': 'Iburelli/protopipe', 'size': '19121'}