In [1]:
import sys
sys.path.append('/home/jbourbeau/cr-composition')
sys.path
Out[1]:
In [2]:
from __future__ import division
import numpy as np
import sys
import matplotlib.pyplot as plt
import seaborn.apionly as sns
# from composition.analysis import load_sim
# import composition.analysis.plotting_functions as plotting
# import composition.analysis.data_functions as data_functions
import composition as comp
%matplotlib inline
In [3]:
sns.set_palette('muted')
sns.set_color_codes()
In [4]:
fig, ax = plt.subplots()
comp_list = ['P', 'He', 'Fe']
for composition in comp_list:
standard_cut_keys = ['lap_reco_success', 'lap_zenith', 'num_hits_1_30', 'IceTopMaxSignalInEdge',
'IceTopMaxSignal', 'IceTopNeighbourMaxSignal',
'StationDensity', 'max_qfrac_1_30', 'lap_containment', 'energy_range_lap']
labels = ['Reco success', 'Reco zenith', 'NChannel/NStation', 'IceTopMaxSignalInEdge',
'IceTopMaxSignal', 'IceTopNeighbourMaxSignal',
'StationDensity', 'InIce charge frac', 'Reco containment', '$6.2 < \log_{10}(\mathrm{E/GeV}) < 8.0$']
# Import ShowerLLH sim reconstructions and cuts to be made
df, cut_dict = comp.load_sim(return_cut_dict=True)
selection_mask = (df.MC_comp == composition).values
# selection_mask = np.array([True]*len(df))
n_tot = len(df)
n_tot_err = np.sqrt(len(df))
initial, initial_err = comp.ratio_error(
n_tot, n_tot_err,
n_tot, n_tot_err)
event_flow = [initial]
event_flow_err = [initial_err]
for cut_key, label in zip(standard_cut_keys, labels):
selection_mask = selection_mask & cut_dict[cut_key]
df_cuts = df[selection_mask]
ratio, ratio_err = comp.ratio_error(
len(df_cuts), np.sqrt(len(df_cuts)),
n_tot, n_tot_err)
event_flow.append(ratio)
event_flow_err.append(ratio_err)
# event_flow.append(len(df_cuts)/n_tot)
# event_flow_err.append(np.sqrt(len(df_cuts))/n_tot)
# fig, ax = plt.subplots()
x = range(len(standard_cut_keys)+1)
print(x)
print(event_flow)
ax.errorbar(x, event_flow, yerr=event_flow_err,
marker='.', linestyle='-.', label=composition, alpha=0.75)
ax.grid()
# ax.set_ylim([1e-1,1e5])
# ax.set_yscale("log", nonposy='clip')
# ax.set_xlabel('True IT Containment Fraction')
plt.xticks(x, ['Coincident'] + labels, rotation='vertical')
ax.set_ylabel('Fraction of events')
ax.set_xlim([-1, len(x)])
ax.set_ylim([0, 1.1])
plt.legend(title='MC composition')
plt.show()
In [ ]: