In [1]:

    

%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
from matplotlib import rc

import numpy as np
import tables

import cascade as cas
import cascade_secs as csx

import earth


    

In [2]:

    

# gamma = 2 #spectral index of incoming neutrino flux
gamma = 'data/phiHGextrap.dat' #this is how you would specify an incoming flux from a file. Needs to be 200x1, on the same energy grid as below
flavor = 3 # 1 = nu_e, 2= nu_mu, 3= nu_tau. Negative for antineutrinos
Na = 6.0221415e23


    

In [3]:

    

def get_avg_attn(flavor,gamma,h5filename = "./data/NuFATECrossSections.h5"):
    w,v,ci,energy_nodes,phi_0 = cas.get_eigs(flavor,gamma,h5filename)
    tlength = 230
    tvec = np.linspace(-1,0,tlength)
    phiv = energy_nodes*0.
    for ctheta in tvec:
        t = earth.get_t_earth(np.arccos(ctheta))*Na # g/ cm^2
        phisol = np.dot(v,(ci*np.exp(w*t)))/phi_0
        phiv = phiv + phisol
    phiv = phiv/tlength
    return energy_nodes, phiv

def get_avg_attn_secs(flavor,gamma,h5filename = "./data/NuFATECrossSections.h5"):
    w,v,ci,energy_nodes,phi_0 = csx.get_eigs(flavor,gamma,h5filename)
    tlength = 230
    tvec = np.linspace(-1,0,tlength)
    phiv = phi_0*0.
    for ctheta in tvec:
        t = earth.get_t_earth(np.arccos(ctheta))*Na # g/ cm^2
        phisol = np.dot(v,(ci*np.exp(w*t)))/phi_0
        phiv = phiv + phisol
    phiv = phiv/tlength
    phiv = phiv[0:200] #only keep non-tau bit
    return energy_nodes, phiv


    

In [4]:

    

flavor = -3
zenith = np.radians(100.)
w,v,ci,energy_nodes,phi_0 = cas.get_eigs(flavor,gamma, "./data/NuFATECrossSections.h5")

t = earth.get_t_earth(zenith)*Na
phisol = np.dot(v,(ci*np.exp(w*t)))/phi_0


    

In [5]:

    

plt.figure(figsize=(6,5))
xsh5 = tables.open_file("./data/NuFATECrossSections.h5","r")
sigma_array = xsh5.root.total_cross_sections.nutaubarxs[:]
xsh5.close()
plt.semilogx(energy_nodes,np.exp(-sigma_array*t),c='b',lw=2)
plt.semilogx(energy_nodes,phisol,c='r',lw=2)

plt.legend(['Exponential suppression', 'Full cascade solution'],
          loc="upper right")
plt.xlabel(r"Neutrino Energy (GeV)")
plt.ylim(0.,1.)
plt.ylabel(r"Attenuation")
plt.grid()


    

In [6]:

    

import glob


    

In [7]:

    

xs_tables=glob.glob("./data/*.h5")
xs_tables_ct=glob.glob("./data/*HERA*.h5")
xs_tables_ct=glob.glob("./data/*CT*.h5")
xs_tables_nn=glob.glob("./data/*NN*.h5")


    

In [8]:

    

from matplotlib import rc, rcParams

rc('text',usetex=True)
rc('font',**{'family':'serif','serif':['Computer Modern'], 'size' : 18})
cols = ['#29A2C6','#FF6D31','#FFCB18','#73B66B','#EF597B', '#333333']

font = {'family' : 'serif',
        'weight' : 'bold',
        'size'   : 18}


    

In [10]:

    

energy_nodes, phim3 = get_avg_attn(-3,gamma)
energy_nodes, phim2 = get_avg_attn(-2,gamma)
energy_nodes, phim1 = get_avg_attn(-1,gamma)

energy_nodes, phi3 = get_avg_attn(3,gamma)
energy_nodes, phi2 = get_avg_attn(2,gamma)
energy_nodes, phi1 = get_avg_attn(1,gamma)


    

    




C:\ProgramData\Anaconda2\lib\site-packages\scipy\integrate\quadpack.py:364: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.
  If increasing the limit yields no improvement it is advised to analyze 
  the integrand in order to determine the difficulties.  If the position of a 
  local difficulty can be determined (singularity, discontinuity) one will 
  probably gain from splitting up the interval and calling the integrator 
  on the subranges.  Perhaps a special-purpose integrator should be used.
  warnings.warn(msg, IntegrationWarning)

In [22]:

    

plt.figure(figsize=(6,5))
plt.semilogx(energy_nodes,phim1,linestyle='--',c='m',label="NuEBar")
plt.semilogx(energy_nodes,phim2,linestyle='--',c='g',label="NuMuBar")
plt.semilogx(energy_nodes,phim3,linestyle='--',c='b',label="NuTauBar")
plt.semilogx(energy_nodes,phi1,c='r',label="NuE")
plt.semilogx(energy_nodes,phi2,c='g',label="NuMu")
plt.semilogx(energy_nodes,phi3,c='b',label="NuTau")
plt.xlim(1e3,1e10)
plt.ylim(0.0,1.)
plt.legend(loc="lower left")
plt.xlabel("Neutrino Energy (GeV)")
plt.ylabel(r"Attenuation")
plt.minorticks_on()
plt.grid()


    

In [11]:

    

gamma = 2
flavor = -1
w,v,ci,energy_nodes,phi_0 = cas.get_eigs(flavor,gamma,"./data/NuFATECrossSections.h5")
tlength = 100
Eindex = 80 #which energy index do you want to see?
d = 0 #at what depth do you want the detection point? default set to 0 km.

tvec = np.linspace(0,1,tlength)
phiv1 = tvec*0
phiv2 = tvec*0
phiv3 = tvec*0
phiv4 = tvec*0
phiv5 = tvec*0
for i in range(0,tlength):
    ctheta = tvec[i]
    t = earth.get_t_earth(np.arccos(-ctheta), d)*Na # g/ cm^2
    phisol = np.dot(v,(ci*np.exp(w*t)))*energy_nodes**(gamma-2.)
    phiv1[i] = phisol[0]
    phiv2[i] = phisol[25]
    phiv3[i] = phisol[50]
    phiv4[i] = phisol[75]
    phiv5[i] = phisol[99]

plt.figure(figsize=(6,5))
plt.plot(tvec,phiv1,lw=2)
plt.plot(tvec,phiv2,lw=2)
plt.plot(tvec,phiv3,lw=2)
plt.plot(tvec,phiv4,lw=2)
plt.plot(tvec,phiv5,lw=2)
plt.legend(['1 TeV','18 TeV','335 TeV','6.1 PeV','10 PeV'],loc="lower left")
plt.xlabel("-cos(zenith)")
plt.ylabel("Attenuation")
plt.title('NuEBar')
plt.minorticks_on()
plt.ylim(0.,1)
plt.grid()


    

    




---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
C:\ProgramData\Anaconda2\lib\site-packages\IPython\core\formatters.pyc in __call__(self, obj)
    332                 pass
    333             else:
--> 334                 return printer(obj)
    335             # Finally look for special method names
    336             method = get_real_method(obj, self.print_method)

C:\ProgramData\Anaconda2\lib\site-packages\IPython\core\pylabtools.pyc in <lambda>(fig)
    238 
    239     if 'png' in formats:
--> 240         png_formatter.for_type(Figure, lambda fig: print_figure(fig, 'png', **kwargs))
    241     if 'retina' in formats or 'png2x' in formats:
    242         png_formatter.for_type(Figure, lambda fig: retina_figure(fig, **kwargs))

C:\ProgramData\Anaconda2\lib\site-packages\IPython\core\pylabtools.pyc in print_figure(fig, fmt, bbox_inches, **kwargs)
    122 
    123     bytes_io = BytesIO()
--> 124     fig.canvas.print_figure(bytes_io, **kw)
    125     data = bytes_io.getvalue()
    126     if fmt == 'svg':

C:\ProgramData\Anaconda2\lib\site-packages\matplotlib\backend_bases.py in print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, **kwargs)
   2198                     orientation=orientation,
   2199                     dryrun=True,
-> 2200                     **kwargs)
   2201                 renderer = self.figure._cachedRenderer
   2202                 bbox_inches = self.figure.get_tightbbox(renderer)

C:\ProgramData\Anaconda2\lib\site-packages\matplotlib\backends\backend_agg.py in print_png(self, filename_or_obj, *args, **kwargs)
    543 
    544     def print_png(self, filename_or_obj, *args, **kwargs):
--> 545         FigureCanvasAgg.draw(self)
    546         renderer = self.get_renderer()
    547         original_dpi = renderer.dpi

C:\ProgramData\Anaconda2\lib\site-packages\matplotlib\backends\backend_agg.py in draw(self)
    462 
    463         try:
--> 464             self.figure.draw(self.renderer)
    465         finally:
    466             RendererAgg.lock.release()

C:\ProgramData\Anaconda2\lib\site-packages\matplotlib\artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
     61     def draw_wrapper(artist, renderer, *args, **kwargs):
     62         before(artist, renderer)
---> 63         draw(artist, renderer, *args, **kwargs)
     64         after(artist, renderer)
     65 

C:\ProgramData\Anaconda2\lib\site-packages\matplotlib\figure.py in draw(self, renderer)
   1142 
   1143             mimage._draw_list_compositing_images(
-> 1144                 renderer, self, dsu, self.suppressComposite)
   1145 
   1146             renderer.close_group('figure')

C:\ProgramData\Anaconda2\lib\site-packages\matplotlib\image.py in _draw_list_compositing_images(renderer, parent, dsu, suppress_composite)
    137     if not_composite or not has_images:
    138         for zorder, a in dsu:
--> 139             a.draw(renderer)
    140     else:
    141         # Composite any adjacent images together

C:\ProgramData\Anaconda2\lib\site-packages\matplotlib\artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
     61     def draw_wrapper(artist, renderer, *args, **kwargs):
     62         before(artist, renderer)
---> 63         draw(artist, renderer, *args, **kwargs)
     64         after(artist, renderer)
     65 

C:\ProgramData\Anaconda2\lib\site-packages\matplotlib\axes\_base.py in draw(self, renderer, inframe)
   2424             renderer.stop_rasterizing()
   2425 
-> 2426         mimage._draw_list_compositing_images(renderer, self, dsu)
   2427 
   2428         renderer.close_group('axes')

C:\ProgramData\Anaconda2\lib\site-packages\matplotlib\image.py in _draw_list_compositing_images(renderer, parent, dsu, suppress_composite)
    137     if not_composite or not has_images:
    138         for zorder, a in dsu:
--> 139             a.draw(renderer)
    140     else:
    141         # Composite any adjacent images together

C:\ProgramData\Anaconda2\lib\site-packages\matplotlib\artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
     61     def draw_wrapper(artist, renderer, *args, **kwargs):
     62         before(artist, renderer)
---> 63         draw(artist, renderer, *args, **kwargs)
     64         after(artist, renderer)
     65 

C:\ProgramData\Anaconda2\lib\site-packages\matplotlib\axis.py in draw(self, renderer, *args, **kwargs)
   1136         ticks_to_draw = self._update_ticks(renderer)
   1137         ticklabelBoxes, ticklabelBoxes2 = self._get_tick_bboxes(ticks_to_draw,
-> 1138                                                                 renderer)
   1139 
   1140         for tick in ticks_to_draw:

C:\ProgramData\Anaconda2\lib\site-packages\matplotlib\axis.py in _get_tick_bboxes(self, ticks, renderer)
   1076         for tick in ticks:
   1077             if tick.label1On and tick.label1.get_visible():
-> 1078                 extent = tick.label1.get_window_extent(renderer)
   1079                 ticklabelBoxes.append(extent)
   1080             if tick.label2On and tick.label2.get_visible():

C:\ProgramData\Anaconda2\lib\site-packages\matplotlib\text.py in get_window_extent(self, renderer, dpi)
    965             raise RuntimeError('Cannot get window extent w/o renderer')
    966 
--> 967         bbox, info, descent = self._get_layout(self._renderer)
    968         x, y = self.get_unitless_position()
    969         x, y = self.get_transform().transform_point((x, y))

C:\ProgramData\Anaconda2\lib\site-packages\matplotlib\text.py in _get_layout(self, renderer)
    351         tmp, lp_h, lp_bl = renderer.get_text_width_height_descent('lp',
    352                                                          self._fontproperties,
--> 353                                                          ismath=False)
    354         offsety = (lp_h - lp_bl) * self._linespacing
    355 

C:\ProgramData\Anaconda2\lib\site-packages\matplotlib\backends\backend_agg.py in get_text_width_height_descent(self, s, prop, ismath)
    228             fontsize = prop.get_size_in_points()
    229             w, h, d = texmanager.get_text_width_height_descent(s, fontsize,
--> 230                                                                renderer=self)
    231             return w, h, d
    232 

C:\ProgramData\Anaconda2\lib\site-packages\matplotlib\texmanager.py in get_text_width_height_descent(self, tex, fontsize, renderer)
    674         else:
    675             # use dviread. It sometimes returns a wrong descent.
--> 676             dvifile = self.make_dvi(tex, fontsize)
    677             dvi = dviread.Dvi(dvifile, 72 * dpi_fraction)
    678             try:

C:\ProgramData\Anaconda2\lib\site-packages\matplotlib\texmanager.py in make_dvi(self, tex, fontsize)
    421                      'string:\n%s\nHere is the full report generated by '
    422                      'LaTeX: \n\n' % repr(tex.encode('unicode_escape')) +
--> 423                      report))
    424             else:
    425                 mpl.verbose.report(report, 'debug')

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







    






<matplotlib.figure.Figure at 0xb16f4a8>

In [24]:

    

energy_nodes,phim2s = get_avg_attn_secs(-2,gamma)
energy_nodes,phim1s = get_avg_attn_secs(-1,gamma)
energy_nodes,phi1s = get_avg_attn_secs(1,gamma)
energy_nodes,phi2s = get_avg_attn_secs(2,gamma)


    

In [25]:

    

plt.figure(figsize=(6,5))
plt.semilogx(energy_nodes,phim1,linestyle='--',c='r')
plt.semilogx(energy_nodes,phim2,linestyle='--',c='g')
plt.semilogx(energy_nodes,phim3,linestyle='--',c='b')
plt.semilogx(energy_nodes,phi1,c='r')
plt.semilogx(energy_nodes,phi2,c='g')
plt.semilogx(energy_nodes,phi3,c='b')
plt.semilogx(energy_nodes,phi1s,c='r',linestyle='-.')
plt.semilogx(energy_nodes,phi2s,c='g',linestyle=':')
plt.semilogx(energy_nodes,phim1s,c='r',linestyle='-.')
plt.semilogx(energy_nodes,phim2s,c='g',linestyle=':')
plt.xlim(1e3,1e10)
plt.ylim(.0,1)
plt.ylabel("Attenuation")
plt.xlabel("Neutrino Energy (GeV)")
plt.minorticks_on()
plt.grid()


    

In [ ]:

    




    

In [ ]:

    




    

In [ ]: