DL2 performances¶

This is an example of how to compute reconstruction performances on this DL2 dataset.
We explicitly note that the products provided are preliminary and do not reflect the final performance of the CTA Observatory.

Note that this dataset does not include energy reconstruction nor gamma/hadron separation. Therefore, only the angular resolution and effective area can be computed.

In [1]:

import ctapipe
print(ctapipe.__version__)

import matplotlib.pyplot as plt
import numpy as np

import ctapipe print(ctapipe.__version__) import matplotlib.pyplot as plt import numpy as np

0.17.0

In [2]:

filename = '../data/gamma-diffuse.dl2.h5'

filename = '../data/gamma-diffuse.dl2.h5'

In [3]:

from ctapipe.io import TableLoader

def read_table(filename, tel_ids=None, stop=None):
    """
    Read data from a file and return it as a pandas DataFrame.

    Parameters
    ----------
    filename : str
        The path to the input file.
    tel_ids : list or None, optional
        A list of telescope IDs to include in the data. If None, all telescopes are included. Default is None.
    stop : int or None, optional
        The number of events to read from the file. If None, all events are read. Default is None.
    Returns
    -------
    pandas.DataFrame
        The data read from the file as a pandas DataFrame.
    """
    loader = TableLoader(
        input_url=filename,
        load_dl1_parameters=True,
        load_dl2=True,
        load_instrument=True,
        load_simulated=True,
        load_true_parameters=True,
    )

    data = loader.read_telescope_events(telescopes=tel_ids, stop=stop)
    data.remove_columns([col for col in data.colnames if len(data[col].shape) > 1])
    return data

from ctapipe.io import TableLoader def read_table(filename, tel_ids=None, stop=None): """ Read data from a file and return it as a pandas DataFrame. Parameters ---------- filename : str The path to the input file. tel_ids : list or None, optional A list of telescope IDs to include in the data. If None, all telescopes are included. Default is None. stop : int or None, optional The number of events to read from the file. If None, all events are read. Default is None. Returns ------- pandas.DataFrame The data read from the file as a pandas DataFrame. """ loader = TableLoader( input_url=filename, load_dl1_parameters=True, load_dl2=True, load_instrument=True, load_simulated=True, load_true_parameters=True, ) data = loader.read_telescope_events(telescopes=tel_ids, stop=stop) data.remove_columns([col for col in data.colnames if len(data[col].shape) > 1]) return data

/home/runner/micromamba/envs/ctaodl2/lib/python3.9/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

In [4]:

# limited stats because of github actions limitations => THE PERFORMANCES ARE NOT REPRESENTATIVE
dl2_params = read_table(filename, stop=100000)

# limited stats because of github actions limitations => THE PERFORMANCES ARE NOT REPRESENTATIVE dl2_params = read_table(filename, stop=100000)

In [5]:

import astropy.units as u
from astropy.coordinates import angular_separation
from pyirf.benchmarks import angular_resolution
from pyirf.binning import create_bins_per_decade

dl2_params['theta'] = angular_separation(dl2_params['true_alt'], dl2_params['true_az'], dl2_params['HillasReconstructor_alt'], dl2_params['HillasReconstructor_az'])

ang_res = angular_resolution(dl2_params, 
    energy_bins=np.array([dl2_params['true_energy'].min(), dl2_params['true_energy'].max()])*dl2_params['true_energy'].unit)['angular_resolution_68'][0]

plt.hist(dl2_params['theta'].to(u.deg)**2, bins=100, range=(0*u.deg**2, 1*u.deg**2))
plt.axvline(ang_res.to_value(u.deg)**2, color='r', label='angular resolution (68%) = {:.2f} deg'.format(ang_res.to_value(u.deg)))
plt.xlabel(r'$\theta^2$ (deg$^2$)')
plt.ylabel('Number of events')
plt.legend()

import astropy.units as u from astropy.coordinates import angular_separation from pyirf.benchmarks import angular_resolution from pyirf.binning import create_bins_per_decade dl2_params['theta'] = angular_separation(dl2_params['true_alt'], dl2_params['true_az'], dl2_params['HillasReconstructor_alt'], dl2_params['HillasReconstructor_az']) ang_res = angular_resolution(dl2_params, energy_bins=np.array([dl2_params['true_energy'].min(), dl2_params['true_energy'].max()])*dl2_params['true_energy'].unit)['angular_resolution_68'][0] plt.hist(dl2_params['theta'].to(u.deg)**2, bins=100, range=(0*u.deg**2, 1*u.deg**2)) plt.axvline(ang_res.to_value(u.deg)**2, color='r', label='angular resolution (68%) = {:.2f} deg'.format(ang_res.to_value(u.deg))) plt.xlabel(r'$\theta^2$ (deg$^2$)') plt.ylabel('Number of events') plt.legend()

Out[5]:

<matplotlib.legend.Legend at 0x7f5a6ed63880>

In [6]:

energy_bins = create_bins_per_decade(1e-2*u.TeV, 1e2*u.TeV, 10)
ang_res = angular_resolution(dl2_params, energy_bins=energy_bins)


plt.figure(figsize=(10, 6))
plt.step(ang_res['true_energy_center'], ang_res['angular_resolution_68'], where='mid')
plt.xscale('log')
plt.yscale('log')
plt.xlabel('True Energy (TeV)')
plt.ylabel('Angular Resolution (deg)')
plt.title('Angular Resolution vs True Energy')
plt.grid(True, which="both", ls="--")
plt.show()

energy_bins = create_bins_per_decade(1e-2*u.TeV, 1e2*u.TeV, 10) ang_res = angular_resolution(dl2_params, energy_bins=energy_bins) plt.figure(figsize=(10, 6)) plt.step(ang_res['true_energy_center'], ang_res['angular_resolution_68'], where='mid') plt.xscale('log') plt.yscale('log') plt.xlabel('True Energy (TeV)') plt.ylabel('Angular Resolution (deg)') plt.title('Angular Resolution vs True Energy') plt.grid(True, which="both", ls="--") plt.show()

In [7]:

from pyirf.simulations import SimulatedEventsInfo

loader = TableLoader(filename)
simulation_table = loader.read_simulation_configuration()
simulation_table[:10]


simulation_info = SimulatedEventsInfo(
    energy_min=simulation_table['energy_range_min'][0] * simulation_table['energy_range_min'].unit,
    energy_max=simulation_table['energy_range_max'][0] * simulation_table['energy_range_max'].unit,
    spectral_index=simulation_table['spectral_index'][0],
    max_impact=simulation_table['max_scatter_range'][0] * simulation_table['max_scatter_range'].unit,
    n_showers=np.sum(simulation_table['n_showers'])*simulation_table['shower_reuse'][0],
    viewcone_min=simulation_table['min_viewcone_radius'][0] * simulation_table['min_viewcone_radius'].unit,
    viewcone_max=simulation_table['max_viewcone_radius'][0] * simulation_table['max_viewcone_radius'].unit,
)

simulation_info

from pyirf.simulations import SimulatedEventsInfo loader = TableLoader(filename) simulation_table = loader.read_simulation_configuration() simulation_table[:10] simulation_info = SimulatedEventsInfo( energy_min=simulation_table['energy_range_min'][0] * simulation_table['energy_range_min'].unit, energy_max=simulation_table['energy_range_max'][0] * simulation_table['energy_range_max'].unit, spectral_index=simulation_table['spectral_index'][0], max_impact=simulation_table['max_scatter_range'][0] * simulation_table['max_scatter_range'].unit, n_showers=np.sum(simulation_table['n_showers'])*simulation_table['shower_reuse'][0], viewcone_min=simulation_table['min_viewcone_radius'][0] * simulation_table['min_viewcone_radius'].unit, viewcone_max=simulation_table['max_viewcone_radius'][0] * simulation_table['max_viewcone_radius'].unit, ) simulation_info

Out[7]:

SimulatedEventsInfo(n_showers=2150000000, energy_min=0.003 TeV, energy_max=330.00 TeV, spectral_index=-2.0, max_impact=1900.00 m, viewcone_min=0.0 degviewcone_max=10.0 deg)

In [8]:

from pyirf.irf import effective_area_per_energy


# Compute the effective area
effective_area = effective_area_per_energy(
    dl2_params,
    simulation_info,
    energy_bins
)

# Plot the effective area
plt.figure(figsize=(10, 6))
plt.step(energy_bins[:-1], effective_area, where='mid')
plt.xscale('log')
plt.yscale('log')
plt.xlabel('Energy (TeV)')
plt.ylabel('Effective Area (m²)')
plt.title('Effective Area vs Energy')
plt.grid(True, which="both", ls="--")
plt.show()

from pyirf.irf import effective_area_per_energy # Compute the effective area effective_area = effective_area_per_energy( dl2_params, simulation_info, energy_bins ) # Plot the effective area plt.figure(figsize=(10, 6)) plt.step(energy_bins[:-1], effective_area, where='mid') plt.xscale('log') plt.yscale('log') plt.xlabel('Energy (TeV)') plt.ylabel('Effective Area (m²)') plt.title('Effective Area vs Energy') plt.grid(True, which="both", ls="--") plt.show()

In [ ]: