In [1]:

    

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

%matplotlib inline


    

    




/home/jbourbeau/.local/lib/python2.7/site-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.
  warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')

In [2]:

    

sns.set_palette('Paired', 10)
# sns.set_palette(sns.color_palette('Blues', 10))
sns.set_color_codes()


    

In [3]:

    

standard_cut_keys = ['reco_exists', 'reco_zenith', 'num_hits', 'IT_signal',
                     'StationDensity', 'max_charge_frac', 'reco_containment', 'energy_range']
labels = ['Reco success', 'Reco zenith', 'NChannel/NStation', 'IT signal cuts',
          'StationDensity', 'InIce charge frac', 'Reco containment', 'Energy range']
# Import ShowerLLH sim reconstructions and cuts to be made
df, cut_dict = load_sim(return_cut_dict=True)
selection_mask = np.array([True]*len(df))
fig, ax = plt.subplots()
containment_bins = np.linspace(0, 6, 100)
IT_containment = df[selection_mask]['IceTop_FractionContainment']
print('IT fraction contained = {}'.format(sum(IT_containment <= 1.0)/len(IT_containment)))
sns.distplot(IT_containment, containment_bins, kde=False,
             hist_kws={"histtype": "step", "linewidth": 1,
                       "alpha": 1}, 
             label='Coincident')
for cut_key, label in zip(standard_cut_keys, labels):
    selection_mask = selection_mask & cut_dict[cut_key]
    df_cuts = df[selection_mask]
    
    IT_containment = df_cuts['IceTop_FractionContainment']
    print('IT fraction contained = {}'.format(sum(IT_containment <= 1.0)/len(IT_containment)))
    sns.distplot(IT_containment, containment_bins, kde=False,
                 hist_kws={"histtype": "step", "linewidth": 1,
                            "alpha": 1}, 
                label='+ ' + label)
# axarr[0].legend()
ax.legend()
ax.set_ylim([1e-1,1e5])
ax.set_yscale("log", nonposy='clip')
ax.set_xlabel('True IT Containment Fraction')
ax.set_ylabel('Counts')


    

    




/home/jbourbeau/composition/analysis/load_sim.py:41: RuntimeWarning: divide by zero encountered in log10
  cut_dict['min_energy'] = (np.log10(df['reco_energy']) > 6.2)
/home/jbourbeau/composition/analysis/load_sim.py:42: RuntimeWarning: divide by zero encountered in log10
  cut_dict['max_energy'] = (np.log10(df['reco_energy']) < 8.0)
/home/jbourbeau/composition/analysis/load_sim.py:63: RuntimeWarning: divide by zero encountered in log10
  df['reco_log_energy'] = np.nan_to_num(np.log10(df['reco_energy']))
/home/jbourbeau/composition/analysis/load_sim.py:64: RuntimeWarning: invalid value encountered in log10
  df['InIce_log_charge'] = np.nan_to_num(np.log10(df['InIce_charge']))






    




IT fraction contained = 0.387007866437
IT fraction contained = 0.456563103363
IT fraction contained = 0.493782516831
IT fraction contained = 0.506261571568
IT fraction contained = 0.997818562427
IT fraction contained = 0.997818471998
IT fraction contained = 0.997817872711
IT fraction contained = 0.999966636194
IT fraction contained = 0.999980728093






    
Out[3]:





<matplotlib.text.Text at 0x57b4510>






    

In [4]:

    

selection_mask = np.array([True]*len(df))
InIce_containment = df[selection_mask]['LLHlap_InIce_containment']
# InIce_containment = df[selection_mask]['InIce_FractionContainment']
print('InIce fraction contained = {}'.format(sum(InIce_containment <= 1.0)/len(InIce_containment)))
fig, ax = plt.subplots()
sns.distplot(InIce_containment, containment_bins, kde=False,
             hist_kws={"histtype": "step", "linewidth": 1,
                       "alpha": 1}, 
             label='Coincident')
for cut_key, label in zip(standard_cut_keys, labels):
    selection_mask = selection_mask & cut_dict[cut_key]
    df_cuts = df[selection_mask]
    
    InIce_containment = df_cuts['InIce_FractionContainment']
    print('InIce fraction contained = {}'.format(sum(InIce_containment <= 1.0)/len(InIce_containment)))
    sns.distplot(InIce_containment, containment_bins, kde=False,
                 hist_kws={"histtype": "step", "linewidth": 1,
                            "alpha": 1}, 
                label='+ ' + label)
    
ax.legend()
ax.set_ylim([1e-1,1e5])
ax.set_yscale("log", nonposy='clip')
ax.set_xlabel('True InIce Containment Fraction')
ax.set_ylabel('Counts')
plt.tight_layout()


    

    




InIce fraction contained = 0.792888614105
InIce fraction contained = 0.716529809532
InIce fraction contained = 0.72037797148
InIce fraction contained = 0.738627429176
InIce fraction contained = 0.775778270602
InIce fraction contained = 0.775774157443
InIce fraction contained = 0.775873239656
InIce fraction contained = 0.974022940953
InIce fraction contained = 0.974869432828






    

In [11]:

    

# Import ShowerLLH sim reconstructions and cuts to be made
df, cut_dict = load_sim(return_cut_dict=True)
selection_mask = np.array([True] * len(df))
standard_cut_keys = ['reco_exists', 'reco_zenith', 'num_hits', 'IceTopMaxSignalInEdge', 'IceTopMaxSignal', 'IceTopNeighbourMaxSignal',
                     'StationDensity', 'max_charge_frac', 'reco_containment', 'min_energy']
for key in standard_cut_keys:
    selection_mask *= cut_dict[key]

df = df[selection_mask]
print('number of events = {}'.format(len(df)))

InIce_containment = df['InIce_FractionContainment']
log_energy = df['reco_log_energy']
MC_log_energy = df['MC_log_energy']
charge = df.InIce_log_charge

containment_bins = np.linspace(0, 1, 75)
energy_bins = np.arange(6.2, 9.51, 0.05)
bins = [containment_bins, energy_bins]

fig, ax = plt.subplots()
im = plotting.histogram_2D(InIce_containment, MC_log_energy, bins, log_counts=True, ax=ax)
# ax.set_title(title)
# plt.xlim([6.2, 9.5])
plt.show()


    

    




number of events = 116310






    




/home/jbourbeau/composition/analysis/plotting_functions.py:59: RuntimeWarning: divide by zero encountered in log10
  h = np.log10(h)






    

In [6]:

    

sns.palplot(sns.color_palette("Paired", 10))


    

In [ ]: