# Remove input cells at runtime (nbsphinx)
import IPython.core.display as d
d.display_html('<script>jQuery(function() {if (jQuery("body.notebook_app").length == 0) { jQuery(".input_area").toggle(); jQuery(".prompt").toggle();}});</script>', raw=True)

Performance poster layout

This notebook produces a nice poster layout about the comparison between pipelines.

It is just a refurbished version of the DL3-level “IRF and sensitivity” notebook.

Latest performance results cannot be shown on this public documentation and are therefore hosted at this RedMine page .

## From the standard library
import os
from pathlib import Path

# From pyirf
import pyirf
from pyirf.binning import bin_center
from pyirf.utils import cone_solid_angle

# From other 3rd-party libraries
from yaml import load, FullLoader
import numpy as np
import astropy.units as u
from astropy.io import fits
from astropy.table import QTable, Table, Column
import uproot

import matplotlib.pyplot as plt
from matplotlib.ticker import ScalarFormatter
%matplotlib inline
plt.rcParams['axes.labelsize'] = 15
plt.rcParams['xtick.labelsize'] = 15
plt.rcParams['ytick.labelsize'] = 15

def load_config(name):
    """Load YAML configuration file."""
    try:
        with open(name, "r") as stream:
            cfg = load(stream, Loader=FullLoader)
    except FileNotFoundError as e:
        print(e)
        raise
    return cfg

# EDIT THIS CELL WITH YOUR LOCAL SETUP INFORMATION
parent_dir = "" # path to 'analyses' folder
analysisName = ""
infile = "performance_protopipe_Prod3b_CTANorth_baseline_full_array_Zd20deg_180deg_Time50.00h.fits.gz"

config_performance = load_config(f"{parent_dir}/{analysisName}/configs/performance.yaml")
obs_time = f'{config_performance["analysis"]["obs_time"]["value"]}{config_performance["analysis"]["obs_time"]["unit"]}'
production = infile.split("protopipe_")[1].split("_Time")[0]
protopipe_file = Path(parent_dir, analysisName, "data/DL3", infile)

# EDIT THIS CELL WITH YOUR LOCAL SETUP INFORMATION

parent_dir_aswg = ""

# MARS performance (available here: https://forge.in2p3.fr/projects/step-by-step-reference-mars-analysis/wiki)
indir_CTAMARS = ""
infile_CTAMARS = "SubarrayLaPalma_4L15M_south_IFAE_50hours_20190630.root"
MARS_label = "CTAMARS (2019)"

# ED performance (available here: https://forge.in2p3.fr/projects/cta_analysis-and-simulations/wiki/Prod3b_based_instrument_response_functions)
indir_ED = ""
infile_ED = "CTA-Performance-North-20deg-S-50h_20181203.root"
ED_label = "EventDisplay (2018)"

MARS_performance = uproot.open(Path(parent_dir_aswg, indir_CTAMARS, infile_CTAMARS))
ED_performance = uproot.open(Path(parent_dir_aswg, indir_ED, infile_ED))

# EDIT THIS CELL WITH YOUR LOCAL SETUP INFORMATION
indir_requirements = ''
site = 'North'
obs_time = '50h'

infile_requirements = {"sens" : f'/{site}-{obs_time}.dat',
                       "AngRes" : f'/{site}-{obs_time}-AngRes.dat',
                       "ERes" : f'/{site}-{obs_time}-ERes.dat'}
requirements = {}

for key in infile_requirements.keys():
    requirements[key] = Table.read(indir_requirements + infile_requirements[key], format='ascii')
requirements['sens'].add_column(Column(data=(10**requirements['sens']['col1']), name='ENERGY'))
requirements['sens'].add_column(Column(data=requirements['sens']['col2'], name='SENSITIVITY'))

# First we check if a _plots_ folder exists already.
# If not, we create it.
Path("./plots").mkdir(parents=True, exist_ok=True)

fig = plt.figure(figsize = (20, 10), constrained_layout=True)

gs = fig.add_gridspec(3, 3, figure=fig)

# ==========================================================================================================
#
#                                       SENSITIVITY
#
# ==========================================================================================================

ax1 = fig.add_subplot(gs[0:-1, 0:-1])

# [1:-1] removes under/overflow bins
sensitivity_protopipe = QTable.read(protopipe_file, hdu='SENSITIVITY')[1:-1]

unit = u.Unit('erg cm-2 s-1')

# Add requirements
ax1.plot(requirements['sens']['ENERGY'],
         requirements['sens']['SENSITIVITY'],
         color='black',
         ls='--',
         lw=2,
         label='Requirements'
)

# protopipe
e = sensitivity_protopipe['reco_energy_center']
w = (sensitivity_protopipe['reco_energy_high'] - sensitivity_protopipe['reco_energy_low'])
s_p = (e**2 * sensitivity_protopipe['flux_sensitivity'])
ax1.errorbar(
    e.to_value(u.TeV),
    s_p.to_value(unit),
    xerr=w.to_value(u.TeV) / 2,
    ls='',
    label='protopipe',
    color='DarkOrange'
)

# ED
s_ED, edges = ED_performance["DiffSens"].to_numpy()
yerr = ED_performance["DiffSens"].errors()
bins = 10**edges
x = bin_center(bins)
width = np.diff(bins)
ax1.errorbar(
    x,
    s_ED,
    xerr=width/2,
    yerr=yerr,
    label=ED_label,
    ls='',
    color='DarkGreen'
)

# MARS
s_MARS, edges = MARS_performance["DiffSens"].to_numpy()
yerr = MARS_performance["DiffSens"].errors()
bins = 10**edges
x = bin_center(bins)
width = np.diff(bins)
ax1.errorbar(
    x,
    s_MARS,
    xerr=width/2,
    yerr=yerr,
    label=MARS_label,
    ls='',
    color='DarkBlue'
)

# Style settings
ax1.set_xscale("log")
ax1.set_yscale("log")
ax1.set_ylabel(rf"$(E^2 \cdot \mathrm{{Flux Sensitivity}}) /$ ({unit.to_string('latex')})")
ax1.grid(which="both")
ax1.legend()

# ==========================================================================================================
#
#                                       SENSITIVITY RATIO
#
# ==========================================================================================================



ax2 = fig.add_subplot(gs[2, 0])

ax2.errorbar(
    e.to_value(u.TeV),
    s_p.to_value(unit) / s_ED,
    xerr=w.to_value(u.TeV)/2,
    ls='',
    label = "",
    color='DarkGreen'
)
ax2.errorbar(
    e.to_value(u.TeV),
    s_p.to_value(unit) / s_MARS,
    xerr=w.to_value(u.TeV)/2,
    ls='',
    label = "",
    color='DarkBlue'
)
ax2.axhline(1, color = 'DarkOrange')

ax2.set_xscale('log')
ax2.set_yscale('log')
ax2.set_xlabel("Reconstructed energy [TeV]")
ax2.set_ylabel('Sensitivity ratio')
ax2.grid()
ax2.yaxis.set_major_formatter(ScalarFormatter())

ax2.set_ylim(0.5, 2.0)
ax2.set_yticks([0.5, 2/3, 1, 3/2, 2])
ax2.set_yticks([], minor=True)

# ==========================================================================================================
#
#                                       EFFECTIVE COLLECTION AREA
#
# ==========================================================================================================

ax3 = fig.add_subplot(gs[0, 2])

# protopipe
# uncomment the other strings to see effective areas
# for the different cut levels. Left out here for better
# visibility of the final effective areas.
suffix =''
#'_NO_CUTS'
#'_ONLY_GH'
#'_ONLY_THETA'

area = QTable.read(protopipe_file, hdu='EFFECTIVE_AREA' + suffix)[0]
ax3.errorbar(
    0.5 * (area['ENERG_LO'] + area['ENERG_HI']).to_value(u.TeV)[1:-1],
    area['EFFAREA'].to_value(u.m**2).T[1:-1, 0],
    xerr=0.5 * (area['ENERG_LO'] - area['ENERG_HI']).to_value(u.TeV)[1:-1],
    ls='',
    label='protopipe ' + suffix,
    color='DarkOrange'
)

# ED
y, edges = ED_performance["EffectiveAreaEtrue"].to_numpy()
yerr = ED_performance["EffectiveAreaEtrue"].errors()
x = bin_center(10**edges)
xerr = 0.5 * np.diff(10**edges)
ax3.errorbar(x,
             y,
             xerr=xerr,
             yerr=yerr,
             ls='',
             label=ED_label,
             color='DarkGreen'
            )

# MARS
y, edges = MARS_performance["EffectiveAreaEtrue"].to_numpy()
yerr = MARS_performance["EffectiveAreaEtrue"].errors()
x = bin_center(10**edges)
xerr = 0.5 * np.diff(10**edges)
ax3.errorbar(x,
             y,
             xerr=xerr,
             yerr=yerr,
             ls='',
             label=MARS_label,
             color='DarkBlue'
            )

# Style settings
ax3.set_xscale("log")
ax3.set_yscale("log")
ax3.set_xlabel("True energy [TeV]")
ax3.set_ylabel("Effective area [m²]")
ax3.grid(which="both")

# ==========================================================================================================
#
#                                       ANGULAR RESOLUTION
#
# ==========================================================================================================


ax4 = fig.add_subplot(gs[2, 1])

# protopipe
ang_res = QTable.read(protopipe_file, hdu='ANGULAR_RESOLUTION')[1:-1]

ax4.errorbar(
    0.5 * (ang_res['reco_energy_low'] + ang_res['reco_energy_high']).to_value(u.TeV),
    ang_res['angular_resolution'].to_value(u.deg),
    xerr=0.5 * (ang_res['reco_energy_high'] - ang_res['reco_energy_low']).to_value(u.TeV),
    ls='',
    label='protopipe',
    color='DarkOrange'
)

# ED
y, edges = ED_performance["AngRes"].to_numpy()
yerr = ED_performance["AngRes"].errors()
x = bin_center(10**edges)
xerr = 0.5 * np.diff(10**edges)
ax4.errorbar(x,
             y,
             xerr=xerr,
             yerr=yerr,
             ls='',
             label=ED_label,
             color='DarkGreen')

# MARS
y, edges = MARS_performance["AngRes"].to_numpy()
yerr = MARS_performance["AngRes"].errors()
x = bin_center(10**edges)
xerr = 0.5 * np.diff(10**edges)
ax4.errorbar(x,
             y,
             xerr=xerr,
             yerr=yerr,
             ls='',
             label=MARS_label,
             color='DarkBlue')

# Requirements
ax4.plot(10**requirements['AngRes']['col1'],
         requirements['AngRes']['col2'],
         color='black',
         ls='--',
         lw=2,
         label='Requirements'
)

# Style settings
ax4.set_xscale("log")
ax4.set_yscale("log")
ax4.set_xlabel("Reconstructed energy [TeV]")
ax4.set_ylabel("Angular resolution [deg]")
ax4.grid(which="both")

None # to remove clutter by mpl objects

# ==========================================================================================================
#
#                                       ENERGY RESOLUTION
#
# ==========================================================================================================


ax5 = fig.add_subplot(gs[2, 2])

# protopipe
bias_resolution = QTable.read(protopipe_file, hdu='ENERGY_BIAS_RESOLUTION')[1:-1]
ax5.errorbar(
    0.5 * (bias_resolution['reco_energy_low'] + bias_resolution['reco_energy_high']).to_value(u.TeV),
    bias_resolution['resolution'],
    xerr=0.5 * (bias_resolution['reco_energy_high'] - bias_resolution['reco_energy_low']).to_value(u.TeV),
    ls='',
    label='protopipe',
    color='DarkOrange'
)


# ED
y, edges = ED_performance["ERes"].to_numpy()
yerr = ED_performance["ERes"].errors()
x = bin_center(10**edges)
xerr = np.diff(10**edges) / 2
ax5.errorbar(x,
             y,
             xerr=xerr,
             yerr=yerr,
             ls='',
             label=ED_label,
             color='DarkGreen'
            )

# MARS
y, edges = MARS_performance["ERes"].to_numpy()
yerr = MARS_performance["ERes"].errors()
x = bin_center(10**edges)
xerr = np.diff(10**edges) / 2
ax5.errorbar(x,
             y,
             xerr=xerr,
             yerr=yerr,
             ls='',
             label=MARS_label,
             color='DarkBlue'
            )

# Requirements
ax5.plot(10**requirements['ERes']['col1'],
         requirements['ERes']['col2'],
         color='black',
         ls='--',
         lw=2,
         label='Requirements'
)

# Style settings
ax5.set_xlabel("Reconstructed energy [TeV]")
ax5.set_ylabel("Energy resolution")
ax5.grid(which="both")
ax5.set_xscale('log')

None # to remove clutter by mpl objects

# ==========================================================================================================
#
#                                       BACKGROUND RATE
#
# ==========================================================================================================



ax6 = fig.add_subplot(gs[1, 2])

from pyirf.utils import cone_solid_angle

# protopipe
rad_max = QTable.read(protopipe_file, hdu='RAD_MAX')[0]
bg_rate = QTable.read(protopipe_file, hdu='BACKGROUND')[0]

reco_bins = np.append(bg_rate['ENERG_LO'], bg_rate['ENERG_HI'][-1])

# first fov bin, [0, 1] deg
fov_bin = 0
rate_bin = bg_rate['BKG'].T[:, fov_bin]

# interpolate theta cut for given e reco bin
e_center_bg = 0.5 * (bg_rate['ENERG_LO'] + bg_rate['ENERG_HI'])
e_center_theta = 0.5 * (rad_max['ENERG_LO'] + rad_max['ENERG_HI'])
theta_cut = np.interp(e_center_bg, e_center_theta, rad_max['RAD_MAX'].T[:, 0])

# undo normalization
rate_bin *= cone_solid_angle(theta_cut)
rate_bin *= np.diff(reco_bins)
ax6.errorbar(
    0.5 * (bg_rate['ENERG_LO'] + bg_rate['ENERG_HI']).to_value(u.TeV)[1:-1],
    rate_bin.to_value(1 / u.s)[1:-1],
    xerr=np.diff(reco_bins).to_value(u.TeV)[1:-1] / 2,
    ls='',
    label='protopipe',
    color='DarkOrange'
)

# ED
y, edges = ED_performance["BGRate"].to_numpy()
yerr = ED_performance["BGRate"].errors()
x = bin_center(10**edges)
xerr = np.diff(10**edges) / 2
ax6.errorbar(x,
             y,
             xerr=xerr,
             yerr=yerr,
             ls='',
             label=ED_label,
             color="DarkGreen")


# MARS
y, edges = MARS_performance["BGRate"].to_numpy()
yerr = MARS_performance["BGRate"].errors()
x = bin_center(10**edges)
xerr = np.diff(10**edges) / 2
ax6.errorbar(x,
             y,
             xerr=xerr,
             yerr=yerr,
             ls='',
             label=MARS_label,
             color="DarkBlue")


# Style settings
ax6.set_xscale("log")
ax6.set_xlabel("Reconstructed energy [TeV]")
ax6.set_ylabel("Background rate [s⁻¹ TeV⁻¹] ")
ax6.grid(which="both")
ax6.set_yscale('log')

fig.suptitle(f'{production} - {obs_time}', fontsize=25)
fig.savefig(f"./plots/protopipe_{production}_{obs_time}.png")
None # to remove clutter by mpl objects