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.load_sim import load_sim
from composition.support_functions.checkdir import checkdir
import composition.analysis.plotting_functions as plotting

%matplotlib inline


    

In [2]:

    

def make_charge_energy_histogram(log_energy, charge, proton_mask, iron_mask, ax, plot_line=False):
    charge_bins = np.linspace(0, 6, 50)
    energy_bins = np.arange(6.2, 8.0, 0.05)
#     energy_bins = np.arange(6.2, 9.51, 0.05)
    energy_midpoints = (energy_bins[1:] + energy_bins[:-1]) / 2
    proton_hist, xedges, yedges = np.histogram2d(log_energy[MC_proton_mask],
                                                 charge[MC_proton_mask],
                                                 bins=[energy_bins,
                                                       charge_bins],
                                                 normed=False)
    proton_hist = np.ma.masked_where(proton_hist == 0, proton_hist)
    iron_hist, xedges, yedges = np.histogram2d(log_energy[MC_iron_mask],
                                               charge[MC_iron_mask],
                                               bins=[energy_bins, charge_bins],
                                               normed=False)

    h = proton_hist / (proton_hist + iron_hist)
    h = np.rot90(h)
    h = np.flipud(h)
    h = np.ma.masked_where(h == 0, h)
    extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]]
    
    def line_fit(array):
        fit = []
        for x in array:
            if x <= 9.0:
                slope = (5.3 - 2.55) / (9.5 - 6.2)
                fit.append(2.55 + slope * (x - 6.2))
            else:
                slope = (5.20 - 4.9) / (9.5 - 9.0)
                fit.append(4.9 + slope * (x - 9.0))
        fit = np.array(fit)
        return fit

    colormap = 'coolwarm'
    im = ax.imshow(h, extent=extent, origin='lower',
               interpolation='none', cmap=colormap,
               aspect='auto', vmin=0, vmax=1)
    x = np.arange(6.2, 9.51, 0.1)
    if plot_line:
        ax.plot(x, line_fit(x), marker='None', linestyle='--',
                 color='k')
    return im


    

In [7]:

    

charge_list = ['','_threequarter', '_half', '_quarter']
title_list = ['Full detector','Top 75\%', 'Top 50\%', 'Top 25\%']

fig, axarr = plt.subplots(2,2)
for charge, title, ax in zip(charge_list, title_list, axarr.flatten()):
    # 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', 'NStations', 'IT_signal',
                         'StationDensity', 'reco_containment', 'energy_range']
    for key in standard_cut_keys:
        selection_mask *= cut_dict[key]
    
#     selection_mask *= (df['LLHlap_InIce_containment'] < 1.0)
#     selection_mask *= (df['NChannels{}'.format(charge)] >= 8)
#     selection_mask *= (df['max_charge_frac{}'.format(charge)] < 0.3)
    df = df[selection_mask]
    print('number of events = {}'.format(len(df)))

    MC_proton_mask = (df.MC_comp == 'P')
    MC_iron_mask = (df.MC_comp == 'Fe')
    log_energy = df['reco_log_energy']
    charge = df['InIce_log_charge{}'.format(charge)]
    im = make_charge_energy_histogram(log_energy, charge, MC_proton_mask, MC_iron_mask, ax=ax)
    ax.set_title(title)

fig.text(0.5, 0.00, '$\log_{10}(E_{\mathrm{reco}}/\mathrm{GeV})$', ha='center')
fig.text(0.00, 0.5, '$\log_{10}(Q_{\mathrm{total}})$', va='center', rotation='vertical')
plt.tight_layout()
fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
fig.colorbar(im, cax=cbar_ax, label='P/(P+Fe) Fraction')
# plt.colorbar(im, label='P/(P+Fe) Fraction')


    

    




number of events = 51612
number of events = 51612
number of events = 51612
number of events = 51612






    
Out[7]:





<matplotlib.colorbar.Colorbar at 0x7586990>






    

In [2]:

    

from pandas.tools.plotting import parallel_coordinates


    

In [3]:

    

df = load_sim()
parallel_coordinates(df, class_column='MC_comp', cols=['InIce_charge', 'MC_energy'])


    

    




/home/jbourbeau/composition/analysis/load_sim.py:70: RuntimeWarning: divide by zero encountered in log10
  df['reco_log_energy'] = np.nan_to_num(np.log10(df['reco_energy']))






    




Cut event flow:
                   reco_exists:  0.645  0.645
                   reco_zenith:  0.809  0.597
                      num_hits:  0.625  0.579
                     IT_signal:  0.227  0.218
                StationDensity:  0.787  0.218
              reco_containment:    0.5  0.169
               max_charge_frac:  0.775  0.169
                  energy_range:  0.447  0.117








    
Out[3]:





<matplotlib.axes._subplots.AxesSubplot at 0x4b50510>






    





RuntimeErrorTraceback (most recent call last)
/home/jbourbeau/.local/lib/python2.7/site-packages/IPython/core/formatters.pyc in __call__(self, obj)
    305                 pass
    306             else:
--> 307                 return printer(obj)
    308             # Finally look for special method names
    309             method = get_real_method(obj, self.print_method)

/home/jbourbeau/.local/lib/python2.7/site-packages/IPython/core/pylabtools.pyc in <lambda>(fig)
    225 
    226     if 'png' in formats:
--> 227         png_formatter.for_type(Figure, lambda fig: print_figure(fig, 'png', **kwargs))
    228     if 'retina' in formats or 'png2x' in formats:
    229         png_formatter.for_type(Figure, lambda fig: retina_figure(fig, **kwargs))

/home/jbourbeau/.local/lib/python2.7/site-packages/IPython/core/pylabtools.pyc in print_figure(fig, fmt, bbox_inches, **kwargs)
    117 
    118     bytes_io = BytesIO()
--> 119     fig.canvas.print_figure(bytes_io, **kw)
    120     data = bytes_io.getvalue()
    121     if fmt == 'svg':

/home/jbourbeau/.local/lib/python2.7/site-packages/matplotlib/backend_bases.pyc in print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, **kwargs)
   2178                     orientation=orientation,
   2179                     dryrun=True,
-> 2180                     **kwargs)
   2181                 renderer = self.figure._cachedRenderer
   2182                 bbox_inches = self.figure.get_tightbbox(renderer)

/home/jbourbeau/.local/lib/python2.7/site-packages/matplotlib/backends/backend_agg.pyc in print_png(self, filename_or_obj, *args, **kwargs)
    525 
    526     def print_png(self, filename_or_obj, *args, **kwargs):
--> 527         FigureCanvasAgg.draw(self)
    528         renderer = self.get_renderer()
    529         original_dpi = renderer.dpi

/home/jbourbeau/.local/lib/python2.7/site-packages/matplotlib/backends/backend_agg.pyc in draw(self)
    472 
    473         try:
--> 474             self.figure.draw(self.renderer)
    475         finally:
    476             RendererAgg.lock.release()

/home/jbourbeau/.local/lib/python2.7/site-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs)
     59     def draw_wrapper(artist, renderer, *args, **kwargs):
     60         before(artist, renderer)
---> 61         draw(artist, renderer, *args, **kwargs)
     62         after(artist, renderer)
     63 

/home/jbourbeau/.local/lib/python2.7/site-packages/matplotlib/figure.pyc in draw(self, renderer)
   1157         dsu.sort(key=itemgetter(0))
   1158         for zorder, a, func, args in dsu:
-> 1159             func(*args)
   1160 
   1161         renderer.close_group('figure')

/home/jbourbeau/.local/lib/python2.7/site-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs)
     59     def draw_wrapper(artist, renderer, *args, **kwargs):
     60         before(artist, renderer)
---> 61         draw(artist, renderer, *args, **kwargs)
     62         after(artist, renderer)
     63 

/home/jbourbeau/.local/lib/python2.7/site-packages/matplotlib/axes/_base.pyc in draw(self, renderer, inframe)
   2322 
   2323         for zorder, a in dsu:
-> 2324             a.draw(renderer)
   2325 
   2326         renderer.close_group('axes')

/home/jbourbeau/.local/lib/python2.7/site-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs)
     59     def draw_wrapper(artist, renderer, *args, **kwargs):
     60         before(artist, renderer)
---> 61         draw(artist, renderer, *args, **kwargs)
     62         after(artist, renderer)
     63 

/home/jbourbeau/.local/lib/python2.7/site-packages/matplotlib/axis.pyc in draw(self, renderer, *args, **kwargs)
   1106         ticks_to_draw = self._update_ticks(renderer)
   1107         ticklabelBoxes, ticklabelBoxes2 = self._get_tick_bboxes(ticks_to_draw,
-> 1108                                                                 renderer)
   1109 
   1110         for tick in ticks_to_draw:

/home/jbourbeau/.local/lib/python2.7/site-packages/matplotlib/axis.pyc in _get_tick_bboxes(self, ticks, renderer)
   1056         for tick in ticks:
   1057             if tick.label1On and tick.label1.get_visible():
-> 1058                 extent = tick.label1.get_window_extent(renderer)
   1059                 ticklabelBoxes.append(extent)
   1060             if tick.label2On and tick.label2.get_visible():

/home/jbourbeau/.local/lib/python2.7/site-packages/matplotlib/text.pyc in get_window_extent(self, renderer, dpi)
    959             raise RuntimeError('Cannot get window extent w/o renderer')
    960 
--> 961         bbox, info, descent = self._get_layout(self._renderer)
    962         x, y = self.get_unitless_position()
    963         x, y = self.get_transform().transform_point((x, y))

/home/jbourbeau/.local/lib/python2.7/site-packages/matplotlib/text.pyc in _get_layout(self, renderer)
    359                 w, h, d = renderer.get_text_width_height_descent(clean_line,
    360                                                         self._fontproperties,
--> 361                                                         ismath=ismath)
    362             else:
    363                 w, h, d = 0, 0, 0

/home/jbourbeau/.local/lib/python2.7/site-packages/matplotlib/backends/backend_agg.pyc in get_text_width_height_descent(self, s, prop, ismath)
    227             fontsize = prop.get_size_in_points()
    228             w, h, d = texmanager.get_text_width_height_descent(s, fontsize,
--> 229                                                                renderer=self)
    230             return w, h, d
    231 

/home/jbourbeau/.local/lib/python2.7/site-packages/matplotlib/texmanager.pyc in get_text_width_height_descent(self, tex, fontsize, renderer)
    673         else:
    674             # use dviread. It sometimes returns a wrong descent.
--> 675             dvifile = self.make_dvi(tex, fontsize)
    676             dvi = dviread.Dvi(dvifile, 72 * dpi_fraction)
    677             try:

/home/jbourbeau/.local/lib/python2.7/site-packages/matplotlib/texmanager.pyc in make_dvi(self, tex, fontsize)
    420                      'string:\n%s\nHere is the full report generated by '
    421                      'LaTeX: \n\n' % repr(tex.encode('unicode_escape')) +
--> 422                      report))
    423             else:
    424                 mpl.verbose.report(report, 'debug')

RuntimeError: LaTeX was not able to process the following string:
'InIce_charge'
Here is the full report generated by LaTeX: 

This is pdfTeXk, Version 3.141592-1.40.3 (Web2C 7.5.6)
 %&-line parsing enabled.
entering extended mode
(./14e874c62740e610704d2634f5fc9462.tex
LaTeX2e <2005/12/01>
Babel <v3.8h> and hyphenation patterns for english, usenglishmax, dumylang, noh
yphenation, arabic, basque, bulgarian, coptic, welsh, czech, slovak, german, ng
erman, danish, esperanto, spanish, catalan, galician, estonian, farsi, finnish,
 french, greek, monogreek, ancientgreek, croatian, hungarian, interlingua, ibyc
us, indonesian, icelandic, italian, latin, mongolian, dutch, norsk, polish, por
tuguese, pinyin, romanian, russian, slovenian, uppersorbian, serbian, swedish, 
turkish, ukenglish, ukrainian, loaded.
(/usr/share/texmf/tex/latex/base/article.cls
Document Class: article 2005/09/16 v1.4f Standard LaTeX document class
(/usr/share/texmf/tex/latex/base/size10.clo))
(/usr/share/texmf/tex/latex/type1cm/type1cm.sty)
(/usr/share/texmf/tex/latex/psnfss/helvet.sty
(/usr/share/texmf/tex/latex/graphics/keyval.sty))
(/usr/share/texmf/tex/latex/psnfss/courier.sty)
(/usr/share/texmf/tex/latex/base/textcomp.sty
(/usr/share/texmf/tex/latex/base/ts1enc.def))
(/usr/share/texmf/tex/latex/geometry/geometry.sty

Package geometry Warning: Over-specification in `h'-direction.
    `width' (5058.9pt) is ignored.


Package geometry Warning: Over-specification in `v'-direction.
    `height' (5058.9pt) is ignored.

)
No file 14e874c62740e610704d2634f5fc9462.aux.
(/usr/share/texmf/tex/latex/base/ts1cmr.fd)
! Missing $ inserted.
<inserted text> 
                $
l.12 ...size{8.330000}{10.412500}{\rmfamily InIce_
                                                  charge}
! Extra }, or forgotten $.
l.12 ...330000}{10.412500}{\rmfamily InIce_charge}
                                                  
! Missing $ inserted.
<inserted text> 
                $
l.13 \end{document}
                   
[1] (./14e874c62740e610704d2634f5fc9462.aux) )
(\end occurred inside a group at level 1)

### simple group (level 1) entered at line 12 ({)
### bottom level
(see the transcript file for additional information)
Output written on 14e874c62740e610704d2634f5fc9462.dvi (1 page, 312 bytes).
Transcript written on 14e874c62740e610704d2634f5fc9462.log.






    






<matplotlib.figure.Figure at 0x4b13d90>

In [ ]: