From c14285f8a78e8df7436424bedd5976a29f8a1088 Mon Sep 17 00:00:00 2001
From: Nicholas Kamp <43788191+nickkamp1@users.noreply.github.com>
Date: Tue, 9 Jan 2024 18:40:54 -0500
Subject: [PATCH] Feature/dark news table (#37)

* reorganizing examples and earthparams files

* starting to add fill interpolation method

* Flush out fill table methods for LIDarkNews

* Minor cleanup
---
 python/LIDarkNews.py                          | 112 +++++++---
 resources/Examples/Example1/DIS_IceCube.py    |  66 ++++++
 .../Examples/Example2/DipolePortal_MINERvA.py |  82 ++++++++
 .../Example2/DipolePortal_MiniBooNE.py        |  37 +---
 .../Examples/Example2/VisualizeOutput.ipynb   | 192 ++++++++++++++++++
 resources/earthparams/densities/CCM.dat       |  47 -----
 resources/earthparams/densities/MINERvA.dat   |  56 +----
 resources/earthparams/densities/MiniBooNE.dat |  46 -----
 .../earthparams/densities/PREM_IceCube.dat    |  19 ++
 .../densities/{ => legacy}/PREM.dat           |   0
 .../densities/{ => legacy}/PREM_mmc.dat       |   0
 .../densities/{ => legacy}/mmc.dat            |   0
 resources/earthparams/materials/IceCube.dat   |  24 +++
 13 files changed, 479 insertions(+), 202 deletions(-)
 create mode 100644 resources/Examples/Example1/DIS_IceCube.py
 create mode 100644 resources/Examples/Example2/DipolePortal_MINERvA.py
 create mode 100644 resources/Examples/Example2/VisualizeOutput.ipynb
 create mode 100644 resources/earthparams/densities/PREM_IceCube.dat
 rename resources/earthparams/densities/{ => legacy}/PREM.dat (100%)
 rename resources/earthparams/densities/{ => legacy}/PREM_mmc.dat (100%)
 rename resources/earthparams/densities/{ => legacy}/mmc.dat (100%)
 create mode 100644 resources/earthparams/materials/IceCube.dat

diff --git a/python/LIDarkNews.py b/python/LIDarkNews.py
index 7f1c58f9..4a4f4ac8 100644
--- a/python/LIDarkNews.py
+++ b/python/LIDarkNews.py
@@ -97,9 +97,16 @@ def __init__(self,
             self.decays.append(PyDarkNewsDecay(dec_case,
                                                table_dir = self.table_dir + table_subdirs))
             
-    def SaveCrossSectionTables(self):
+    def SaveCrossSectionTables(self,fill_tables_at_exit=True):
+        if not fill_tables_at_exit:
+            print('WARNING: Saving tables without filling PyDarkNewsCrossSection interpolation tables. Future updates to DarkNews can lead to inconsistent behavior if new entries are ever added to this table')
         for cross_section in self.cross_sections:
+            if fill_tables_at_exit:
+                print('Filling cross section table at %s'%cross_section.table_dir)
+                num = cross_section.FillInterpolationTables()
+                print('Added %d points'%num)
             cross_section.SaveInterpolationTables()
+            
         
 
 
@@ -112,7 +119,7 @@ def __init__(self,
                  table_dir=None, # table to store 
                  tolerance=1e-6, # supposed to represent machine epsilon
                  interp_tolerance=5e-2, # relative interpolation tolerance
-                 interpolate_differential = False):
+                 interpolate_differential=False):
         
         DarkNewsCrossSection.__init__(self) # C++ constructor
         
@@ -170,23 +177,22 @@ def _redefine_interpolation_objects(self,total=False,diff=False):
                 # If we only have two energy points, don't try to construct interpolator
                 if len(np.unique(self.differential_cross_section_table[:,0])) <= 2: return
                 self.differential_cross_section_interpolator = CloughTocher2DInterpolator(self.differential_cross_section_table[:,:2],
-                                                                                        self.differential_cross_section_table[:,2],
-                                                                                        rescale=True)
+                                                                                          self.differential_cross_section_table[:,2],
+                                                                                          rescale=True)
 
     # Check whether we have close-enough entries in the intrepolation tables
-    def _query_interpolation_table(self,inputs,mode):
-        # 
-        # returns:
-        # 0 if we are not close enough to any points in the interpolation table
-        # otherwise, returns the desired interpolated value
-        
+    def _interpolation_flags(self,inputs,mode):
+        #
+        # returns UseSinglePoint,Interpolate,closest_idx
+        # UseSinglePoint: whether to use a single point in table
+        # Interpolate: whether to interpolate bewteen different points
+        # closest_idx: index of closest point in table (for UseSinglePoint)
+
         # Determine which table we are using
         if mode=='total':
             interp_table = self.total_cross_section_table
-            interpolator = self.total_cross_section_interpolator
         elif mode=='differential':
             interp_table = self.differential_cross_section_table
-            interpolator = self.differential_cross_section_interpolator
         else:
             print('Invalid interpolation table mode %s'%mode)
             exit(0)
@@ -197,20 +203,15 @@ def _query_interpolation_table(self,inputs,mode):
         # bools to keep track of whether to use a single point or interpolate
         UseSinglePoint = True
         Interpolate = True
-        prev_closest_idx = None
+        # First check whether we have a close-enough single point
+        closest_idx = np.argmin(np.sum(np.abs(interp_table[:,:-1] - inputs)))
+        diff = (interp_table[closest_idx,:-1] - inputs)/inputs
+        if np.all(np.abs(diff)<self.tolerance): 
+            UseSinglePoint = True
         for i,input in enumerate(inputs):
             closest_idx = np.argmin(np.abs(interp_table[:,i] - input))
             diff = (interp_table[closest_idx,i] - input)/input # relative difference
-            # First check if we are close enough to use a single point
-            if np.abs(diff) >= self.tolerance:
-                # We are not close enough to use one existing table point
-                UseSinglePoint = False
-            elif UseSinglePoint:
-                # Check whether the closest_idx found here is the same as the previous
-                if i>0 and closest_idx != prev_closest_idx:
-                    UseSinglePoint = False
-                prev_closest_idx = closest_idx
-            # Now check if we are close enough to interpolate
+            # Check if we are close enough to interpolate
             if np.abs(diff) >= self.interp_tolerance:
                 Interpolate = False
             else:
@@ -239,6 +240,28 @@ def _query_interpolation_table(self,inputs,mode):
                         # check if the node above is also within the interpolation tolerance
                         if diff_above>=self.interp_tolerance:
                             Interpolate = False
+        return UseSinglePoint, Interpolate, closest_idx
+    
+    # return entries in interpolation table if we have inputs 
+    def _query_interpolation_table(self,inputs,mode):
+        # 
+        # returns:
+        # 0 if we are not close enough to any points in the interpolation table
+        # otherwise, returns the desired interpolated value
+        
+        # Determine which table we are using
+        if mode=='total':
+            interp_table = self.total_cross_section_table
+            interpolator = self.total_cross_section_interpolator
+        elif mode=='differential':
+            interp_table = self.differential_cross_section_table
+            interpolator = self.differential_cross_section_interpolator
+        else:
+            print('Invalid interpolation table mode %s'%mode)
+            exit(0)
+
+        UseSinglePoint, Interpolate, closest_idx = self._interpolation_flags(inputs,mode)
+
         if UseSinglePoint:
             return interp_table[closest_idx,-1]
         elif Interpolate:
@@ -246,6 +269,49 @@ def _query_interpolation_table(self,inputs,mode):
         else:
             return 0    
     
+    # Fills the total and differential cross section tables within interp_tolerance
+    def FillInterpolationTables(self,total=True,diff=True):
+        Emin = (1.0+self.tolerance)*self.ups_case.Ethreshold
+        Emax = np.max(self.total_cross_section_table[:,0])
+        num_added_points = 0
+        if total:
+            E = Emin
+            while E<Emax:
+                UseSinglePoint,Interpolate,_ = self._interpolation_flags([E],'total')
+                if not (UseSinglePoint or Interpolate):
+                    xsec = self.ups_case.scalar_total_xsec(E)
+                    self.total_cross_section_table = np.append(self.total_cross_section_table,[[E,xsec]],axis=0)
+                    num_added_points+=1
+                E *= (1+self.interp_tolerance)
+        if diff:
+            # interaction record to calculate Q2 bounds
+            interaction = LI.dataclasses.InteractionRecord()
+            interaction.signature.primary_type = self.GetPossiblePrimaries()[0] # only one primary
+            interaction.signature.target_type = self.GetPossibleTargets()[0] # only one target
+            interaction.target_mass = self.ups_case.MA
+            E = Emin
+            zmin,zmax = self.tolerance,1
+            z = zmin
+            while E < Emax:
+                interaction.primary_momentum = [E,0,0,0]
+                Q2min = self.Q2Min(interaction)
+                Q2max = self.Q2Max(interaction)
+                z = zmin
+                while z < zmax:
+                    Q2 = Q2min + z*(Q2max-Q2min)
+                    UseSinglePoint,Interpolate,_ = self._interpolation_flags([E,z],'differential')
+                    if not (UseSinglePoint or Interpolate):
+                        dxsec = self.ups_case.diff_xsec_Q2(E, Q2).item()
+                        self.differential_cross_section_table = np.append(self.differential_cross_section_table,[[E,z,dxsec]],axis=0)
+                        num_added_points+=1
+                    z *= (1+self.interp_tolerance)
+                E *= (1+self.interp_tolerance)
+        self._redefine_interpolation_objects(total=total,diff=diff)
+        return num_added_points
+
+
+
+    
     # Saves the tables for the scipy interpolation objects
     def SaveInterpolationTables(self,total=True,diff=True):
         if total:
diff --git a/resources/Examples/Example1/DIS_IceCube.py b/resources/Examples/Example1/DIS_IceCube.py
new file mode 100644
index 00000000..924e6d98
--- /dev/null
+++ b/resources/Examples/Example1/DIS_IceCube.py
@@ -0,0 +1,66 @@
+import numpy as np
+import os
+
+import leptoninjector as LI
+import sys
+sys.path.insert(1,'/n/holylfs05/LABS/arguelles_delgado_lab/Everyone/nkamp/LIV2/sources/LeptonInjector/python')
+from LIController import LIController
+ 
+LI_SRC = os.environ.get('LEPTONINJECTOR_SRC')
+
+# Number of events to inject
+events_to_inject = 1000
+
+# Expeirment to run
+experiment = 'IceCube'
+
+# Define the controller
+controller = LIController(events_to_inject,
+                          experiment)
+
+# Particle to inject
+primary_type = LI.dataclasses.Particle.ParticleType.NuMu
+
+# Cross Section Model
+xsfiledir = LI_SRC+'resources/CrossSectionTables/DISSplines/'
+target_type = LI.dataclasses.Particle.ParticleType.Nucleon
+
+DIS_xs = LI.crosssections.DISFromSpline(xsfiledir+'test_xs.fits',
+                                        xsfiledir+'test_xs_total.fits',
+                                        [primary_type],
+                                        [target_type])
+
+primary_xs = LI.crosssections.CrossSectionCollection(primary_type, [DIS_xs])
+controller.SetCrossSections(primary_xs)
+
+# Primary distributions
+primary_injection_distributions = {}
+primary_physical_distributions = {}
+
+# energy distribution
+edist = LI.distributions.PowerLaw(2,1e3,1e6)
+primary_injection_distributions['energy'] = edist
+primary_physical_distributions['energy'] = edist
+
+# direction distribution
+direction_distribution = LI.distributions.IsotropicDirection()
+primary_injection_distributions['direction'] = direction_distribution
+primary_physical_distributions['direction'] = direction_distribution
+
+# position distribution
+muon_range_func = LI.distributions.LeptonDepthFunction()
+position_distribution = LI.distributions.ColumnDepthPositionDistribution(600, 600.,
+                                                                         muon_range_func,
+                                                                         set(controller.GetEarthModelTargets()[0]))
+primary_injection_distributions['position'] = position_distribution
+
+# SetProcesses
+controller.SetProcesses(primary_type,
+                        primary_injection_distributions,
+                        primary_physical_distributions)
+
+controller.Initialize()
+
+events = controller.GenerateEvents()
+
+controller.SaveEvents('output/MINERvA_Dipole_M%2.2f_mu%2.2e_example.hdf5'%(model_kwargs['m4'],model_kwargs['mu_tr_mu4']))
\ No newline at end of file
diff --git a/resources/Examples/Example2/DipolePortal_MINERvA.py b/resources/Examples/Example2/DipolePortal_MINERvA.py
new file mode 100644
index 00000000..b13f7fbf
--- /dev/null
+++ b/resources/Examples/Example2/DipolePortal_MINERvA.py
@@ -0,0 +1,82 @@
+import numpy as np
+import os
+
+import leptoninjector as LI
+import sys
+sys.path.insert(1,'/n/holylfs05/LABS/arguelles_delgado_lab/Everyone/nkamp/LIV2/sources/LeptonInjector/python')
+from LIController import LIController
+ 
+LI_SRC = os.environ.get('LEPTONINJECTOR_SRC')
+
+# Define a DarkNews model 
+model_kwargs = {
+    'm4': 0.47,#0.140,
+    'mu_tr_mu4': 2.5e-6, #1e-6, # GeV^-1
+    'UD4': 0,
+    'Umu4': 0,
+    'epsilon': 0.0,
+    'gD': 0.0,
+    'decay_product':'photon',
+    'noHC':True,
+    'HNLtype':"dirac",
+}
+
+# Number of events to inject
+events_to_inject = 1000
+
+# Expeirment to run
+experiment = 'MINERvA'
+
+# Define the controller
+controller = LIController(events_to_inject,
+                          experiment)
+
+# Particle to inject
+primary_type = LI.dataclasses.Particle.ParticleType.NuMu
+
+# Define DarkNews Model
+table_dir = LI_SRC+'/resources/CrossSectionTables/DarkNewsTables/Dipole_M%2.2f_mu%2.2e/'%(model_kwargs['m4'],model_kwargs['mu_tr_mu4'])
+controller.InputDarkNewsModel(primary_type,
+                              table_dir,
+                              model_kwargs)
+
+# Primary distributions
+primary_injection_distributions = {}
+primary_physical_distributions = {}
+
+# energy distribution
+flux_file = LI_SRC + '/resources/Fluxes/NUMI_Flux_Tables/ME_FHC_numu.txt'
+edist = LI.distributions.TabulatedFluxDistribution(flux_file, True)
+edist_gen = LI.distributions.TabulatedFluxDistribution(1.05*model_kwargs['m4'], 10, flux_file, False)
+primary_injection_distributions['energy'] = edist_gen
+primary_physical_distributions['energy'] = edist
+
+# Flux normalization: go from cm^-2 to m^-2
+flux_units = LI.distributions.NormalizationConstant(1e4)
+primary_physical_distributions['normalization'] = flux_units
+
+# direction distribution
+direction_distribution = LI.distributions.FixedDirection(LI.math.Vector3D(0, 0, 1.0))
+primary_injection_distributions['direction'] = direction_distribution
+primary_physical_distributions['direction'] = direction_distribution
+
+# position distribution
+decay_range_func = LI.distributions.DecayRangeFunction(model_kwargs['m4'],
+                                                       controller.DN_min_decay_width,
+                                                       3,
+                                                       240)
+position_distribution = LI.distributions.RangePositionDistribution(1.24, 5.,
+                                                                   decay_range_func,
+                                                                   set(controller.GetEarthModelTargets()[0]))
+primary_injection_distributions['position'] = position_distribution
+
+# SetProcesses
+controller.SetProcesses(primary_type,
+                        primary_injection_distributions,
+                        primary_physical_distributions)
+
+controller.Initialize()
+
+events = controller.GenerateEvents()
+
+controller.SaveEvents('output/MINERvA_Dipole_M%2.2f_mu%2.2e_example.hdf5'%(model_kwargs['m4'],model_kwargs['mu_tr_mu4']))
\ No newline at end of file
diff --git a/resources/Examples/Example2/DipolePortal_MiniBooNE.py b/resources/Examples/Example2/DipolePortal_MiniBooNE.py
index a30669c0..d175b5ab 100644
--- a/resources/Examples/Example2/DipolePortal_MiniBooNE.py
+++ b/resources/Examples/Example2/DipolePortal_MiniBooNE.py
@@ -23,7 +23,7 @@
 }
 
 # Number of events to inject
-events_to_inject = 10
+events_to_inject = 1000
 
 # Expeirment to run
 experiment = 'MiniBooNE'
@@ -80,38 +80,5 @@
 
 events = controller.GenerateEvents()
 
-controller.SaveEvents('Dipole_M%2.2f_mu%2.2e_example.hdf5'%(model_kwargs['m4'],model_kwargs['mu_tr_mu4']))
-
-
-# Analysis
-
-upscattering_vertices = {}
-targets = []
-decay_vertices = []
-for tree in events:
-    for datum in tree.tree:
-        if datum.record.signature.primary_type == LI.dataclasses.Particle.ParticleType.NuMu:
-            target = datum.record.signature.target_type
-            if target not in upscattering_vertices.keys():
-                upscattering_vertices[target] = []
-            upscattering_vertices[target].append(datum.record.interaction_vertex)
-        else:
-            decay_vertices.append(datum.record.interaction_vertex)
-
-for k,vertices in upscattering_vertices.items():
-    vertices = np.array(vertices)
-    plt.scatter(vertices[:,2],vertices[:,0],label=k,alpha=0.3)
-plt.xlabel('Z [m]')
-plt.ylabel('X [m]')
-plt.legend()
-plt.show()
-plt.savefig('Figures/%s_upscattering.png'%experiment,dpi=100)
-
-decay_vertices = np.array(decay_vertices)
-plt.scatter(decay_vertices[:,2],decay_vertices[:,0],alpha=0.3)
-plt.xlabel('Z [m]')
-plt.ylabel('X [m]')
-plt.legend()
-plt.show()
-plt.savefig('Figures/%s_decay.png'%experiment,dpi=100)
+controller.SaveEvents('output/MiniBooNE_Dipole_M%2.2f_mu%2.2e_example.hdf5'%(model_kwargs['m4'],model_kwargs['mu_tr_mu4']))
 
diff --git a/resources/Examples/Example2/VisualizeOutput.ipynb b/resources/Examples/Example2/VisualizeOutput.ipynb
new file mode 100644
index 00000000..a3d86f07
--- /dev/null
+++ b/resources/Examples/Example2/VisualizeOutput.ipynb
@@ -0,0 +1,192 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import h5py\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "\n",
+    "POT = 1.875e+21"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 26,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "f = h5py.File('output/MINERvA_Dipole_M0.47_mu2.50e-06_example.hdf5','r')\n",
+    "four_momenta = {'nu':[],\n",
+    "                'HNL':[],\n",
+    "                'Gamma':[],\n",
+    "                'NuOut':[]}\n",
+    "vertex = {'upscattering':[],\n",
+    "          'decay':[]}\n",
+    "weights = []\n",
+    "for i in range(f.attrs['num_events']):\n",
+    "    event_group = f['event%d'%i]\n",
+    "    weights.append(event_group.attrs['event_weight'])\n",
+    "    for j in range(event_group.attrs['num_interactions']):\n",
+    "        int_group = event_group['interaction%d'%j]\n",
+    "        if int_group.attrs['primary_type'] == 'ParticleType.NuMu':\n",
+    "            four_momenta['nu'].append(list(int_group['primary_momentum']))\n",
+    "            vertex['upscattering'].append(list(int_group['vertex']))\n",
+    "        else:\n",
+    "            four_momenta['HNL'].append(list(int_group['primary_momentum']))\n",
+    "            vertex['decay'].append(list(int_group['vertex']))\n",
+    "            for isec in range(2):\n",
+    "                if int_group.attrs['secondary_type%d'%isec] == 'ParticleType.Gamma':\n",
+    "                    four_momenta['Gamma'].append(list(int_group['secondary_momentum%d'%isec]))\n",
+    "                else:\n",
+    "                    four_momenta['NuOut'].append(list(int_group['secondary_momentum%d'%isec]))\n",
+    "\n",
+    "for k in four_momenta.keys():\n",
+    "    four_momenta[k] = np.array(four_momenta[k])\n",
+    "for k in vertex.keys():\n",
+    "    vertex[k] = np.array(vertex[k])\n",
+    "weights = np.where(np.isnan(weights),0,weights)\n",
+    "\n",
+    "f.close()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 27,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl0AAAGdCAYAAAAogsYCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA02klEQVR4nO3dB3gUdf7H8W8gJNQQSiAgvSid0IsNJAKCBY9+gMBRpJ5IEzwEBU8QpBeRQ4oKUk5FD5De/kdVygkIHEiXEhQh0kv2/3x/d7O3GxKKJr/sZt+v55lnszO/nczsJJlPfm2DXC6XSwAAAJCs0iTv7gEAAEDoAgAAsISaLgAAAAsIXQAAABYQugAAACwgdAEAAFhA6AIAALCA0AUAAGBBsI1vEiji4uLk1KlTkiVLFgkKCkrpwwEAAPdB54n/9ddfJW/evJImTfLVRxG6kpAGrvz58yflLgEAgCUnTpyQfPnyJdv+CV1JSGu4nIsWFhaWlLsGAADJJDY21lSaOPfx5ELoSkJOk6IGLkIXAAD+JSiZuwbRkR4AAMACQhcAAIAFhC4AAAAL6NMFAEASTj1w69YtuX37Nu+pD0mbNq0EBwen+HROhC4AAJLAjRs35PTp03LlyhXeTx+UMWNGyZMnj4SEhKTYMRC6AABIgsmxjxw5YmpUdIJNvbGndK0K/lf7qIH43Llz5hoVL148WSdAvRtCFwAAv5Pe1DV46VxPWqMC35IhQwZJly6dHDt2zFyr9OnTp8hx0JEeAICkuqmmUA0K/OPapPwRAAAABABCFwAAgAX06QIAIBmNXflva+/vq08/bO174cFR0wUAAGABoQsAAMACmhcBeJmya0qyvyPdorrxrgM+olatWlKuXDkzjcL06dPNHGNdunSRN998U44ePSqFCxeWnTt3SlRUlCl/4cIFyZYtm6xdu9a8FvePmi4AAALc7NmzJVOmTLJ161YZOXKkDB06VFauXJnSh5XqUNMFAECA05quIUOGmK91xvZJkybJ6tWrzddIOtR0AQAQ4DR0edLPKIyJiUmx40mtCF0AAAQ4/YgcT/q5kfqxRs4s7vr5hY6bN29aP77UgtAFAAASFBERYR5Pnz7tXrdr1y7erd+IPl0AACDRD4quXr26jBgxwoxi1CbHQYMG8W79RoQuAACSkb/PEj9jxgzp0KGDVKpUSR555BEzurFu3bopfVh+idAFAEAAW7du3R3rFi1a5P66ZMmSsmnTJq/tnn28cP/o0wUAAGABoQsAAMACQhcAAIAFhC4AAAALCF0AAACBELr0U8x15lvPpUSJEu7t165dk+7du0uOHDkkc+bM0rhxYzl79qzXPo4fPy4NGzaUjBkzSq5cuaRfv35y69atO0ZnVKxYUUJDQ6VYsWIya9asO45l8uTJUqhQIfNJ69WqVZNt27Yl45kDAIBA4hNTRpQuXVpWrVrlfh4c/L/DevXVV2XJkiWycOFCyZo1q/To0UP+8Ic/yMaNG83227dvm8AVGRlphrTqrLkvvfSS+UiDd955x5Q5cuSIKdOlSxeZM2eO+RDPjh07ms+Wqlevnikzf/586d27t0ydOtUErnHjxpltBw4cMEEOABIyZdeUZH9jukV1480HUoEUr+lyQpaGJmfJmTOnWX/x4kX58MMPZcyYMfLUU0+ZidlmzpxpwtWWLVtMmRUrVsj3338vn3zyiURFRckzzzwjw4YNM7VWN27cMGU0SOlMuqNHjzbzjWhwa9KkiYwdO9Z9DPo9OnXqJO3bt5dSpUqZ12jNmU4KBwAAkCpC18GDByVv3rxSpEgRadWqlWkuVNu3bzcfrBkdHe0uq02PBQoUkM2bN5vn+li2bFnJnTu3u4zWUMXGxsrevXvdZTz34ZRx9qHhTL+XZxn9kE997pRJyPXr18338VwAAAB8snlRm/K0f5V+tIA2Db711lvy+OOPy549e+TMmTMSEhIi4eHhXq/RgKXblD56Bi5nu7PtbmU0JF29elV++eUX00yZUJn9+/cneuzDhw83xwsAQKLWDrf35tQe+MAvadeunVy4cMFrFnqnL3Tt2rXNPVI/5Fq/1pag7777TtKmTesup/do7ZKj+1HaN7pXr15mgY/VdGlzYNOmTaVcuXKm9mnp0qXm4i9YsEB83cCBA00TqLOcOHEipQ8JAIBkc/jwYfnoo494h/01dMWnifnhhx+WQ4cOmf5d2vSnIcyTjl7UbUof449mdJ7fq0xYWJj5BHXtQ6apPaEyzj4SoiMhdR+eCwAAqVXPnj1lyJAhpnsNUkHounTpkvzwww9mZKF2nNdRiDra0KGjCbXPV40aNcxzfdy9e7fExMS4y6xcudIEIK0Gdcp47sMp4+xDmzD1e3mWiYuLM8+dMgAABDptMtQpmSZOnJjSh+KXUjx09e3bV9avXy9Hjx41oxJffPFFU+vUsmVLM0VEhw4dzFQOa9euNZ3ddXShBqHq1aub19etW9eEqzZt2si//vUvWb58uQwaNMjM7aU1UUqnitAq0f79+5s+WlOmTDHNlzodhUO/x9/+9jeZPXu27Nu3T7p27SqXL1823w8AgNRs8eLFZi5Mz0W7/8Sno/q1pkv7NGu3GvhZ6Dp58qQJWNqRvlmzZmYSVJ0OIiIiwmzXaR2effZZMynqE088YZr7Pv/8c/frNaDpD4s+ahhr3bq1madr6NCh7jI6XYTO9aW1W+XLlzdTR0yfPt09R5dq3ry5vPfeezJ48GAz9YR2Gly2bNkdnesBAEhttJO83vc8F71PJkQrQ/Re/e6771o/Tn+X4qMX582bd9ftOju8zrmlS2IKFixoOuDfTa1atWTnzp13LaPzd+kCAEAgyZQpk/m0lviVIonNrfnXv/7VjFbknulnNV0AAMC/6KwD+mkyTJvkZzVdAADA/4wYMcKrm46nH3/80TRRxm+VypYtmwQyaroAAMAD04/n00VHM8anfaQrVKjgtSxZsiTg32VqugAA8LFZ4m3ST4VJrC+0y+W642tPOmNAfDobARJGTRcAAIAFhC4AAAALCF0AAAAWELoAAAAsIHQBAABYQOgCAACwgNAFAABgAaELAADAAkIXAACABYQuAAAAC/gYIAAAktGUXVOsvb/dorr9ptedOXNGhg8fbj4f8eTJk5I1a1YpVqyYtG7dWtq2bSsZM2ZM8mMNRIQuAAAC2OHDh+XRRx+V8PBweeedd6Rs2bISGhoqu3fvlmnTpslDDz0kzz//fEofZqpA8yIAAAGsW7duEhwcLN9++600a9ZMSpYsKUWKFJEXXnjB1Hw999xzptyYMWNMIMuUKZPkz5/fvO7SpUteH5ytwW3x4sXyyCOPmNqxJk2ayJUrV2T27NlSqFAhyZYtm/z5z3+W27dvu1+n699++2156aWXJHPmzFKwYEH56quv5Ny5c+YYdF25cuXM8Tl+/vlnadmypQmE+n30uD799FPxdYQuAAAClIaXFStWSPfu3U2YSkhQUJB5TJMmjUyYMEH27t1rQtSaNWukf//+XmU1YGmZefPmybJly2TdunXy4osvytKlS83y8ccfywcffCB///vfvV43duxYU9u2c+dOadiwobRp08aEMG3e3LFjhxQtWtQ8d7lcpvy1a9ekUqVKJhTu2bNHOnfubF6zbds28WWELgAAAtShQ4dMkNGaKU85c+Y0NUy6vPbaa2Zdr169pHbt2qZm6qmnnjK1UwsWLPB63c2bN+X999+XChUqyBNPPGFquv75z3/Khx9+KKVKlZJnn33W7GPt2rVer2vQoIG8/PLLUrx4cRk8eLDExsZKlSpVpGnTpvLwww+bY9i3b5+cPXvWlNcarr59+0pUVJSplevZs6fUr1//juPxNfTpAgAAXrTGKC4uTlq1aiXXr18361atWmU62+/fv9+Eolu3bpkaJ63dcjra66PWSjly585tQpqGN891MTEx4kmbDz23K20yjL9OXxcZGWmaJ7X/mYasH3/8UW7cuGGO09c7/FPTBQBAgNIRitp8eODAAa/1Wnuk2zJkyGCeHz161NRSaTj67LPPZPv27TJ58mSzTQOPI126dF770X0ntC4uLs5rnWcZpzkzoXXO60aNGiXjx483NWBaa7Zr1y6pV6+e17H4IkIXAAABKkeOHPL000/LpEmT5PLly4mW05ClgWf06NFSvXp10+R36tQpSSkbN240ney1z1f58uVNSPz3v/8tvo7QBQBAAJsyZYppKqxcubLMnz/f9J3Smq9PPvnENCWmTZvW1Hppf62JEyeaKSa0Q/zUqVNT7JiLFy8uK1eulE2bNpnj1f5gTn8vX0boAgAggGkfLB01GB0dLQMHDjQ1RxrANGBpZ/Vhw4aZdTplxLvvvitlypSROXPmmP5dKWXQoEFSsWJF06RYq1Yt08+rUaNG4uuCXM74S/xu2rFQZ/G9ePGihIWF8Y7CL9mYPfu3zprti3i/oLRD+ZEjR6Rw4cKSPn163hQ/u0axlu7f1HQBAABYQOgCAACwgNAFAABgAaELAADAAkIXAACABYQuAAAACwhdAAAAFhC6AAAALCB0AQAAWEDoAgAAsCDYxjcBACBQnZs4ydr3iujZ44Ff065dO5k9e7b5LMUBAwa41y9atEhefPFFeZBPCzx//rwMHTpUvvjiCzl9+rTkzJlT6tevL2+++aYUKFDggY4rKCjI7McfPlPxflHTBQBAgNPPItQPs/7ll19+8z40cFWvXl1WrVolU6dOlUOHDsm8efPMY5UqVeTw4cMS6AhdAAAEuOjoaImMjDS1XQnRmqqoqCivdePGjZNChQq5n//lL3+RU6dOmdD1zDPPmJqtJ554QpYvXy7p0qWT7t27u8vq6/T1nnT/+n2c7Upr2rTGy/P7+DNCFwAAAS5t2rTyzjvvyMSJE+XkyZMP/Pq4uDhTq9WqVSsT3jxlyJBBunXrZsKX1obdj2+++cY8zpw50zRTOs/9XYqHLk3VWu2YJUsWyZUrl2m7PXDggFeZWrVqmaTruXTp0sWrzPHjx6Vhw4aSMWNGs59+/frJrVu3vMqsW7dOKlasKKGhoVKsWDGZNWvWHcczefJkk6i1qrVatWqybdu2ZDpzAAB8h9YqaW3TkCFDHvi1586dkwsXLkjJkiUT3K7rXS6XaWq8HxEREeYxPDzchDjnub9L8dC1fv16U+W4ZcsWWblypdy8eVPq1q0rly9f9irXqVMnk3adZeTIke5tt2/fNoHrxo0bsmnTJtMhUAPV4MGD3WWOHDliytSuXVt27dolvXr1ko4dO5rk7Zg/f7707t3b/MDt2LFDypcvL/Xq1ZOYmBhL7wYAAClH+3XpPXTfvn2/6fUP0uk+EKV46Fq2bJkZOVG6dGkTcjQsaa3V9u3bvcppDZamXWcJCwtzb1uxYoV8//338sknn5iUrm3Jw4YNM7VWGsSUduorXLiwjB492iTuHj16SJMmTWTs2LHu/YwZM8aEu/bt20upUqXMa/T7zpgxw+I7AgBAytA+WFrZMHDgQK/1adKkuSNQaSWJQ2uitFYqsbCm64OCgkwr0/3sL7VK8dAV38WLF81j9uzZvdbPmTPHDD0tU6aM+WG4cuWKe9vmzZulbNmykjt3bvc6/aGJjY2VvXv3ustoR0FPWkbXKw1nGvQ8y+gPhT53ysR3/fp18z08FwAA/NmIESPkH//4h9e9T0PVmTNnvIKSthp53i+bNWsmc+fONeU8Xb16VaZMmWLuuc69XfenrVYOvX9qi5Qn7XyvLVmpiU+FLu2Ip81+jz76qAlXjj/+8Y+mFmvt2rUmcH388cfSunVr93a9wJ6BSznPnYufWBm90PoD8dNPP5mLm1CZ+D9Anv3RsmbN6l7y58+fBO8CAAApRysxtEP8hAkTvPpWa78t7drzww8/mJakr7/+2ut12hFfW6Kefvpps+3EiROyYcMGE7a0Fmvy5Mnusk899ZS5l//f//2f7N69W9q2bWs683vS/tWrV6829+DfM5WFL/Gp0KV9u/bs2WNGQHjq3LmzuWjOD8JHH31kJkzTC5+SNABqzZyz6A8YAAD+Tic41YoQh3bL0doqDU7aFUgHmfXt29frNTly5DD9s7Xv9MsvvyxFixY1tV/6qKMPixQp4nX/fPLJJ+XZZ581/a11EJ2W86TdgbSvt1ZoVKhQQVIDn5mRXvtYLV682KTifPny3bWsjipUOgpCL5Im6/ijDM+ePWsenaGr+uis8yyjfcN0OKsmbF0SKhN/+KtDR0HqAgBAUs4Sb1NCI/m1lkm70HjSWQPizxzw+uuvez3XbkBaQ+ZZS5aQsLCwOypYtLbL03PPPWeW1CTFa7q0fVgDl9ZcrVmzxnR2vxenHTlPnjzmsUaNGqZ60nOUoaZjvajaId4po9WUnrSMrlchISFSqVIlrzKa8vW5UwYAAMBva7q0SVE73n355Zdmri6n/5T2kdIaKG1C1O0NGjQwVZffffedvPrqq2aERbly5UxZnWJCw1WbNm1Me7PuY9CgQWbfTk2UpvNJkyZJ//795U9/+pMJeAsWLJAlS5a4j0Wni9CkXblyZalataqZLVenrtDRjAAAAH4dut5//313Jz1POgutTiWhNVD6kQJOANK23caNG5tQ5dBmQW2a7Nq1q6mVypQpkwlP2ibt0Bo0DVga2MaPH2+aMKdPn276ijmaN29uOgrq/F4a3HT6CZ3SIn7negAAAL8LXfeaSE1Dlk6gei8FCxaUpUuX3rWMBrudO3fetYw2deoCAACQqvp0AQAABAJCFwAASYSPwfFdLh/4iCJCFwAAv5POnq48Py0FvuXKf6+Nc60Csk8XAAD+Tgd06WcPOlMX6ef26mcNwjdquK5cuWKujV6j+DPf20ToAgAgCTgTaXvOGQnfER4enuhk57YQugAASAJas6WTdufKlct81iB8R7p06VK0hstB6AIAIAk5HysHxEdHegAAAAsIXQAAABYQugAAACygTxcAAEgVpuya8pted/XSVbGBmi4AAAALCF0AAAAWELoAAAAsIHQBAABYQOgCAACwgNAFAABgAaELAADAAkIXAACABYQuAAAACwhdAAAAFhC6AAAALCB0AQAAWEDoAgAAsIDQBQAAYAGhCwAAwAJCFwAAgAWELgAAAAsIXQAAABYQugAAACwgdAEAAFhA6AIAALCA0AUAAGABoQsAAMACQhcAAIAFhC4AAAALCF0AAAAWELoAAAAsIHQBAABYQOiKZ/LkyVKoUCFJnz69VKtWTbZt22bjOgAAgFSO0OVh/vz50rt3bxkyZIjs2LFDypcvL/Xq1ZOYmJiUu0IAACBVIHR5GDNmjHTq1Enat28vpUqVkqlTp0rGjBllxowZKXeFAABAqkDo+q8bN27I9u3bJTo6+n9vTpo05vnmzZtT6voAAIBUIjilD8BX/PTTT3L79m3JnTu313p9vn///gRfc/36dbM4YmNjk/04AQCAfyJ0/Q7Dhw+Xt956S1KrcxMnSWoQ0bNHSh+CX+kW1S2lD8Gv8H7B3/8Op6a/kd1+498vrTTpK30ludG8+F85c+aUtGnTytmzZ73eIH0eGRmZ4Js3cOBAuXjxons5ceJEsl8wAADgnwhd/xUSEiKVKlWS1atXu9+cuLg487xGjRoJvnmhoaESFhbmtQAAACSE5kUPOl1E27ZtpXLlylK1alUZN26cXL582YxmBAAA+D0IXR6aN28u586dk8GDB8uZM2ckKipKli1bdkfnegAAgAdF6IqnR48eZgEAAEhK9OkCAACwgNAFAABgAaELAADAAkIXAACABYQuAAAACwhdAAAAFhC6AAAALCB0AQAAWEDoAgAAsIDQBQAAYAGhCwAAwAJCFwAAgAWELgAAAAsIXQAAABYQugAAACwgdAEAAFhA6AIAALCA0AUAAGABoQsAAMACQhcAAIAFhC4AAAALCF0AAAAWELoAAAAsIHQBAABYQOgCAACwgNAFAABgAaELAADAAkIXAACABYQuAAAACwhdAAAAFhC6AAAALCB0AQAAWEDoAgAAsIDQBQAAYAGhCwAAwAJCFwAAgAWELgAAAAsIXQAAABYQugAAAFJz6Dp69Kh06NBBChcuLBkyZJCiRYvKkCFD5MaNG15lgoKC7li2bNnita+FCxdKiRIlJH369FK2bFlZunSp13aXyyWDBw+WPHnymO8VHR0tBw8e9Cpz/vx5adWqlYSFhUl4eLg5tkuXLiXzuwAAAAJFioWu/fv3S1xcnHzwwQeyd+9eGTt2rEydOlVef/31O8quWrVKTp8+7V4qVark3rZp0yZp2bKlCUk7d+6URo0amWXPnj3uMiNHjpQJEyaY/W/dulUyZcok9erVk2vXrrnLaODS41i5cqUsXrxYNmzYIJ07d7bwTgAAgEAQ5NJqIB8xatQoef/99+Xw4cPumi6tCdMwFRUVleBrmjdvLpcvXzZByVG9enVTXkOWnl7evHmlT58+0rdvX7P94sWLkjt3bpk1a5a0aNFC9u3bJ6VKlZJvvvlGKleubMosW7ZMGjRoICdPnjSvvx+xsbGSNWtWs3+tMfN35yZOktQgomePlD4EAPDZv8P8jRRr92+f6tOlJ5s9e/Y71j///POSK1cueeyxx+Srr77y2rZ582bTXOhJa7F0vTpy5IicOXPGq4y+sdWqVXOX0UdtUnQCl9LyadKkMTVjAAAAv1ew+IhDhw7JxIkT5b333nOvy5w5s4wePVoeffRRE4A+++wz03S4aNEiE8SUBiqttfKkz3W9s91Zd7cyGuo8BQcHmwDolEnI9evXzeKZlAEAAKzUdA0YMCDBzu+ei/bn8vTjjz9K/fr1pWnTptKpUyf3+pw5c0rv3r1NrVSVKlVkxIgR0rp1a9MM6QuGDx9uas2cJX/+/Cl9SAAAIFBqurTvVLt27e5apkiRIu6vT506JbVr15aaNWvKtGnT7rl/DWDa2d0RGRkpZ8+e9Sqjz3W9s91Zp6MXPcs4/cS0TExMjNc+bt26ZUY0Oq9PyMCBA00o9KzpIngBAAAroSsiIsIs90NruDRw6WjEmTNnmibEe9m1a5dXeKpRo4asXr1aevXq5V6noUzXK+2Ir8FJyzghS8OR9tXq2rWrex8XLlyQ7du3u0dGrlmzxoyu1JCXmNDQULMAAAD4bJ8uDVy1atWSggULmn5c586dc29zapdmz54tISEhUqFCBfP8888/lxkzZsj06dPdZV955RV58sknTd+vhg0byrx58+Tbb79115ppc6YGsrfffluKFy9uQtgbb7xhRiRq/zBVsmRJ07ypTZs64vHmzZvSo0cPM7LxfkcuAgAA+GTo0too7TyvS758+by2ec5iMWzYMDl27Jjp2K4ToM6fP1+aNGni3q7NknPnzpVBgwaZOb40WGlH+zJlyrjL9O/f30wrofNuaY2WjoLUKSF0MlXHnDlzTNCqU6eOqXFr3LixmdsLAAAg1c3T5e+Yp8s3MQcNAH/FPF12BOQ8XQAAAKkVoQsAAIDQBQAAkDpQ0wUAAGABoQsAAMACQhcAAIAFhC4AAAALCF0AAAAWELoAAAAsIHQBAABYQOgCAACwgNAFAABgAaELAADAAkIXAACABYQuAAAACwhdAAAAFhC6AAAALCB0AQAAWEDoAgAAsIDQBQAAYAGhCwAAwAJCFwAAgAWELgAAAAsIXQAAABYQugAAACwgdAEAAFhA6AIAALCA0AUAAGABoQsAAMACQhcAAIAFhC4AAAALCF0AAAAWELoAAAAsIHQBAABYQOgCAACwgNAFAABgAaELAADAAkIXAACABYQuAACA1B66ChUqJEFBQV7LiBEjvMp899138vjjj0v69Oklf/78MnLkyDv2s3DhQilRooQpU7ZsWVm6dKnXdpfLJYMHD5Y8efJIhgwZJDo6Wg4ePOhV5vz589KqVSsJCwuT8PBw6dChg1y6dCmZzhwAAASaFK/pGjp0qJw+fdq99OzZ070tNjZW6tatKwULFpTt27fLqFGj5M0335Rp06a5y2zatElatmxpQtLOnTulUaNGZtmzZ4+7jAa1CRMmyNSpU2Xr1q2SKVMmqVevnly7ds1dRgPX3r17ZeXKlbJ48WLZsGGDdO7c2eI7AQAAUrMUD11ZsmSRyMhI96KByDFnzhy5ceOGzJgxQ0qXLi0tWrSQP//5zzJmzBh3mfHjx0v9+vWlX79+UrJkSRk2bJhUrFhRJk2a5K7lGjdunAwaNEheeOEFKVeunHz00Udy6tQpWbRokSmzb98+WbZsmUyfPl2qVasmjz32mEycOFHmzZtnygEAAPh96NLmxBw5ckiFChVMTdatW7fc2zZv3ixPPPGEhISEuNdpDdWBAwfkl19+cZfR5kJPWkbXqyNHjsiZM2e8ymTNmtWEK6eMPmqTYuXKld1ltHyaNGlMzRgAAMDvFSwpSGuttFYqe/bspplw4MCBponRqcnSsFS4cGGv1+TOndu9LVu2bObRWedZRtc75Txfl1iZXLlyeW0PDg42x+WUScj169fN4tkcCgAAYKWma8CAAXd0jo+/7N+/35Tt3bu31KpVyzT5denSRUaPHm2a9TyDjC8bPny4qTVzFu3oDwAAYKWmq0+fPtKuXbu7lilSpEiC67XJT5sXjx49Ko888ojp43X27FmvMs5z3eY8JlTGc7uzTkcvepaJiopyl4mJifHahx6Hjmh0Xp8QrZnT4OhZ00XwAgAAVkJXRESEWX6LXbt2mX5UTlNfjRo15C9/+YvcvHlT0qVLZ9bp6EINZNq06JRZvXq19OrVy70fLaPrlTZPanDSMk7I0nCkfbW6du3q3seFCxfMCMlKlSqZdWvWrJG4uDgTBBMTGhpqFgAAAJ/tSK+d13VU4b/+9S85fPiwGan46quvSuvWrd2B6o9//KPpRK/TQeh0DvPnzzejFT1rl1555RUz8lCbJrXZUqeU+Pbbb6VHjx5muzZnaiB7++235auvvpLdu3fLSy+9JHnz5jVTSygd9agjIDt16iTbtm2TjRs3mtfraEktBwAA4Lcd6bWGSKdk0JCkfbi0RkpDl2eg0n5SK1askO7du5saqJw5c5pJTj3nz6pZs6bMnTvXTAnx+uuvS/Hixc1UEGXKlHGX6d+/v1y+fNm8Tmu0dEoIDWo6mapDQ58GrTp16pjatsaNG5u5vQAAAJJCkEsnskKS0GZLDYoXL140M9v7u3MT/zPXmb+L6PmfWk8A8Dc2/g7zN1Ks3b9TfJ4uAACAQEDoAgAAsIDQBQAAYAGhCwAAwAJCFwAAgAWELgAAAAsIXQAAABYQugAAACwgdAEAAFhA6AIAALCA0AUAAGABoQsAAMACQhcAAIAFhC4AAAALCF0AAAAWELoAAAAsCLbxTeCfInr2SPbvcW7ipGT/HgAA+AJqugAAACygpgupvjYNAABfQE0XAACABYQuAAAACwhdAAAAFhC6AAAALCB0AQAAWEDoAgAAsIDQBQAAYAGhCwAAwAJCFwAAgAWELgAAAAsIXQAAABYQugAAACwgdAEAAFhA6AIAALCA0AUAAGABoQsAAMACQhcAAIAFhC4AAAALCF0AAAAWELoAAAAsIHQBAACk5tC1bt06CQoKSnD55ptvTJmjR48muH3Lli1e+1q4cKGUKFFC0qdPL2XLlpWlS5d6bXe5XDJ48GDJkyePZMiQQaKjo+XgwYNeZc6fPy+tWrWSsLAwCQ8Plw4dOsilS5csvBMAACAQpFjoqlmzppw+fdpr6dixoxQuXFgqV67sVXbVqlVe5SpVquTetmnTJmnZsqUJSTt37pRGjRqZZc+ePe4yI0eOlAkTJsjUqVNl69atkilTJqlXr55cu3bNXUYD1969e2XlypWyePFi2bBhg3Tu3NnSuwEAAFK7IJdWA/mAmzdvykMPPSQ9e/aUN954w13TpSFMw1RUVFSCr2vevLlcvnzZBCVH9erVTXkNWXp6efPmlT59+kjfvn3N9osXL0ru3Lll1qxZ0qJFC9m3b5+UKlXK1LA5gW/ZsmXSoEEDOXnypHn9/YiNjZWsWbOa/WuNGQAAv8e5iZOS/Q2M6NlDAl2spfu3z/Tp+uqrr+Tnn3+W9u3b37Ht+eefl1y5csljjz1mynnavHmzaS70pLVYul4dOXJEzpw541VG39hq1aq5y+ijNil61rBp+TRp0piascRcv37dXCjPBQAAwKdD14cffmjCUr58+dzrMmfOLKNHjzZ9tpYsWWJClzYdegYvDVRaa+VJn+t6Z7uz7m5lNNR5Cg4OluzZs7vLJGT48OEmwDlL/vz5f9d7AAAAUq8kD10DBgxItIO8s+zfv9/rNdqEt3z5ctMvy1POnDmld+/eplaqSpUqMmLECGndurWMGjVKfMHAgQNNVaSznDhxIqUPCQAA+KjgpN6h9p1q167dXcsUKVLE6/nMmTMlR44cphnxXjSAaWd3R2RkpJw9e9arjD7X9c52Z52OXvQs4/QT0zIxMTFe+7h165YZ0ei8PiGhoaFmAQAAsB66IiIizHK/tKO7hq6XXnpJ0qVLd8/yu3bt8gpPNWrUkNWrV0uvXr3c6zSU6XqlHfE1OGkZJ2Rp3yvtq9W1a1f3Pi5cuCDbt293j4xcs2aNxMXFmZAHAADgc6HrQWm40c7uOl1EfLNnz5aQkBCpUKGCef7555/LjBkzZPr06e4yr7zyijz55JOm71fDhg1l3rx58u2338q0adPMdm3O1ED29ttvS/HixU0I09GROiJR+4epkiVLSv369aVTp05mxKOOpOzRo4cZ2Xi/IxcBAAB8OnRpB3qds0snN03IsGHD5NixY6Zju5aZP3++NGnSxL1dXzt37lwZNGiQvP766yZYLVq0SMqUKeMu079/fzOthM67pTVa2iFfp4TQyVQdc+bMMUGrTp06ZtRi48aNzdxeAAAAqWqertSAeboAAEmJebrsCLh5ugAAAFIzQhcAAIAFhC4AAAALCF0AAAAWELoAAAAsIHQBAABYQOgCAACwgNAFAABgAaELAADAAkIXAACABYQuAAAACwhdAAAAFhC6AAAALCB0AQAAWEDoAgAAsIDQBQAAYAGhCwAAwAJCFwAAgAWELgAAAAsIXQAAABYQugAAACwgdAEAAFhA6AIAALCA0AUAAGABoQsAAMACQhcAAIAFhC4AAAALCF0AAAAWELoAAAAsIHQBAABYQOgCAACwgNAFAABgAaELAADAAkIXAACABYQuAAAACwhdAAAAFhC6AAAALCB0AQAAWEDoAgAA8OfQ9de//lVq1qwpGTNmlPDw8ATLHD9+XBo2bGjK5MqVS/r16ye3bt3yKrNu3TqpWLGihIaGSrFixWTWrFl37Gfy5MlSqFAhSZ8+vVSrVk22bdvmtf3atWvSvXt3yZEjh2TOnFkaN24sZ8+efeBjAQAA8LnQdePGDWnatKl07do1we23b982IUfLbdq0SWbPnm0C1eDBg91ljhw5YsrUrl1bdu3aJb169ZKOHTvK8uXL3WXmz58vvXv3liFDhsiOHTukfPnyUq9ePYmJiXGXefXVV+Uf//iHLFy4UNavXy+nTp2SP/zhDw90LAAAAL9HkMvlckky0vCiYenChQte67/++mt59tlnTQDKnTu3WTd16lR57bXX5Ny5cxISEmK+XrJkiezZs8f9uhYtWph9LVu2zDzXmq0qVarIpEmTzPO4uDjJnz+/9OzZUwYMGCAXL16UiIgImTt3rjRp0sSU2b9/v5QsWVI2b94s1atXv69juR+xsbGSNWtW8z3DwsKS6B0EAASqcxP/c29LThE9e0igi7V0/06xPl0aeMqWLesOOUprqPTE9+7d6y4THR3t9Toto+uV1kxt377dq0yaNGnMc6eMbr9586ZXmRIlSkiBAgXcZe7nWBJy/fp1U8ZzAQAA8KnQdebMGa+Qo5znuu1uZTTcXL16VX766SfTNJhQGc99aE1V/H5l8cvc61gSMnz4cJOMnUVr2AAAAH536NLmuqCgoLsu2nQXKAYOHGiqIp3lxIkTKX1IAADARwU/SOE+ffpIu3bt7lqmSJEi97WvyMjIO0YZOiMKdZvzGH+UoT7X9tYMGTJI2rRpzZJQGc99aDOk9gPzrO2KX+Zex5IQHVGpCwAAQJLWdGmHdO0Pdbflfjud16hRQ3bv3u01ynDlypUmUJUqVcpdZvXq1V6v0zK6Xun3qlSpklcZ7Uivz50yuj1dunReZQ4cOGCmiHDK3M+xAAAAWKvpehAaas6fP28etd+VTvmgdK4tnSurbt26JtC0adNGRo4cafpODRo0yMyn5dQedenSxYxK7N+/v/zpT3+SNWvWyIIFC8yIRodOF9G2bVupXLmyVK1aVcaNGyeXL1+W9u3bm+3a16pDhw6mXPbs2U2Q0pGNGrR05KK6n2MBAADwydClc1zpfFeOChUqmMe1a9dKrVq1TLPg4sWLzTxeGoAyZcpkwtPQoUPdrylcuLAJWDrP1vjx4yVfvnwyffp0M7LQ0bx5czOtg34/DUtRUVFmOgnPjvFjx441oxp1UlQdcaivnzJlinv7/RwLAACAT8/TFUiYpwsAkJSYp8uOVD9PFwAAQCAhdAEAAFhA6AIAALCA0AUAAGABoQsAAMACQhcAAIAFhC4AAAALCF0AAAAWELoAAAAsIHQBAAD482cvAgCA3yeiZw/ewlSEmi4AAAALCF0AAAAWELoAAAAsIHQBAABYQOgCAACwgNAFAABgAaELAADAAkIXAACABYQuAAAACwhdAAAAFhC6AAAALCB0AQAAWEDoAgAAsIDQBQAAYAGhCwAAwIJgG98kULhcLvMYGxub0ocCAADuk3Pfdu7jyYXQlYR+/vln85g/f/6k3C0AALB0H8+aNWuy7Z/QlYSyZ89uHo8fP56sF80X/0PQoHnixAkJCwuTQMF5c70DAT/n/JwHgosXL0qBAgXc9/HkQuhKQmnS/KeLnAauQAofDj1nzjtwcL0DC9c7sATq9U7z3/t4su0/WfcOAAAAg9AFAABgAaErCYWGhsqQIUPMYyDhvLnegYCfc37OAwE/56HJ+v4GuZJ7fCQAAACo6QIAALCB5kUAAAALCF0AAAAWELoAAAAsIHTdw+TJk6VQoUKSPn16qVatmmzbtu2u5RcuXCglSpQw5cuWLStLly712q7jFgYPHix58uSRDBkySHR0tBw8eFD8+bz/9re/yeOPPy7ZsmUzi55T/PLt2rWToKAgr6V+/friz+c9a9asO85JX5far3etWrXuOG9dGjZs6DfXe8OGDfLcc89J3rx5zbEtWrTonq9Zt26dVKxY0YzuKlasmLn+v/fvha+f9+effy5PP/20REREmIkya9SoIcuXL/cq8+abb95xrfVvoD+ft17rhH7Gz5w5k6qvd0K/t7qULl3ar6738OHDpUqVKpIlSxbJlSuXNGrUSA4cOHDP19m4fxO67mL+/PnSu3dvMw3Ejh07pHz58lKvXj2JiYlJsPymTZukZcuW0qFDB9m5c6e50Lrs2bPHXWbkyJEyYcIEmTp1qmzdulUyZcpk9nnt2jXx1/PWP1B63mvXrpXNmzebjwSqW7eu/Pjjj17l9KZ7+vRp9/Lpp5+KL3nQ81Z6I/I8p2PHjnltT43XW2/EnuesP99p06aVpk2b+s31vnz5sjlPvWnejyNHjphQWbt2bdm1a5f06tVLOnbs6BVAfsvPj6+ft960NXTpzWf79u3m/PUmrn/fPOlN2fNa//Of/xRf8qDn7dAbted56Q08NV/v8ePHe52vfrSbfixO/N9tX7/e69evl+7du8uWLVtk5cqVcvPmTXNP0vcjMdbu3zplBBJWtWpVV/fu3d3Pb9++7cqbN69r+PDhCZZv1qyZq2HDhl7rqlWr5nr55ZfN13Fxca7IyEjXqFGj3NsvXLjgCg0NdX366ad+e97x3bp1y5UlSxbX7Nmz3evatm3reuGFF1y+7EHPe+bMma6sWbMmur9Aud5jx4411/vSpUt+db0d+mfwiy++uGuZ/v37u0qXLu21rnnz5q569eol2fvoi+edkFKlSrneeust9/MhQ4a4ypcv7/IX93Pea9euNeV++eWXRMsEwvXW8kFBQa6jR4/67fVWMTEx5vzXr1/vSoyt+zc1XYm4ceOG+c9Oqw89P5NJn2ttTkJ0vWd5pSnYKa//LWv1tGcZ/ZxGrZZObJ/+cN7xXblyxfxnEf+DQ7VGTP9TfOSRR6Rr167m09x9xW8970uXLknBggVN7d4LL7wge/fudW8LlOv94YcfSosWLcx/ff5yvR/UvX63k+J99AdxcXHy66+/3vG7rU0s2oRVpEgRadWqlRw/flxSg6ioKNOUpLV9GzdudK8PlOutv9t6Tvo3zp+v98WLF83j3T7M2tb9m9CViJ9++klu374tuXPn9lqvz+O36zt0/d3KO48Psk9/OO/4XnvtNfML6fnDqU1NH330kaxevVreffddU/37zDPPmO/lr+etYWLGjBny5ZdfyieffGJuSDVr1pSTJ08GzPXWPixa/a5NbZ58/Xo/qMR+t2NjY+Xq1atJ8nvjD9577z3zj0azZs3c6/Smo/3bli1bJu+//765OWkfTw1n/kqDljYhffbZZ2bRf6q0L6M2I6pAuN6nTp2Sr7/++o7fbX+73nFxcaY7wKOPPiplypRJtJyt+3fwAx4/cFcjRoyQefPmmVoOz07lWhPi0A6K5cqVk6JFi5pyderU8ct3VTsV6+LQwFWyZEn54IMPZNiwYRII9D9hvZ5Vq1b1Wp8ar3egmzt3rrz11lvmnwzPvk0aph16nfWmrDUjCxYsMP1j/JH+Q6WL5+/2Dz/8IGPHjpWPP/5YAsHs2bMlPDzc9Gvy5G/Xu3v37uYfQ1/pd0ZNVyJy5sxpOgefPXvWa70+j4yMTPA1uv5u5Z3HB9mnP5y353/BGrpWrFhhfhnvRqul9XsdOnRI/P28HenSpZMKFSq4zym1X2/tlKoB+37+0Pra9X5Qif1u60AKHcWUFD8/vkyvs9Z46I01fhNMfHqjfvjhh/32WidG/7Fwzim1X2/tAqa1+G3atJGQkBC/vd49evSQxYsXm0Fe+fLlu2tZW/dvQlci9AetUqVKpnnEs5pSn3vWbnjS9Z7llY6ccMoXLlzYXBzPMto8oaMgEtunP5y3M6pDa3e0yrly5cr3/D7aBKd9fLQa35/P25M2N+zevdt9Tqn5ejvDq69fvy6tW7f2u+v9oO71u50UPz++Skedtm/f3jx6TguSGG1+1Fohf73WidFRq845pebrrbQ7gIao+/mHyhevt8vlMoHriy++kDVr1pi/xfdi7f59313uA9C8efPMyIRZs2a5vv/+e1fnzp1d4eHhrjNnzpjtbdq0cQ0YMMBdfuPGja7g4GDXe++959q3b58Z5ZEuXTrX7t273WVGjBhh9vHll1+6vvvuOzPCq3Dhwq6rV6+6/PW89ZxCQkJcf//7312nT592L7/++qvZro99+/Z1bd682XXkyBHXqlWrXBUrVnQVL17cde3aNZe/nreO4Fq+fLnrhx9+cG3fvt3VokULV/r06V179+5N1dfb8dhjj5kRfPH5w/XWY9y5c6dZ9M/gmDFjzNfHjh0z2/V89bwdhw8fdmXMmNHVr18/87s9efJkV9q0aV3Lli277/fRH897zpw55m+anq/n77aO2nL06dPHtW7dOnOt9W9gdHS0K2fOnGbEmL+et47IXbRokevgwYPm7/crr7ziSpMmjflZTs3X29G6dWszci8h/nC9u3btakaW63F6/txeuXLFXSal7t+ErnuYOHGiq0CBAiZU6BDhLVu2uLc9+eSTZmi8pwULFrgefvhhU16HmC9ZssRruw47feONN1y5c+c2v7B16tRxHThwwOXP512wYEHzCx1/0R9apT/odevWdUVERJgfYi3fqVMnn/rj9FvOu1evXu6yej0bNGjg2rFjR6q/3mr//v3mGq9YseKOffnD9XamBIi/OOepj3re8V8TFRVl3qMiRYqYKUMe5H30x/PWr+9WXmnwzpMnjznnhx56yDw/dOiQy5/P+91333UVLVrU/BOVPXt2V61atVxr1qxJ9ddbaaDOkCGDa9q0aQnu0x+utyRwzrp4/s6m1P076L8HCAAAgGREny4AAAALCF0AAAAWELoAAAAsIHQBAABYQOgCAAAgdAEAAKQO1HQBAABYQOgCAACwgNAFAABgAaELAADAAkIXAACABYQuAAAASX7/D2SfqM7jXcxyAAAAAElFTkSuQmCC",
+      "text/plain": [
+       "<Figure size 640x480 with 1 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "162.66083680289623\n"
+     ]
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGiCAYAAADNzj2mAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAknklEQVR4nO3dCZhU1YEv8MO+yaIoiwaQuGJciGAQMRklGDKuGKNxokYzjCYjYoBJVOKCGgQkLgwKGo0KGg3GTHAUDYooOkbcwzw0ihpccAHMKKAge73vnPeqp7ttg0h1c7r79/u+a3XdulV1OHXt+ve5Z2lQKBQKAQAgIw23dgEAACoTUACA7AgoAEB2BBQAIDsCCgCQHQEFAMiOgAIAZEdAAQCyI6AAANkRUACA2h9QHnvssXDUUUeFHXfcMTRo0CDcfffdFR6PM+dfdNFFoXPnzqFFixZhwIAB4dVXX61wzAcffBBOOumk0KZNm9CuXbswePDg8PHHH2/5vwYAqJ8BZeXKlWG//fYLkyZNqvLx8ePHh4kTJ4brr78+PPXUU6FVq1Zh4MCBYfXq1WXHxHDy4osvhlmzZoUZM2ak0HPGGWds2b8EAKgzGmzJYoGxBWX69Olh0KBB6X58qdiy8m//9m/hpz/9adq3fPny0LFjxzBlypRw4oknhpdeeinstdde4Zlnngm9e/dOx8ycOTMcfvjh4e23307PBwDqt8alfLHXX389LF68OF3WKWrbtm3o06dPmDt3bgoo8TZe1imGkyge37Bhw9Ticuyxx37qddesWZO2oo0bN6bLRO3bt08hCQDIX2zI+Oijj1JjRPzer7GAEsNJFFtMyov3i4/F2w4dOlQsROPGYbvttis7prKxY8eGSy65pJRFBQC2kkWLFoUvfelLNRdQqsvIkSPDiBEjyu7Hy0Zdu3ZN/8DY0bbkHrtyk4f87S/bbPKY7X/8oxIViLrsxv9z4yaPOXbanzZ5zPZjppSoRNTlcylyPrG1rFixInTp0iW0bt16k8eWNKB06tQp3S5ZsiSN4imK93v27Fl2zNKlSys8b/369emSTfH5lTVr1ixtlcVwUi0BpVXzTR6ypkWLTR5TLWWjzmmxzabPpdbNmmzyGOcbn+dccj6Rg8/TPaOk86B07949hYzZs2dXSEuxb0nfvn3T/Xi7bNmy8Nxzz5Ud8/DDD6d+JbGvCgDAZregxPlKXnvttQodY+fNm5f6kMTLLsOGDQujR48Ou+22WwosF154YeoMUxzp06NHj/Dtb387nH766Wko8rp168JZZ52VOtAawQMAfKGA8uyzz4ZDDz207H6xb8ipp56ahhKfc845aa6UOK9JbCk5+OCD0zDi5s3/97LJ7bffnkLJN7/5zdSL97jjjktzpwAAfKGAcsghh6RhQn/vutKll16ats8SW1vuuOMOnwAAX0j8Hor9Fzds2KAGM9KoUaM0MrcUU4DUilE8AFC0du3a8N5774VVq1aplAy1bNkyDZRp2rTpFr2OgAJArREHVMS+j/Ev9dhvMX4JmrAzn1atGB7ff//99BnFvqibmozt7xFQAKg14hdgDClxLo34lzp5iYsEN2nSJLz55pvpsyrf/3RzlXSYMQDUhC35y5za8dn4hAGA7AgoAEB29EEBoE64etYrNfZeww/bvcbeq77SggIAZEdAAQCyI6AAQA2IM7GfffbZaUmYOKN6XFz34osvTo+98cYbaT6XuLZd0bJly9K+OXPm1MvPR0ABgBoyderU0KpVq/DUU0+F8ePHp2VhZs2apf6roJMsANSQfffdN4waNSr9HGdavfbaa8Ps2bPTz1SkBQUAajCglBfXrFm6dKn6r4KAAgA1JE4DX17sYxKn7i/OvhrXsylat25dvf5cBBQA2Mp22GGHdBtXaS6aV67DbH2kDwoAZLDI3oEHHhjGjRsXunfvni77XHDBBaE+E1AAqBNq++yuN998cxg8eHDo1atX2GOPPdIon29961uhvhJQAKAGVDWfyd133132c48ePcITTzxR4fFCuT4p9Y0+KABAdgQUACA7AgoAkB0BBQDIjoACAGRHQAEAsiOgAADZEVAAgOwIKABAdswkC0Dd8MjYmnuvQ0fW3HvVU1pQAKAGnHbaaWHQoEFVToHfoEGDsGzZsrKfv/KVr4QNGzZUOK5du3ZhypQpZfd33nnnMGHChDr72QkoAJCZhQsXhltvvTXUZwIKAGRm6NChYdSoUWHNmjWhvhJQACAzw4YNC+vXrw/XXHNNqK90kgWAGjJjxoywzTbbVNhXua9J1LJly9SC8vOf/zycfvrpoW3btvXuM9KCAgA15NBDDw3z5s2rsP3617+u8tjBgweH9u3bh8svv7xefj5aUACghrRq1SrsuuuuFfa9/fbbVR7buHHjcNlll6XRP2eddVaob7SgAECmjj/++DTk+JJLLgn1jRYUAMjYuHHjwsCBA6t87J133kmXicrr1q1b2HbbbUNtJ6AAUDfU0dld+/fvn7YHH3zwU49dccUVaSvvtttuCyeffHKo7QQUAKgB5WeBLe+QQw4JhULhUz+X98ADD4TK3njjjVCX6YMCAGRHQAEAsiOgAADZEVAAgOwIKABAdgQUACA7AgoAkB0BBQDIjoACAGTHTLIA1AmT502usfc6s+eZNfZe9ZUWFACoIYsXLw4/+clPwq677hqaN28eOnbsGPr16xeuu+66sGrVKp9DOVpQAKAGLFy4MIWRdu3ahTFjxoR99tknNGvWLMyfPz/ccMMNYaeddgpHH320z+L/04ICADXgzDPPDI0bNw7PPvtsOOGEE0KPHj3Cl7/85XDMMceE++67Lxx11FHpuKuuuiqFl1atWoUuXbqk53388ccVFh2MIWfGjBlhjz32CC1btgzf/e53UwvM1KlTw8477xy23XbbcPbZZ4cNGzaUPS/uHz16dPjBD34Qttlmm9CtW7dwzz33hPfffz+VIe7bd999U/mK/ud//if80z/9UwpP8X1iuX7729/WyPkioABANYtf9A8++GAYMmRICh5VadCgwf/7Ym7YMEycODG8+OKLKXA8/PDD4ZxzzqlwbAwj8Zhp06aFmTNnhjlz5oRjjz023H///Wm77bbbwq9+9avw+9//vsLzrr766tSK8+c//zkcccQR4ZRTTkmB5eSTTw7PP/982GWXXdL94orKq1evDr169UoB6oUXXghnnHFGes7TTz8dqpuAAgDV7LXXXktf+rHFo7ztt98+tVzE7dxzz037hg0bFg499NDU4tG/f//U6vG73/2uwvPWrVuX+q189atfDd/4xjdSC8rjjz8ebrrpprDXXnuFI488Mr3GI488UuF5hx9+ePjRj34Udtttt3DRRReFFStWhAMOOCAcf/zxYffdd09leOmll8KSJUvS8bHl5Kc//Wno2bNnau0ZOnRo+Pa3v/2p8lQHfVAAYCuJLREbN24MJ510UlizZk3a99BDD4WxY8eGl19+OQWI9evXp5aM2GoSL7NE8Ta2dhTFzrYx0MSgU37f0qVLK7xfvIRT/vEoXrapvC8+r1OnTukSUewvEwPJO++8E9auXZvKWSxHddKCAgDVLI7aiZdwFixYUGF/bJWIj7Vo0SLdf+ONN1Lrx7777hv+4z/+Izz33HNh0qRJ6bEYDoqaNGlS4XXia1e1L4af8sofU7ykVNW+4vN++ctfhn//939PLSuxNWbevHlh4MCBFcpSXQQUAKhm7du3D4cddli49tprw8qVKz/zuBhIYji48sorw4EHHpguu7z77rtb7fP505/+lDrQxj4q++23XwpUr7zySo28t4ACADVg8uTJ6XJN7969w5133pn6esQWld/85jfpck6jRo1Sa0rsX3LNNdekYcmxs+v111+/1T6f2Fdl1qxZ4Yknnkjljf1Xiv1Tqps+KADUCbnP7hr7jMTRM7FPx8iRI8Pbb7+d5kGJnVpjR9Q4nDj27YjDjC+//PJ0TOwAG/ujxJE1W8MFF1yQglK8rBPLFkfxDBo0KCxfvrza37tBoTiWqBaJnYbatm2bKqhNmzalf4NHxm7ykPdfaL3JY3YYelaJCkR9n577+Nse3eQxO1x5Z4lKRF2f6r02n0+xs+jrr78eunfvnmZipXZ9Rpvz/e0SDwCQHQEFAMiOgAIA1P2AEid1ufDCC9O1pziuO3YK+sUvflE2bW4Uf44z2HXu3DkdM2DAgPDqq6+WuigAQC1V8oASex7H6XfjWO84JCneHz9+fBoyVRTvxzUE4tCpp556Kq1LEHsIx441AAAlH2Ycx0rHSV3iIkRRnHo3rnxYXFgotp5MmDAhDV2Kx0W33nprml737rvvDieeeKJPBQDquZK3oBx00EFh9uzZZTPN/fd//3dawOgf//Ef0/049Gjx4sXpsk5RHHLUp0+fMHfu3CpfM877H4cmld8AgLqr5C0o5513XgoQe+65Z5oVL/ZJueyyy9JCSFEMJ+UXJCqK94uPVRYnqbnkkktKXVQAoL60oMQVD2+//fZwxx13hOeffz5MnTo1XHHFFen2i4qz6cVJXYrbokWLSlpmAKCOt6D87Gc/S60oxb4kcRnnN998M7WCnHrqqWn55ijO5R9H8RTF+z179qzyNeNUwHEDgM/y/jXX1ljlmCm8FragrFq1KjRsWPFl46We4tLNcfhxDCmxn0pRvCQUR/P07du31MUBgCycdtppoUGDBmHcuHEV9scBInH/5vjggw/CsGHDQrdu3ULTpk3DjjvuGP75n/85vPXWW5tdrvjesQx1PqAcddRRqc/JfffdF954440wffr0tPDRscceW1YRsVJHjx4d7rnnnjB//vy0CFKs3LgAEQDUVXFtmjj9xocffviFX+ODDz4IBx54YHjooYfSdB2vvfZamDZtWro94IAD0uJ+dUHJA0qc7+S73/1uWpWxR48eaYXGuDxznKyt6JxzzglDhw5NqyLGyvz444/DzJkzLfwEQJ0WR7DGqwix20NVLr744k91d5gwYUKasqPo/PPPD++++24KKHGEbNeuXdOqxw888EBo0qRJGDJkSNmx8Xnx+eXF14/vU3w8io0IsQGh/PvUuYDSunXrVBmx38knn3wS/vrXv6bWktgEVRQr4dJLL02jduLkbLGSd99991IXBQCyErs8jBkzJv0x//bbb2/28zdu3JhaS+LI2GKfzqI4M3tsHIhBJbayfB7PPPNMur3lllvCe++9V3Y/B9biAYAaFFsrYivGqFGjNvu577//fli2bFm6QlGVuD9OiBov93weO+ywQ7pt165dCjzF+zkQUACghsV+KHH6jbgkzBdRKLe+XV0loABADYt9RuIadHGer/LiKNjK4WPdunVlP8cWjtja8VnBJu6P3Sh23XXXz/V6ORNQAGAriMON77333grLvMQAEvtnFsqFinnz5pX9HAPHCSeckCZDrTz7euz3OXny5BR8tttuu7LXi31Lyk/rEZecKS92rI2zvudGQAGArSBOZBo7u06cOLFs3yGHHJL6mYwfPz4NMpk0aVL44x//WOF5sZNt7C9y2GGHpcfi7OqPPfZYCiaxdSQ+p6h///7htttuC//1X/+VpvWIE6bGjrrlxZE7cW6yGHi2ZPhz9jPJAsDWUBtnd40jWu+8884KnVxjK8iYMWPS9BzHHXdcmq7jhhtuKDumffv24cknn0zPjdN4xGARW0zikOPf/OY3adhxUbyEFFtMjjzyyLQwb3zNyi0oV155ZRgxYkS48cYbw0477ZTmMMtBg0It7GkTm6hiRcd1edq0aVP6N3ik6vHp5b3/Qus6+T8LNW/yvMmbPOb42x7d5DE7XPm/v+Sonz7PuVTbz6c4NUX8go2zksdJz6hdn9HmfH+7xAMAZEdAAQCyI6AAANkRUACA7AgoANQ6tXB8R71RKNFnI6AAUGvEScWiVatWbe2i8BmKn03xs/qizIMCQK0RJxmLU70vXbo03W/ZsmWa2p08Wk5iOImfTfyMKk8It7kEFABqlTiLalQMKeSluDLylhJQAKhVYotJ586dQ4cOHWrNwnf1RZMmTba45aRIQAGgVopfhKX6MiQ/OskCANkRUACA7AgoAEB2BBQAIDsCCgCQHQEFAMiOgAIAZEdAAQCyI6AAANkRUACA7AgoAEB2BBQAIDsCCgCQHQEFAMiOgAIAZEdAAQCyI6AAANkRUACA7AgoAEB2BBQAIDsCCgCQHQEFAMiOgAIAZEdAAQCyI6AAANkRUACA7AgoAEB2BBQAIDsCCgCQHQEFAMiOgAIAZEdAAQCyI6AAANkRUACA7AgoAEB2BBQAIDsCCgCQHQEFAMiOgAIAZEdAAQCyI6AAANkRUACA7AgoAEB2BBQAIDsCCgCQHQEFAKgfAeWdd94JJ598cmjfvn1o0aJF2GeffcKzzz5b9nihUAgXXXRR6Ny5c3p8wIAB4dVXX62OogAAtVDJA8qHH34Y+vXrF5o0aRL++Mc/hr/85S/hyiuvDNtuu23ZMePHjw8TJ04M119/fXjqqadCq1atwsCBA8Pq1atLXRwAoBZqXOoXvPzyy0OXLl3CLbfcUrave/fuFVpPJkyYEC644IJwzDHHpH233npr6NixY7j77rvDiSeeWOoiAQD1vQXlnnvuCb179w7HH3986NChQ/jqV78abrzxxrLHX3/99bB48eJ0Waeobdu2oU+fPmHu3LlVvuaaNWvCihUrKmwAQN1V8oCycOHCcN1114XddtstPPDAA+Ff//Vfw9lnnx2mTp2aHo/hJIotJuXF+8XHKhs7dmwKMcUtttAAAHVXyQPKxo0bw/777x/GjBmTWk/OOOOMcPrpp6f+Jl/UyJEjw/Lly8u2RYsWlbTMAEAdDyhxZM5ee+1VYV+PHj3CW2+9lX7u1KlTul2yZEmFY+L94mOVNWvWLLRp06bCBgDUXSUPKHEEz4IFCyrse+WVV0K3bt3KOszGIDJ79uyyx2Ofkjiap2/fvqUuDgBQC5V8FM/w4cPDQQcdlC7xnHDCCeHpp58ON9xwQ9qiBg0ahGHDhoXRo0enfioxsFx44YVhxx13DIMGDSp1cQCAWqjkAeWAAw4I06dPT/1GLr300hRA4rDik046qeyYc845J6xcuTL1T1m2bFk4+OCDw8yZM0Pz5s1LXRwAoBYqeUCJjjzyyLR9ltiKEsNL3AAAKrMWDwCQHQEFAMiOgAIAZEdAAQCyI6AAANkRUACA7AgoAEB2BBQAIDsCCgCQHQEFAMiOgAIAZEdAAQCyI6AAANkRUACA7AgoAEB2BBQAIDsCCgCQHQEFAMiOgAIAZEdAAQCyI6AAANkRUACA7AgoAEB2BBQAIDsCCgCQHQEFAMiOgAIAZEdAAQCyI6AAANkRUACA7AgoAEB2BBQAIDsCCgCQHQEFAMiOgAIAZEdAAQCyI6AAANkRUACA7AgoAEB2BBQAIDsCCgCQHQEFAMiOgAIAZEdAAQCyI6AAANkRUACA7AgoAEB2BBQAIDsCCgCQHQEFAMiOgAIAZEdAAQCyI6AAANkRUACA7AgoAEB2BBQAIDsCCgCQHQEFAMiOgAIAZEdAAQCyI6AAANkRUACA7AgoAED9Cyjjxo0LDRo0CMOGDSvbt3r16jBkyJDQvn37sM0224TjjjsuLFmypLqLAgDUEtUaUJ555pnwq1/9Kuy7774V9g8fPjzce++94a677gqPPvpoePfdd8N3vvOd6iwKAFCLVFtA+fjjj8NJJ50UbrzxxrDtttuW7V++fHm46aabwlVXXRX69+8fevXqFW655ZbwxBNPhCeffLLK11qzZk1YsWJFhQ0AqLuqLaDESzhHHHFEGDBgQIX9zz33XFi3bl2F/XvuuWfo2rVrmDt3bpWvNXbs2NC2bduyrUuXLtVVbACgrgaUadOmheeffz4Fi8oWL14cmjZtGtq1a1dhf8eOHdNjVRk5cmRqeSluixYtqo5iAwCZaFzqF4zh4Sc/+UmYNWtWaN68eUles1mzZmkDAOqHkregxEs4S5cuDfvvv39o3Lhx2mJH2IkTJ6afY0vJ2rVrw7Jlyyo8L47i6dSpU6mLAwDUQiVvQfnmN78Z5s+fX2HfD3/4w9TP5Nxzz039R5o0aRJmz56dhhdHCxYsCG+99Vbo27dvqYsDANRCJQ8orVu3DnvvvXeFfa1atUpznhT3Dx48OIwYMSJst912oU2bNmHo0KEpnBx44IGlLg4AUAuVPKB8HldffXVo2LBhakGJQ4gHDhwYJk+evDWKAgDU14AyZ86cCvdj59lJkyalDQCgMmvxAADZEVAAgOwIKABAdgQUACA7AgoAkB0BBQDIjoACAGRHQAEAsiOgAADZEVAAgOwIKABAdgQUACA7AgoAkB0BBQDIjoACAGRHQAEAsiOgAADZEVAAgOwIKABAdgQUACA7AgoAkB0BBQDIjoACAGRHQAEAsiOgAADZEVAAgOwIKABAdgQUACA7AgoAkB0BBQDIjoACAGRHQAEAsiOgAADZEVAAgOwIKABAdgQUACA7AgoAkB0BBQDIjoACAGRHQAEAsiOgAADZEVAAgOwIKABAdgQUACA7AgoAkB0BBQDIjoACAGRHQAEAsiOgAADZEVAAgOwIKABAdgQUACA7AgoAkB0BBQDIjoACAGRHQAEAsiOgAADZEVAAgOwIKABAdgQUACA7AgoAkB0BBQDIjoACANT9gDJ27NhwwAEHhNatW4cOHTqEQYMGhQULFlQ4ZvXq1WHIkCGhffv2YZtttgnHHXdcWLJkSamLAgDUUiUPKI8++mgKH08++WSYNWtWWLduXfjWt74VVq5cWXbM8OHDw7333hvuuuuudPy7774bvvOd75S6KABALdW41C84c+bMCvenTJmSWlKee+658I1vfCMsX7483HTTTeGOO+4I/fv3T8fccsstoUePHinUHHjggZ96zTVr1qStaMWKFaUuNgBQn/qgxEASbbfdduk2BpXYqjJgwICyY/bcc8/QtWvXMHfu3M+8bNS2bduyrUuXLtVdbACgrgaUjRs3hmHDhoV+/fqFvffeO+1bvHhxaNq0aWjXrl2FYzt27Jgeq8rIkSNT0CluixYtqs5iAwB17RJPebEvygsvvBAef/zxLXqdZs2apQ0AqB+qrQXlrLPOCjNmzAiPPPJI+NKXvlS2v1OnTmHt2rVh2bJlFY6Po3jiYwAAJQ8ohUIhhZPp06eHhx9+OHTv3r3C47169QpNmjQJs2fPLtsXhyG/9dZboW/fvj4RAKD0l3jiZZ04Quc///M/01woxX4lsXNrixYt0u3gwYPDiBEjUsfZNm3ahKFDh6ZwUtUIHgCg/il5QLnuuuvS7SGHHFJhfxxKfNppp6Wfr7766tCwYcM0QVscPjxw4MAwefLkUhcFAKilGlfHJZ5Nad68eZg0aVLaAAAqsxYPAJAdAQUAyI6AAgBkR0ABALIjoAAA2RFQAIDsCCgAQHYEFAAgOwIKAJAdAQUAyI6AAgBkR0ABALIjoAAA2RFQAIDsCCgAQHYEFAAgOwIKAJAdAQUAyI6AAgBkR0ABALIjoAAA2RFQAIDsCCgAQHYEFAAgOwIKAJAdAQUAyI6AAgBkR0ABALIjoAAA2RFQAIDsCCgAQHYEFAAgOwIKAJAdAQUAyI6AAgBkR0ABALIjoAAA2RFQAIDsCCgAQHYEFAAgOwIKAJAdAQUAyI6AAgBkR0ABALIjoAAA2RFQAIDsCCgAQHYEFAAgOwIKAJAdAQUAyI6AAgBkR0ABALIjoAAA2RFQAIDsCCgAQHYEFAAgOwIKAJAdAQUAyI6AAgBkR0ABALIjoAAA2RFQAIDsCCgAQHa2akCZNGlS2HnnnUPz5s1Dnz59wtNPP701iwMA1PeAcuedd4YRI0aEUaNGheeffz7st99+YeDAgWHp0qVbq0gAQCYab603vuqqq8Lpp58efvjDH6b7119/fbjvvvvCzTffHM4777wKx65ZsyZtRcuXL0+3K1asqJ7CrVy9yUM++mTTVdesuspHnfLJx59s8piP1qzb5DHONz7PueR8Ymsqfm8XCoVNH1zYCtasWVNo1KhRYfr06RX2/+AHPygcffTRnzp+1KhR8V9iUwfOAeeAc8A54BwItb8OFi1atMmssFVaUP72t7+FDRs2hI4dO1bYH++//PLLnzp+5MiR6XJQ0caNG8MHH3wQ2rdvHxo0aFDydNelS5ewaNGi0KZNm5K+Nuq5pjmf1XNd4nyu/XUdW04++uijsOOOO+Z7iWdzNGvWLG3ltWvXrlrfM34gAkr1U881Qz2r57rE+Vy767pt27b5dpLdfvvtQ6NGjcKSJUsq7I/3O3XqtDWKBABkZKsElKZNm4ZevXqF2bNnV7hsE+/37dt3axQJAMjIVrvEE/uUnHrqqaF3797ha1/7WpgwYUJYuXJl2aierSVeSopDnytfUkI910bOZ/Vclzif61ddN4g9ZbfWm1977bXhl7/8ZVi8eHHo2bNnmDhxYpqwDQCo37ZqQAEAqIq1eACA7AgoAEB2BBQAIDsCCgCQnXoZUCZNmhR23nnn0Lx58zRq6Omnn/67x991111hzz33TMfvs88+4f7776+xstaXer7xxhvD17/+9bDtttumbcCAAZv8XNj8ei5v2rRpaamIQYMGqcoSn8/RsmXLwpAhQ0Lnzp3TUM3dd9/d745qqOc4RcUee+wRWrRokaZmHz58eFi9etMLvtZnjz32WDjqqKPSdPPxd8Ddd9+9yefMmTMn7L///ulc3nXXXcOUKVOqv6CFembatGmFpk2bFm6++ebCiy++WDj99NML7dq1KyxZsqTK4//0pz+lhQ3Hjx9f+Mtf/lK44IILCk2aNCnMnz+/xstel+v5+9//fmHSpEmFP//5z4WXXnqpcNpppxXatm1bePvtt2u87HW5notef/31wk477VT4+te/XjjmmGNqrLz1pZ7jgqi9e/cuHH744YXHH3881fecOXMK8+bNq/Gy1+V6vv322wvNmjVLt7GOH3jggULnzp0Lw4cPr/Gy1yb3339/4fzzzy/84Q9/SAv3VV64t7KFCxcWWrZsWRgxYkT6HrzmmmvS9+LMmTOrtZz1LqB87WtfKwwZMqTs/oYNGwo77rhjYezYsVUef8IJJxSOOOKICvv69OlT+NGPflTtZa1P9VzZ+vXrC61bty5MnTq1GktZP+s51u1BBx1U+PWvf1049dRTBZRqqOfrrruu8OUvf7mwdu3azftA67nNred4bP/+/Svsi1+i/fr1q/ay1hXhcwSUc845p/CVr3ylwr7vfe97hYEDB1Zr2erVJZ61a9eG5557Ll0+KGrYsGG6P3fu3CqfE/eXPz4aOHDgZx7PF6vnylatWhXWrVsXtttuO1VawvM5uvTSS0OHDh3C4MGD1W011fM999yTlu2Il3jiKu177713GDNmTFrFndLV80EHHZSeU7wMtHDhwnQZ7fDDD1fNJbS1vgdrxWrGpfK3v/0t/YKIvzDKi/dffvnlKp8TZ7mt6vi4n9LVc2Xnnntuuj5a+X8KtqyeH3/88XDTTTeFefPmqcpqrOf4Rfnwww+Hk046KX1hvvbaa+HMM89MoTtOH05p6vn73/9+et7BBx8crwaE9evXhx//+Mfh5z//uSouoc/6HlyxYkX45JNPUv+f6lCvWlCoHcaNG5c6cE6fPj11lKM0Pvroo3DKKaekDslxRXGqT1z8NLZS3XDDDWlh1O9973vh/PPPD9dff71qL6HYcTO2TE2ePDk8//zz4Q9/+EO47777wi9+8Qv1XAfUqxaU+Eu5UaNGYcmSJRX2x/udOnWq8jlx/+Yczxer56IrrrgiBZSHHnoo7LvvvqqzhOfzX//61/DGG2+k3vvlv0ijxo0bhwULFoRddtlFnW9hPUdx5E6TJk3S84p69OiR/hKNlzLiiu5seT1feOGFKXT/y7/8S7ofR1nGRWfPOOOMFAjjJSK23Gd9D7Zp06baWk+ievXpxV8K8a+Z2bNnV/gFHe/H68VVifvLHx/NmjXrM4/ni9VzNH78+PSXz8yZM9Mq15T2fI5D5efPn58u7xS3o48+Ohx66KHp5zhEky2v56hfv37psk4xAEavvPJKCi7CSWnO52JftcohpBgKLTNXOlvte7BQD4exxWFpU6ZMScOlzjjjjDSMbfHixenxU045pXDeeedVGGbcuHHjwhVXXJGGv44aNcow42qo53HjxqXhhb///e8L7733Xtn20Ucflf4kqMf1XJlRPNVTz2+99VYahXbWWWcVFixYUJgxY0ahQ4cOhdGjR2/hJ163bW49x9/HsZ5/+9vfpqGwDz74YGGXXXZJoy/5bPH3apzSIW4xBlx11VXp5zfffDM9Hus41nXlYcY/+9nP0vdgnBLCMONqEsdwd+3aNX0hxmFtTz75ZNlj//AP/5B+aZf3u9/9rrD77run4+NQq/vuu6+6ilZv67lbt27pf5TKW/wFROnquTIBpXrO5+iJJ55IUxLEL9w45Piyyy5LQ7wpXT2vW7eucPHFF6dQ0rx580KXLl0KZ555ZuHDDz9UzX/HI488UuXv22LdxttY15Wf07Nnz/S5xPP5lltuKVS3BvE/1dtGAwCweepVHxQAoHYQUACA7AgoAEB2BBQAIDsCCgCQHQEFAMiOgAIAZEdAAQCyI6AAANkRUACA7AgoAEDIzf8Fm66ViP5rtgMAAAAASUVORK5CYII=",
+      "text/plain": [
+       "<Figure size 640x480 with 1 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "w = weights * POT\n",
+    "fid_mask = np.sqrt(np.sum(vertex['decay']**2,axis=-1))<5\n",
+    "w[np.logical_not(fid_mask)] = 0\n",
+    "\n",
+    "# MB efficiency\n",
+    "w *= 0.29 - 0.12*four_momenta['Gamma'][:,0]\n",
+    "w *= four_momenta['Gamma'][:,0]>0.14\n",
+    "\n",
+    "bins = np.linspace(0,5,50)\n",
+    "for k,v in four_momenta.items():\n",
+    "    #if k!='nu': continue\n",
+    "    n,_,_ = plt.hist(four_momenta[k][:,0],bins=bins,alpha=0.5,weights=w,label=k)\n",
+    "plt.xlim(0,2)\n",
+    "plt.legend()\n",
+    "plt.show()\n",
+    "\n",
+    "\n",
+    "bins = np.linspace(0,1,50)\n",
+    "for k,v in four_momenta.items():\n",
+    "    if k == 'HNL':\n",
+    "        print(np.sqrt(np.sum(four_momenta[k][:,1:]**2)))\n",
+    "        CosTheta = four_momenta[k][:,3] / np.sqrt(np.sum(four_momenta[k][:,1:]**2,axis=-1))\n",
+    "    else:\n",
+    "        CosTheta = four_momenta[k][:,3]/four_momenta[k][:,0]\n",
+    "    plt.hist(CosTheta,bins=bins,alpha=0.5,weights=w,label=k)\n",
+    "plt.legend()\n",
+    "#plt.semilogy()\n",
+    "plt.ylim(0,100)\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 28,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAisAAAGdCAYAAADT1TPdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAZeUlEQVR4nO3de4xU9d348e8CsoDCCgoCZeVmqxI0WlAeLzVQiZfYWhND+gdVMfyoGKxa+VXZaL20lrVqrZFaqjaiaVWIbdXGu/FSkwrVotZiw1ZF6hYE9dHuItEF2fPknCe7DwuLXbCzfJx5vZLTdWbOzpn9dvbMmzPnO1uVZVmWAACC6rG7HwAAwKcRKwBAaGIFAAhNrAAAoYkVACA0sQIAhCZWAIDQxAoAEFqvFFhra2tau3Zt6t+/f6qqqtrdDwcA6IL882Y3bNiQhg8fnnr06FHesZKHSm1t7e5+GADALmhsbEwjRoxIZR0r+RGV3PBr61KPvn1298MBALqg9aOP09qL69tfx8s6Vtre+slDRawAwOfLf+oUDifYAgChiRUAIDSxAgCEJlYAgNDECgAQmlgBAEITKwBAaGIFAAhNrAAAoYkVACA0sQIAhCZWAIDQxAoAEJpYAQBCEysAQGhiBQAITawAAKGJFQAgNLECAFR2rKxZsyZ961vfSvvss0/q27dvOuSQQ9Kf//znUm8WACgTvUp55x988EE65phj0pQpU9IjjzySBg8enF577bU0cODAUm4WACgjJY2VH//4x6m2tjYtWrSo/brRo0eXcpMAQJkp6dtAv//979PEiRPTtGnT0pAhQ9Lhhx+ebrvtth2u39LSkpqbmzssAEBlK2msrFq1Ki1cuDB98YtfTI899lg699xz0/nnn5/uvPPOTtevr69PNTU17Ut+VAYAqGxVWZZlpbrz3r17F0dWnnvuufbr8lh54YUX0tKlSzs9spIvbfIjK3mwjFhwVerRt0+pHiYA8B/U+tHH6Z/fuSI1NTWlAQMGxD6yMmzYsDRu3LgO1x188MHprbfe6nT96urq4ofaegEAKltJYyWfCdTQ0NDhur///e9p5MiRpdwsAFBGShor3/3ud9OyZcvS/Pnz0+uvv57uvvvudOutt6Y5c+aUcrMAQBkpaawcccQR6b777kv33HNPGj9+fPrhD3+YbrzxxjR9+vRSbhYAKCMl/ZyV3Ne+9rViAQDYFf42EAAQmlgBAEITKwBAaGIFAAhNrAAAoYkVACA0sQIAhCZWAIDQxAoAEJpYAQBCEysAQGhiBQAITawAAKGJFQAgNLECAIQmVgCA0MQKABCaWAEAQhMrAEBoYgUACE2sAAChiRUAIDSxAgCEJlYAgNDECgAQmlgBAEITKwBAaGIFAAhNrAAAoYkVACA0sQIAhCZWAIDQxAoAEJpYAQBCEysAQGhiBQAITawAAKGJFQAgNLECAIQmVgCA0MQKABCaWAEAQuu2WLnmmmtSVVVVuvDCC7trkwBAGeiWWHnhhRfSLbfckg499NDu2BwAUEZKHisffvhhmj59errtttvSwIEDS705AKDMlDxW5syZk0455ZQ0derUf7tuS0tLam5u7rAAAJWtVynvfPHixenFF18s3gbqivr6+nTVVVeV8iEBAJ8zJTuy0tjYmC644IJ01113pT59+nTpe+rq6lJTU1P7kt8HAFDZSnZkZfny5emdd95JX/7yl9uv27JlS3r22WfTz372s+Itn549e3b4nurq6mIBACh5rBx//PHpr3/9a4frzj777HTQQQelSy65ZLtQAQDo1ljp379/Gj9+fIfr9txzz7TPPvtsdz0AwI74BFsAoHJnA23rmWee6c7NAQBlwJEVACA0sQIAhCZWAIDQxAoAEJpYAQBCEysAQGhiBQAITawAAKGJFQAgNLECAIQmVgCA0MQKABCaWAEAQhMrAEBoYgUACE2sAAChiRUAIDSxAgCEJlYAgNDECgAQmlgBAEITKwBAaGIFAAhNrAAAoYkVACA0sQIAhCZWAIDQxAoAEJpYAQBCEysAQGhiBQAITawAAKGJFQAgNLECAIQmVgCA0MQKABCaWAEAQhMrAEBoYgUACE2sAAChiRUAIDSxAgBUbqzU19enI444IvXv3z8NGTIknXbaaamhoaGUmwQAykxJY+UPf/hDmjNnTlq2bFl64okn0ubNm9MJJ5yQNm7cWMrNAgBlpFcp7/zRRx/tcPmOO+4ojrAsX748HXfccaXcNABQJkoaK9tqamoqvg4aNKjT21taWoqlTXNzc7c9NgCgwk+wbW1tTRdeeGE65phj0vjx43d4jktNTU37Ultb210PDwCo9FjJz11ZsWJFWrx48Q7XqaurK46+tC2NjY3d9fAAgEp+G+i8885LDz74YHr22WfTiBEjdrhedXV1sQAAdEusZFmWvvOd76T77rsvPfPMM2n06NGl3BwAUIZ6lfqtn7vvvjs98MADxWetrFu3rrg+Px+lb9++pdw0AFAmSnrOysKFC4tzTyZPnpyGDRvWvixZsqSUmwUAykjJ3wYCAPgs/G0gACA0sQIAhCZWAIDQxAoAEJpYAQBCEysAQGhiBQAITawAAKGJFQAgNLECAIQmVgCA0MQKABCaWAEAQhMrAEBoYgUACE2sAAChiRUAIDSxAgCEJlYAgNB67e4HAADRfGnWC11a7++3HdGl9Ubs/99d3vY/39qny+tWCkdWAIDQxAoAEJpYAQBCEysAQGhiBQAITawAAKGJFQAgNJ+zAgC7+PkppfjslJ4be3ZpvS17bkmVwpEVACA0sQIAhCZWAIDQxAoAEJpYAQBCEysAQGimLgNAiY27ck2X1/3blV8o6WP5PHJkBQAITawAAKGJFQAgNLECAIQmVgCA0MQKABBat8TKzTffnEaNGpX69OmTJk2alJ5//vnu2CwAUAZK/jkrS5YsSRdddFH6xS9+UYTKjTfemE488cTU0NCQhgwZUurNA8ButzOfnTLu+g+6dp//f2CqFCU/snLDDTekWbNmpbPPPjuNGzeuiJZ+/fql22+/vdSbBgDKQEljZdOmTWn58uVp6tSp/7fBHj2Ky0uXLt1u/ZaWltTc3NxhAQAqW0lj5b333ktbtmxJ++23X4fr88vr1q3bbv36+vpUU1PTvtTW1pby4QEAnwOhZgPV1dWlpqam9qWxsXF3PyQAoJxPsN13331Tz5490/r16ztcn18eOnTodutXV1cXCwBAtxxZ6d27d5owYUJ68skn269rbW0tLh911FGl3DQAUCZKPnU5n7Z81llnpYkTJ6YjjzyymLq8cePGYnYQAFSCvV7bo8vrVtKU5DCx8s1vfjO9++676fLLLy9Oqj3ssMPSo48+ut1JtwAAuyVWcuedd16xAAB8rmcDAQBsS6wAAKGJFQAgNLECAIQmVgCA0LplNhAAlKNxP3y7S+v97fvDSv5YypkjKwBAaGIFAAhNrAAAoYkVACA0sQIAhCZWAIDQTF0GgF30jxtrurTem//1yy7f5+iH/p//P7bhyAoAEJpYAQBCEysAQGhiBQAITawAAKGJFQAgNLECAITmc1YAYBdt/O9+XVrPZ6d8No6sAAChiRUAIDSxAgCEJlYAgNDECgAQmlgBAEITKwBAaGIFAAhNrAAAoYkVACA0sQIAhCZWAIDQxAoAEJpYAQBCEysAQGhiBQAITawAAKGJFQAgNLECAIQmVgCA0MQKAFCZsbJ69eo0c+bMNHr06NS3b980duzYdMUVV6RNmzaVapMAQBnqVao7XrlyZWptbU233HJLOuCAA9KKFSvSrFmz0saNG9P1119fqs0CAGWmZLFy0kknFUubMWPGpIaGhrRw4UKxAgDs/ljpTFNTUxo0aNAOb29paSmWNs3Nzd30yACAVOkn2L7++utpwYIF6ZxzztnhOvX19ammpqZ9qa2t7a6HBwCUS6zMmzcvVVVVfeqSn6+ytTVr1hRvCU2bNq04b2VH6urqiqMvbUtjY+Ou/VQAQOW+DTR37tw0Y8aMT10nPz+lzdq1a9OUKVPS0UcfnW699dZP/b7q6upiAQDY5VgZPHhwsXRFfkQlD5UJEyakRYsWpR49fKwLABDkBNs8VCZPnpxGjhxZzP559913228bOnRoqTYLAJSZksXKE088UZxUmy8jRozocFuWZaXaLABQZkr2vkx+XkseJZ0tAABd5SQSACA0sQIAhCZWAIDQxAoAEJpYAQBCEysAQGhiBQAITawAAKGJFQAgNLECAIQmVgCA0MQKABCaWAEAQhMrAEBoYgUACE2sAAChiRUAIDSxAgCEJlYAgNDECgAQmlgBAEITKwBAaGIFAAhNrAAAoYkVACA0sQIAhCZWAIDQxAoAEJpYAQBCEysAQGhiBQAITawAAKGJFQAgNLECAIQmVgCA0MQKABCaWAEAQhMrAEBoYgUACE2sAAChiRUAIDSxAgCE1i2x0tLSkg477LBUVVWVXn755e7YJABQJrolVi6++OI0fPjw7tgUAFBmSh4rjzzySHr88cfT9ddfX+pNAQBlqFcp73z9+vVp1qxZ6f7770/9+vXr0ttF+dKmubm5lA8PAKjkIytZlqUZM2ak2bNnp4kTJ3bpe+rr61NNTU37UltbW6qHBwCUa6zMmzevOFH205aVK1emBQsWpA0bNqS6urou33e+blNTU/vS2Ni4sw8PAKj0t4Hmzp1bHDH5NGPGjElPPfVUWrp0aaquru5wW36UZfr06enOO+/c7vvydbddHwCobDsdK4MHDy6Wf+emm25KV199dfvltWvXphNPPDEtWbIkTZo0aecfKQBQkUp2gu3+++/f4fJee+1VfB07dmwaMWJEqTYLAJQZn2ALAFTu1OWtjRo1qpghBACwMxxZAQBCEysAQGhiBQAITawAAKGJFQAgNLECAIQmVgCA0MQKABCaWAEAQhMrAEBoYgUACE2sAAChiRUAIDSxAgCEJlYAgNDECgAQmlgBAEITKwBAaGIFAAhNrAAAoYkVACA0sQIAhCZWAIDQxAoAEJpYAQBCEysAQGhiBQAITawAAKGJFQAgNLECAIQmVgCA0MQKABCaWAEAQhMrAEBoYgUACE2sAAChiRUAIDSxAgCEJlYAgNDECgAQmlgBAEITKwBA5cbKQw89lCZNmpT69u2bBg4cmE477bRSbg4AKEO9SnXHv/3tb9OsWbPS/Pnz01e/+tX0ySefpBUrVpRqcwBAmSpJrORhcsEFF6TrrrsuzZw5s/36cePGlWJzAEAZK8nbQC+++GJas2ZN6tGjRzr88MPTsGHD0sknn/xvj6y0tLSk5ubmDgsAUNlKEiurVq0qvl555ZXpsssuSw8++GBxzsrkyZPT+++/v8Pvq6+vTzU1Ne1LbW1tKR4eAPA5slOxMm/evFRVVfWpy8qVK1Nra2ux/qWXXppOP/30NGHChLRo0aLi9nvvvXeH919XV5eampral8bGxs/+EwIAlXPOyty5c9OMGTM+dZ0xY8akt99+e7tzVKqrq4vb3nrrrR1+b75OvgAA7FKsDB48uFj+nfxISh4dDQ0N6dhjjy2u27x5c1q9enUaOXLkzmwSAKhwJZkNNGDAgDR79ux0xRVXFOed5IGSzwzKTZs2rRSbBADKVMk+ZyWPk169eqUzzjgjffTRR8WHwz311FPFibYAALs9VvbYY490/fXXFwsAwK7yt4EAgNDECgAQmlgBAEITKwBAaGIFAAhNrAAAoYkVACA0sQIAhCZWAIDQxAoAEJpYAQBCEysAQGhiBQAITawAAKGJFQAgNLECAIQmVgCA0MQKABBarxRYlmXF19aPPt7dDwUA6KK21+221/HPqir7T91TCaxatSqNHTt2dz8MAGAXvPHGG2nMmDGprI+sDBo0qPj61ltvpZqamt39cEJpbm5OtbW1qbGxMQ0YMGB3P5xQjI2x8bzxO2Vfs3s1NTWl/fffv/11vKxjpUeP/z2lJg8VL8idy8fF2BibneV5Y2x2heeNcdnV1/HPygm2AEBoYgUACC10rFRXV6crrrii+Iqx8bzxO2V/Y18cideo7hub0LOBAABCH1kBABArAEBoYgUACE2sAAChhY6Vhx56KE2aNCn17ds3DRw4MJ122mkdbs8/2faUU05J/fr1S0OGDEnf+9730ieffJLK3ahRo1JVVVWH5ZprrumwziuvvJK+8pWvpD59+hSfdHvttdemStLS0pIOO+ywYmxefvnlDrdV6ticeuqpxSdK5j/3sGHD0hlnnJHWrl2bKn1sVq9enWbOnJlGjx5d7GvyP/GRz2LYtGlTqvSx+dGPfpSOPvroYh+79957d7pOpe6HczfffHOxP86fE/lr1fPPP58qzbPPPpu+/vWvp+HDhxf72/vvv7/D7fkcnssvv7zY5+S/X1OnTk2vvfbazm8oC+o3v/lNNnDgwGzhwoVZQ0ND9uqrr2ZLlixpv/2TTz7Jxo8fn02dOjV76aWXsocffjjbd999s7q6uqzcjRw5MvvBD36Qvf322+3Lhx9+2H57U1NTtt9++2XTp0/PVqxYkd1zzz1Z3759s1tuuSWrFOeff3528skn5zPdiudHm0oemxtuuCFbunRptnr16uyPf/xjdtRRRxVLpY/NI488ks2YMSN77LHHsjfeeCN74IEHsiFDhmRz587NKn1sLr/88uJ5c9FFF2U1NTXb3V7J++HFixdnvXv3zm6//fbi9WnWrFnZ3nvvna1fvz6rJA8//HB26aWXZr/73e+K/e19993X4fZrrrmmeO7cf//92V/+8pfs1FNPzUaPHp199NFHO7WdkLGyefPm7Atf+EL2y1/+8lMHqEePHtm6devar8vDZsCAAVlLS0tW7rHy05/+dIe3//znPy9Cb+txuOSSS7IDDzwwqwT5c+Oggw4qdiDbxkqlj83W8hflqqqqbNOmTcVlY/N/rr322mKH2qbSx2bRokWdxkol74ePPPLIbM6cOe2Xt2zZkg0fPjyrr6/PKlXaJlZaW1uzoUOHZtddd137df/617+y6urqIvh3Rsi3gV588cW0Zs2a4m8KHH744cXho5NPPjmtWLGifZ2lS5emQw45JO23337t15144onFH7F79dVXU7nL3/bZZ599ivG57rrrOhx2zcfmuOOOS7179+4wNg0NDemDDz5I5Wz9+vVp1qxZ6Ve/+lVxWHpblTw2W3v//ffTXXfdVRzi32OPPYrrjE3HP8K29R9gMzadq9T9cP4W4fLly4u3NNrkr1f55XxM+F9vvvlmWrduXYdxyv/WX/6W2c6OU8hYWbVqVfH1yiuvTJdddll68MEHi3NWJk+eXOxkc/kAbP0Lkmu7nN9Wzs4///y0ePHi9PTTT6dzzjknzZ8/P1188cXtt1fq2ORhP2PGjDR79uw0ceLETtep1LFpc8kll6Q999yzCN38XIMHHnig/bZKH5s2r7/+elqwYEHxu9XG2HSuUsflvffeS1u2bOn0Zy/nn3tntY3Ff2KcujVW5s2bt92JodsuK1euTK2trcX6l156aTr99NPThAkT0qJFi4rb77333lSOujo2uYsuuqgIt0MPPbR4Yf7JT35S7Fzzk0oreWzyMdiwYUOqq6tLlWJnnje5/OTHl156KT3++OOpZ8+e6cwzzywirxzt7Njk8iO6J510Upo2bVpxhK4c7cq4wO7Wqzs3Nnfu3OJfvp9mzJgx6e233y7+e9y4ce3X539fIL8t/9dgbujQodudeZ2/BdB22+dNV8emM/khtfxtoHxWw4EHHlj8/G1jUUlj89RTTxWHFrf9WxT5UZbp06enO++8s2LHps2+++5bLF/60pfSwQcfXMxqWbZsWTrqqKMqfmzymVFTpkwp3hq79dZbO4xhOT1vPsu+Zlvlth/uqvx3KI/9zp4T5fxz76y2scjHJT+do01+OZ+tGTZWBg8eXCz/Tn4kJX/Byc8jOPbYY4vrNm/eXLwYjxw5sric71zzaXXvvPNOMV0u98QTT6QBAwZ0iJzPi66OTWfyqbn5+6Vt45CPTX5UKh+ztvMR8rHJQyZ/O61cx+amm25KV199dYcXn/z98yVLlhRBV8lj05m2I5htR+QqeWzyIyp5qLQdxc1/n7ZWTmPzWZ4z2yq3/XBX5ee85c+VJ598sv0jNfLfp/zyeeedt7sfXhj5xwHkwZKPS1uc5Ocz/elPf0rnnnvuzt1ZFtQFF1xQzAjKpxOuXLkymzlzZjGd8P333+8wZe6EE07IXn755ezRRx/NBg8eXPZT5p577rliJlD+M+fTLH/9618XP/eZZ57Z4WzrfJrlGWecUUyzzKfY9evXr+ynWW7rzTff3G42UKWOzbJly7IFCxYUY5FPXX7yySezo48+Ohs7dmz28ccfV/TY/POf/8wOOOCA7Pjjjy/+e+uPBGhTqWPzj3/8o3jOXHXVVdlee+1V/He+bNiwoaL3w7n8OZDParnjjjuyv/3tb9m3v/3tYury1jOjKsGGDRvanxf5/jaf6p7/d/7caZu6nI9LPvvwlVdeyb7xjW+Uz9TlXD6dMv+cgzxQ+vfvX8zjz3cSW8t3uvlnaeSfd5DP7c/Xz6c9l7Ply5dnkyZNKqYR9unTJzv44IOz+fPnt7/gtMnnsx977LHFL1MeffkTptJ0FiuVOjb5TmLKlCnZoEGDip971KhR2ezZs4sX50ofm3xabv486Wyp9LE566yzOh2Xp59+uqL3w23yfwDsv//+xeet5FOZ838UVJqnn3660+dI/txpm778/e9/v4j9/Hcn/0dB/tlpO6sq/5/SHAACAPjsQk5dBgBoI1YAgNDECgAQmlgBAEITKwBAaGIFAAhNrAAAoYkVACA0sQIAhCZWAIDQxAoAEJpYAQBSZP8Dp3wmF8cXDEIAAAAASUVORK5CYII=",
+      "text/plain": [
+       "<Figure size 640x480 with 1 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAisAAAGdCAYAAADT1TPdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAb9ElEQVR4nO3dC3BU5f038CcQCKAkgnItkZutyqAvFpR6qQOV8fK3VWcc2ulQFYeXioNVq1Mlo/XSC7FirVNqKdoRnVaFaqu0eB+UOv8K1aLUYgeqRUoKgtcm6NhwyXnnnL5JCTc3ups82f18Zg7J7h72nH04e/bLc57fs2VJkiQBACBSXTp6BwAA9kdYAQCiJqwAAFETVgCAqAkrAEDUhBUAIGrCCgAQNWEFAIhaeYhYU1NT2LRpU+jdu3coKyvr6N0BAHKQzje7devWMHjw4NClS5fiDitpUKmuru7o3QAAPoa6urowZMiQUNRhJe1RSZ0U/ieUh24dvTsAQA52hO3hf8OjLZ/jRR1Wmi/9pEGlvExYAYBO4f9/62C+hnAYYAsARE1YAQCiJqwAAFETVgCAqAkrAEDUhBUAIGrCCgAQNWEFAIiasAIARE1YAQCiJqwAAFETVgCAqAkrAEDUhBUAIGrCCgAQNWEFAIiasAIARE1YAQCiJqwAAKUdVjZu3Bi+9rWvhYMPPjj07NkzHHXUUeFPf/pToTcLABSJ8kI++XvvvRdOPPHEMHHixPDYY4+Ffv36hVdffTX06dOnkJsFAIpIQcPKD37wg1BdXR0WLFjQct/w4cMLuUkAoMgU9DLQb3/72zBu3LgwefLk0L9//3DMMceEO++8c5/rNzY2hoaGhlYLAFDaChpW1q1bF+bNmxc+/elPhyeeeCJcfPHF4dJLLw333HPPXtevra0NVVVVLUvaKwMAlLayJEmSQj159+7ds56V5557ruW+NKy88MILYfny5XvtWUmXZmnPShpYJoSzQ3lZt0LtJgCQRzuS7WFZWBzq6+tDZWVl3D0rgwYNCqNGjWp135FHHhk2bNiw1/UrKiqyF7XrAgCUtoKGlbQSaO3ata3u+9vf/haGDh1ayM0CAEWkoGHlm9/8ZlixYkWYPXt2eO2118J9990X7rjjjjBz5sxCbhYAKCIFDSvHHntseOihh8L9998fRo8eHb773e+G2267LUyZMqWQmwUAikhB51lJffGLX8wWAICPw3cDAQBRE1YAgKgJKwBA1IQVACBqwgoAEDVhBQCImrACAERNWAEAoiasAABRE1YAgKgJKwBA1IQVACBqwgoAEDVhBQCImrACAERNWAEAoiasAABRE1YAgKgJKwBA1IQVACBqwgoAEDVhBQCImrACAERNWAEAoiasAABRE1YAgKgJKwBA1IQVACBqwgoAEDVhBQCImrACAERNWAEAoiasAABRE1YAgKgJKwBA1IQVACBqwgoAEDVhBQCImrACAERNWAEAoiasAABRE1YAgKi1W1i56aabQllZWbj88svba5MAQBFol7DywgsvhPnz54ejjz66PTYHABSRgoeV999/P0yZMiXceeedoU+fPoXeHABQZAoeVmbOnBnOPPPMMGnSpI9ct7GxMTQ0NLRaAIDSVl7IJ1+4cGF48cUXs8tAuaitrQ033nhjIXcJAOhkCtazUldXFy677LJw7733hh49euT0d2pqakJ9fX3Lkj4HAFDaCtazsnLlyvDmm2+Gz372sy337dy5Mzz77LPhJz/5SXbJp2vXrq3+TkVFRbYAABQ8rJxyyinhL3/5S6v7LrzwwnDEEUeEq6++eo+gAgDQrmGld+/eYfTo0a3uO+CAA8LBBx+8x/0AAPtiBlsAoHSrgXa3bNmy9twcAFAE9KwAAFETVgCAqAkrAEDUhBUAIGrCCgAQNWEFAIiasAIARE1YAQCiJqwAAFETVgCAqAkrAEDUhBUAIGrCCgAQNWEFAIiasAIARE1YAQCiJqwAAFETVgCAqAkrAEDUhBUAIGrCCgAQNWEFAIiasAIARE1YAQCiJqwAAFETVgCAqAkrAEDUhBUAIGrCCgAQNWEFAIiasAIARE1YAQCiJqwAAFETVgCAqAkrAEDUhBUAIGrCCgAQNWEFAIiasAIARE1YAQCiJqwAAFETVgCA0g0rtbW14dhjjw29e/cO/fv3D+ecc05Yu3ZtITcJABSZgoaV3//+92HmzJlhxYoV4amnngrbt28Pp556avjggw8KuVkAoIiUF/LJH3/88Va377777qyHZeXKleHkk08u5KYBgCJR0LCyu/r6+uxn37599/p4Y2NjtjRraGhot30DAEp8gG1TU1O4/PLLw4knnhhGjx69zzEuVVVVLUt1dXV77R4AUOphJR27snr16rBw4cJ9rlNTU5P1vjQvdXV17bV7AEApXwa65JJLwpIlS8Kzzz4bhgwZss/1KioqsgUAoF3CSpIk4Rvf+EZ46KGHwrJly8Lw4cMLuTkAoAiVF/rSz3333RcWL16czbWyefPm7P50PErPnj0LuWkAoEgUdMzKvHnzsrEnEyZMCIMGDWpZFi1aVMjNAgBFpOCXgQAAPgnfDQQARE1YAQCiJqwAAFETVgCAqLXrdwMRt52njM1pva5LVxZ8XwCgmZ4VACBqwgoAEDVhBQCImrACAERNWAEAoiasAABRE1YAgKiV3Dwruc4lUorziXSG19uWf798vubOMAdNvvdx3f1jclpvxFdX5bQe8XE+pLPQswIARE1YAQCiJqwAAFETVgCAqAkrAEDUhBUAIGqdonR554Qxoay8R7uXjG66+oSc1hv8g+fy+nwD/tRYNGW3HbXtfLfNlnEVOa03eGkoGrGXJOf6firEe6pYlNrr7chzV1uO11w/U0qJnhUAIGrCCgAQNWEFAIiasAIARE1YAQCiJqwAAFETVgCAqHWKeVaO+/7KUHFgt/2u8/yYrvmfV2Bc7nXxMc8Rkm/FNF9Arm2d7/lTimmOkI6aS6ctx1ZHvVfy3Tb5fh3mWencx2vXEponR88KABA1YQUAiJqwAgBETVgBAKImrAAAURNWAICodYrS5SUPnhC6VvTY7zqDQ25lYevuH5Pzdkd89bkOKSeMvcyyM5Tg5VoanOtryfe/Sa7lyIWQ62vZMq4ir68l3/8mhSihj700uJRKVYtNIaYrKCV6VgCAqAkrAEDUhBUAIGrCCgAQNWEFAIiasAIARK1dwsrtt98ehg0bFnr06BHGjx8fnn/++fbYLABQBAo+z8qiRYvCFVdcEX72s59lQeW2224Lp512Wli7dm3o379/Ts8x6Ed/DOVl3fKyPyO+uqogdfEdUTvfGeZcyHUf8z0fy+ClIfq5P3LVlnlt8jl/Sr7nJsn1+O+oeYYK8Z7K95w2HTmvTGc431C8Ct6zcuutt4bp06eHCy+8MIwaNSoLLb169Qp33XVXoTcNABSBgoaVbdu2hZUrV4ZJkyb9d4NdumS3ly9fvsf6jY2NoaGhodUCAJS2goaVt99+O+zcuTMMGDCg1f3p7c2bN++xfm1tbaiqqmpZqqurC7l7AEAnEFU1UE1NTaivr29Z6urqOnqXAIBiHmB7yCGHhK5du4YtW7a0uj+9PXDgwD3Wr6ioyBYAgHbpWenevXsYO3ZsWLr0v6UZTU1N2e3jjz++kJsGAIpEWZIkSaFLly+44IIwf/78cNxxx2Wly7/61a/CmjVr9hjLsrt0gG06dmVCOPsjS5fzXfrakXItf821FDTfr3nd/WNyXrfHql5RfyV6ZzgecpXv90C+S4g7Q1vnu9Q432Xn+f43aUvJdL5fS6l5YtOfc173tMH/J3R2O5LtYVlYnA3pqKysjH+ela985SvhrbfeCtddd102qHbMmDHh8ccf/8igAgDQLmEldckll2QLAECnrgYCANidsAIARE1YAQCiJqwAAFETVgCA0p5n5ZNonmfl8xOuD+XlPUrm6+I7w3wUuYp9ro7Y57TpyH3Mdbu5zr+R72OhM7yX//F/d+a03tCfdy25cwPFbUee51nRswIARE1YAQCiJqwAAFETVgCAqAkrAEDUhBUAIGrt8kWG7aEQJX2xlxrHvn+doW3y/bX3+S73LURJcr7bJt/HYSGO63w/Z67r9RiX2/HQdWnHlH8rhY5vaoGUf5c96VkBAKImrAAAURNWAICoCSsAQNSEFQAgasIKABA1YQUAiFpZkiRJiFRDQ0OoqqoKE8LZobysW4hV7POdxL5/HTl/REe1Tb5fb2eYm6EzHIf5njMj38fhlnEVed2/XJ+vEHMSlZq2vOdjfy/nYkeyPSwLi0N9fX2orKwMn5SeFQAgasIKABA1YQUAiJqwAgBETVgBAKImrAAAUVO6DCWmM5QQA53bDqXLAEApcRkIAIiasAIARE1YAQCiJqwAAFETVgCAqAkrAEDUyjt6B4D2Zf4UoLPRswIARE1YAQCiJqwAAFETVgCAqAkrAEDUhBUAoDTDyvr168O0adPC8OHDQ8+ePcPIkSPD9ddfH7Zt21aoTQIARahg86ysWbMmNDU1hfnz54fDDjssrF69OkyfPj188MEH4ZZbbinUZgGAIlOWJEnSXhubM2dOmDdvXli3bl1O6zc0NISqqqowIZwdysu6FXz/AIBPbkeyPSwLi0N9fX2orKzsXDPYpjvdt2/ffT7e2NiYLbuGFQCgtLXbANvXXnstzJ07N1x00UX7XKe2tjbrSWleqqur22v3AIBiCSuzZs0KZWVl+13S8Sq72rhxYzj99NPD5MmTs3Er+1JTU5P1vjQvdXV1H+9VAQClO2blrbfeCu+8885+1xkxYkTo3r179vumTZvChAkTwuc+97lw9913hy5dcs9HxqwAQOfT4WNW+vXrly25SHtUJk6cGMaOHRsWLFjQpqACAFDQAbZpUEl7VIYOHZqVKqc9Ms0GDhyo9QGAjg0rTz31VDaoNl2GDBnS6rF2rJYGADq5gl2XmTp1ahZK9rYAAOTKIBIAIGrCCgAQNWEFAIiasAIARE1YAQCiJqwAAFETVgCAqAkrAEDUhBUAIGrCCgAQNWEFAIiasAIARE1YAQCiJqwAAFETVgCAqAkrAEDUhBUAIGrCCgAQNWEFAIiasAIARE1YAQCiJqwAAFETVgCAqAkrAEDUhBUAIGrCCgAQNWEFAIiasAIARE1YAQCiJqwAAFETVgCAqAkrAEDUhBUAIGrCCgAQNWEFAIiasAIARE1YAQCiJqwAAFETVgCAqAkrAEDUhBUAIGrCCgAQtXYJK42NjWHMmDGhrKwsrFq1qj02CQAUiXYJK1dddVUYPHhwe2wKACgyBQ8rjz32WHjyySfDLbfcUuhNAQBFqLyQT75ly5Ywffr08PDDD4devXrldLkoXZo1NDQUcvcAgFLuWUmSJEydOjXMmDEjjBs3Lqe/U1tbG6qqqlqW6urqQu0eAFCsYWXWrFnZQNn9LWvWrAlz584NW7duDTU1NTk/d7pufX19y1JXV9fW3QMAikxZknaBtMFbb70V3nnnnf2uM2LEiPDlL385/O53v8vCS7OdO3eGrl27hilTpoR77rnnI7eVXgZKe1gmhLNDeVm3tuwmANBBdiTbw7KwOOt4qKysbP+wkqsNGza0GnOyadOmcNppp4UHH3wwjB8/PgwZMuQjn0NYAYDOJ99hpWADbA899NBWtw888MDs58iRI3MKKgAAKTPYAgClW7q8q2HDhmUVQgAAbaFnBQCImrACAERNWAEAoiasAABRE1YAgKgJKwBA1IQVACBqwgoAEDVhBQCImrACAERNWAEAoiasAABRE1YAgKgJKwBA1IQVACBqwgoAEDVhBQCImrACAERNWAEAoiasAABRE1YAgKgJKwBA1IQVACBqwgoAEDVhBQCImrACAERNWAEAoiasAABRE1YAgKgJKwBA1IQVACBqwgoAEDVhBQCImrACAERNWAEAoiasAABRE1YAgKgJKwBA1IQVACBqwgoAEDVhBQCImrACAJRuWHnkkUfC+PHjQ8+ePUOfPn3COeecU8jNAQBFqLxQT/zrX/86TJ8+PcyePTt84QtfCDt27AirV68u1OYAgCJVkLCSBpPLLrsszJkzJ0ybNq3l/lGjRhVicwBAESvIZaAXX3wxbNy4MXTp0iUcc8wxYdCgQeGMM874yJ6VxsbG0NDQ0GoBAEpbQcLKunXrsp833HBDuPbaa8OSJUuyMSsTJkwI77777j7/Xm1tbaiqqmpZqqurC7F7AEAn0qawMmvWrFBWVrbfZc2aNaGpqSlb/5prrgnnnntuGDt2bFiwYEH2+AMPPLDP56+pqQn19fUtS11d3Sd/hQBA6YxZufLKK8PUqVP3u86IESPCG2+8sccYlYqKiuyxDRs27PPvpuukCwDAxwor/fr1y5aPkvakpKFj7dq14aSTTsru2759e1i/fn0YOnRoWzYJAJS4glQDVVZWhhkzZoTrr78+G3eSBpS0Mig1efLkQmwSAChSBZtnJQ0n5eXl4bzzzgsffvhhNjnc008/nQ20BQDIVVmSJEmIVFq6nFYFTQhnh/Kybh29OwBADnYk28OysDgrlkmvtnxSvhsIAIiasAIARE1YAQCiJqwAAFETVgCAqAkrAEDUhBUAIGrCCgAQNWEFAIiasAIARE1YAQCiJqwAAFETVgCAqAkrAEDUhBUAIGrCCgAQNWEFAIiasAIARK08RCxJkuznjrA9hP/8CgBELvvc3uVzvKjDyjvvvJP9/N/waEfvCgDwMT7Hq6qqQlGHlb59+2Y/N2zYkJcXW0waGhpCdXV1qKurC5WVlR29O1HRNtrGceM95VzTserr68Ohhx7a8jle1GGlS5f/DKlJg4oP5L1L20XbaJu2ctxom4/DcaNdPu7n+CdlgC0AEDVhBQCIWtRhpaKiIlx//fXZT7SN48Z7yvnGuTgmPqPar23KknzVFQEAlFrPCgCAsAIARE1YAQCiJqwAAFGLOqw88sgjYfz48aFnz56hT58+4Zxzzmn1eDqz7Zlnnhl69eoV+vfvH771rW+FHTt2hGI3bNiwUFZW1mq56aabWq3z8ssvh89//vOhR48e2Uy3N998cygljY2NYcyYMVnbrFq1qtVjpdo2Z511VjajZPq6Bw0aFM4777ywadOmUOpts379+jBt2rQwfPjw7FwzcuTIrIph27ZtodTb5vvf/3444YQTsnPsQQcdtNd1SvU8nLr99tuz83F6TKSfVc8//3woNc8++2z40pe+FAYPHpydbx9++OFWj6c1PNddd112zknfX5MmTQqvvvpq2zeUROrBBx9M+vTpk8ybNy9Zu3Zt8sorrySLFi1qeXzHjh3J6NGjk0mTJiUvvfRS8uijjyaHHHJIUlNTkxS7oUOHJt/5zneSN954o2V5//33Wx6vr69PBgwYkEyZMiVZvXp1cv/99yc9e/ZM5s+fn5SKSy+9NDnjjDPSSrfs+GhWym1z6623JsuXL0/Wr1+f/OEPf0iOP/74bCn1tnnssceSqVOnJk888UTy97//PVm8eHHSv3//5Morr0xKvW2uu+667Li54oorkqqqqj0eL+Xz8MKFC5Pu3bsnd911V/b5NH369OSggw5KtmzZkpSSRx99NLnmmmuS3/zmN9n59qGHHmr1+E033ZQdOw8//HDy5z//OTnrrLOS4cOHJx9++GGbthNlWNm+fXvyqU99Kvn5z3++3wbq0qVLsnnz5pb70mBTWVmZNDY2JsUeVn70ox/t8/Gf/vSnWdDbtR2uvvrq5PDDD09KQXpsHHHEEdkJZPewUupts6v0Q7msrCzZtm1bdlvb/NfNN9+cnVCblXrbLFiwYK9hpZTPw8cdd1wyc+bMlts7d+5MBg8enNTW1ialKuwWVpqampKBAwcmc+bMabnvX//6V1JRUZEF/raI8jLQiy++GDZu3Jh9p8AxxxyTdR+dccYZYfXq1S3rLF++PBx11FFhwIABLfeddtpp2ZfYvfLKK6HYpZd9Dj744Kx95syZ06rbNW2bk08+OXTv3r1V26xduza89957oZht2bIlTJ8+PfziF7/IuqV3V8pts6t333033HvvvVkXf7du3bL7tE3rL2Hb9QvYtM3elep5OL1EuHLlyuySRrP08yq9nbYJ//H666+HzZs3t2qn9Lv+0ktmbW2nKMPKunXrsp833HBDuPbaa8OSJUuyMSsTJkzITrKptAF2fYOkmm+njxWzSy+9NCxcuDA888wz4aKLLgqzZ88OV111Vcvjpdo2abCfOnVqmDFjRhg3btxe1ynVtml29dVXhwMOOCALuulYg8WLF7c8Vupt0+y1114Lc+fOzd5bzbTN3pVqu7z99tth586de33txfy626q5LfLRTu0aVmbNmrXHwNDdlzVr1oSmpqZs/WuuuSace+65YezYsWHBggXZ4w888EAoRrm2TeqKK67IgtvRRx+dfTD/8Ic/zE6u6aDSUm6btA22bt0aampqQqloy3GTSgc/vvTSS+HJJ58MXbt2Deeff34W8opRW9smlfbonn766WHy5MlZD10x+jjtAh2tvD03duWVV2b/892fESNGhDfeeCP7fdSoUS33p98vkD6W/m8wNXDgwD1GXqeXAJof62xybZu9SbvU0stAaVXD4Ycfnr3+5rYopbZ5+umns67F3b+LIu1lmTJlSrjnnntKtm2aHXLIIdnymc98Jhx55JFZVcuKFSvC8ccfX/Jtk1ZGTZw4Mbs0dscdd7Rqw2I6bj7JuWZ3xXYezlX6HkrD/t6OiWJ+3W3V3BZpu6TDOZqlt9NqzWjDSr9+/bLlo6Q9KekHTjqO4KSTTsru2759e/ZhPHTo0Ox2enJNy+refPPNrFwu9dRTT4XKyspWIaezyLVt9iYtzU2vlza3Q9o2aa9U2mbN4xHStkmDTHo5rVjb5sc//nH43ve+1+rDJ71+vmjRoizQlXLb7E1zD2Zzj1wpt03ao5IGleZe3PT9tKtiaptPcszsrtjOw7lKx7ylx8rSpUtbptRI30/p7UsuuaSjdy8a6XQAaWBJ26U5nKTjmf74xz+Giy++uG1PlkTqsssuyyqC0nLCNWvWJNOmTcvKCd99991WJXOnnnpqsmrVquTxxx9P+vXrV/Qlc88991xWCZS+5rTM8pe//GX2us8///xWo63TMsvzzjsvK7NMS+x69epV9GWWu3v99df3qAYq1bZZsWJFMnfu3Kwt0tLlpUuXJieccEIycuTI5N///ndJt80///nP5LDDDktOOeWU7PddpwRoVqpt849//CM7Zm688cbkwAMPzH5Pl61bt5b0eTiVHgNpVcvdd9+d/PWvf02+/vWvZ6XLu1ZGlYKtW7e2HBfp+TYtdU9/T4+d5tLltF3S6sOXX345Ofvss4undDmVllOm8xykAaV3795ZHX96kthVetJN59JI5ztIa/vT9dOy52K2cuXKZPz48VkZYY8ePZIjjzwymT17dssHTrO0nv2kk07K3kxp6EsPmFKzt7BSqm2TniQmTpyY9O3bN3vdw4YNS2bMmJF9OJd626Rluelxsrel1Nvmggsu2Gu7PPPMMyV9Hm6W/gfg0EMPzeZbSUuZ0/8UlJpnnnlmr8dIeuw0ly9/+9vfzsJ++t5J/1OQzp3WVmXpH4XpAAIA+OSiLF0GAGgmrAAAURNWAICoCSsAQNSEFQAgasIKABA1YQUAiJqwAgBETVgBAKImrAAAURNWAICoCSsAQIjZ/wN70cEMuktjCgAAAABJRU5ErkJggg==",
+      "text/plain": [
+       "<Figure size 640x480 with 1 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "Nbins = 50\n",
+    "zmin = -60\n",
+    "zmax = 10\n",
+    "xbound = 7\n",
+    "\n",
+    "zbins = np.linspace(zmin,zmax,Nbins)\n",
+    "xbins = np.linspace(-xbound,xbound,Nbins)\n",
+    "\n",
+    "plt.hist2d(vertex['upscattering'][:,2],vertex['upscattering'][:,0],bins=(zbins,xbins),weights=w)\n",
+    "plt.show()\n",
+    "\n",
+    "plt.hist2d(vertex['decay'][:,2],vertex['decay'][:,0],bins=(zbins,xbins))#,weights=w)\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "lienv",
+   "language": "python",
+   "name": "lienv"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/resources/earthparams/densities/CCM.dat b/resources/earthparams/densities/CCM.dat
index f1576cde..e79103fc 100644
--- a/resources/earthparams/densities/CCM.dat
+++ b/resources/earthparams/densities/CCM.dat
@@ -1,50 +1,3 @@
-#
-# crust model file
-# PREM + mmc crust setup
-# crust data is consistent with MMC mediadef.
-#
-#   format as a context free grammar:
-#   [line] -> [entry] [comment]
-#   [comment] ->
-#   [comment] -> #(string comment)
-#   [entry] ->
-#   [entry] -> [object]
-#   [entry] -> [detector]
-#   [detector] -> detector (double x) (double y) (double z)
-#   [object] -> object [shape] [label] [material name] [density distribution]
-#   [shape] -> [sphere]
-#   [shape] -> [box]
-#   [sphere] -> sphere [shape coords] (double radius)
-#   [box] -> box [shape coords] (double dx) (double dy) (double dz)
-#   [shape coords] -> (double x) (double y) (double z)
-#   [label] -> (string label)
-#   [material name] -> (string material)
-#   [density distribution] -> [constant]
-#   [density distribution] -> [radial polynomial]
-#   [constant] -> constant (double density)
-#   [radial polynomial] -> radial_polynomial [axis coords] [polynomial]
-#   [axis coords] -> (double x) (double y) (double z)
-#   [polynomial] -> (int n_params) [polynomial {n_params}]
-#   [polynomial {n_params}] -> (double param_n) [polynomial {n_params-1}]
-#   [polynomial {1}] -> (double param_1)
-#
-#   format notes:
-#    - no space before first column
-#    - delimiter must be whitespaces
-#    - for overlapping objects the object defined later in the file is given precidence
-#    - MediumType must be consistent with EarthModelService::MediumType
-#    - number_of_parameters is the number of params for polynominal density function
-#    - parameters : ex.
-#       for 3 parameters p0, p1, p2 and distance from Earth's center x :
-#       density = p0 + p1*x + p2*x*x
-#
-#       CAUTION : the unit of density must be [g/cm3] (CGS unit)
-#
-#    - lines start from '#' or whitespace are treated
-#      as comments
-#
-
-
 # Detector Hall
 object box    0 0 0           0 0 0      100 100 100    surr_air           AIR         constant    0.001225   # 0.673atm x 1.205e-3(g/cm3 for 1atm)
 # adjust floor z pos if necessary: this is a guess
diff --git a/resources/earthparams/densities/MINERvA.dat b/resources/earthparams/densities/MINERvA.dat
index 944f8951..cc515459 100644
--- a/resources/earthparams/densities/MINERvA.dat
+++ b/resources/earthparams/densities/MINERvA.dat
@@ -1,49 +1,3 @@
-#
-# crust model file
-# PREM + mmc crust setup
-# crust data is consistent with MMC mediadef.
-#
-#   format as a context free grammar:
-#   [line] -> [entry] [comment]
-#   [comment] ->
-#   [comment] -> #(string comment)
-#   [entry] ->
-#   [entry] -> [object]
-#   [entry] -> [detector]
-#   [detector] -> detector (double x) (double y) (double z)
-#   [object] -> object [shape] [label] [material name] [density distribution]
-#   [shape] -> [sphere]
-#   [shape] -> [box]
-#   [sphere] -> sphere [shape coords] (double radius)
-#   [box] -> box [shape coords] (double dx) (double dy) (double dz)
-#   [shape coords] -> (double x) (double y) (double z)
-#   [label] -> (string label)
-#   [material name] -> (string material)
-#   [density distribution] -> [constant]
-#   [density distribution] -> [radial polynomial]
-#   [constant] -> constant (double density)
-#   [radial polynomial] -> radial_polynomial [axis coords] [polynomial]
-#   [axis coords] -> (double x) (double y) (double z)
-#   [polynomial] -> (int n_params) [polynomial {n_params}]
-#   [polynomial {n_params}] -> (double param_n) [polynomial {n_params-1}]
-#   [polynomial {1}] -> (double param_1)
-#
-#   format notes:
-#    - no space before first column
-#    - delimiter must be whitespaces
-#    - for overlapping objects the object defined later in the file is given precidence
-#    - MediumType must be consistent with EarthModelService::MediumType
-#    - number_of_parameters is the number of params for polynominal density function
-#    - parameters : ex.
-#       for 3 parameters p0, p1, p2 and distance from Earth's center x :
-#       density = p0 + p1*x + p2*x*x
-#
-#       CAUTION : the unit of density must be [g/cm3] (CGS unit)
-#
-#    - lines start from '#' or whitespace are treated
-#      as comments
-#
-
 # Dirt
 object box       0 0 -100            0 0 0      100 100 300    surr_earth           ROCK         constant    2.900 
 
@@ -54,12 +8,12 @@ object sphere    0 0 0            0 0 0      10    surr_air           AIR
 object box       0 0 -1           0 0 0      2 2 0.075     veto_wall    STEEL     constant   7.83   
 
 # HCAL and ECAL
-#object extr			 0 0 2.0672       0 0 0      6 0.000 0.98150 0.85000 0.49075 0.85000 -0.49075 0.000 -0.98150 -0.85000 -0.49075 -0.85000 0.49075		2 -2.0672 0 0 1 2.0672 0 0 1   HCal  PLANE    constant      1.043
-#object extr			 0 0 2.0672       0 0 0      6 0.000 1.23553 1.07000 0.61776 1.07000 -0.61776 0.000 -1.23553 -1.07000 -0.61776 -1.07000 0.61776		2 -2.0672 0 0 1 2.0672 0 0 1   Ecal  ECAL    constant      2.244
-object extr			 0 0 2.0672       0 0 0      6 0.000 1.23553 1.07000 0.61776 1.07000 -0.61776 0.000 -1.23553 -1.07000 -0.61776 -1.07000 0.61776		2 -2.07 0 0 1 2.07 0 0 1   Ecal  ECAL    constant      2.244
+#object extr	  0 0 2.0672       0 0 0      6 0.000 0.98150 0.85000 0.49075 0.85000 -0.49075 0.000 -0.98150 -0.85000 -0.49075 -0.85000 0.49075		2 -2.0672 0 0 1 2.0672 0 0 1   HCal  PLANE    constant      1.043
+#object extr	  0 0 2.0672       0 0 0      6 0.000 1.23553 1.07000 0.61776 1.07000 -0.61776 0.000 -1.23553 -1.07000 -0.61776 -1.07000 0.61776		2 -2.0672 0 0 1 2.0672 0 0 1   Ecal  ECAL    constant      2.244
+object extr		  0 0 2.0672       0 0 0      6 0.000 1.23553 1.07000 0.61776 1.07000 -0.61776 0.000 -1.23553 -1.07000 -0.61776 -1.07000 0.61776		2 -2.07 0 0 1 2.07 0 0 1   Ecal  ECAL    constant      2.244
 
 # OLD 85 cm apothem
-#object extr			 0 0 2.0672       0 0 0      6 0.000 0.98150 0.85000 0.49075 0.85000 -0.49075 0.000 -0.98150 -0.85000 -0.49075 -0.85000 0.49075		2 -2.0672 0 0 1 2.0672 0 0 1   tracker  PLANE    constant      1.043
+#object extr	  0 0 2.0672       0 0 0      6 0.000 0.98150 0.85000 0.49075 0.85000 -0.49075 0.000 -0.98150 -0.85000 -0.49075 -0.85000 0.49075		2 -2.0672 0 0 1 2.0672 0 0 1   tracker  PLANE    constant      1.043
 #object extr      0 0 0.136        0 0 0      5 -0.2475 0.8386 0.6025 -0.6336 0.000 -0.98150 -0.85000 -0.49075 -0.85000 0.49075    2 0 0 0 1 0.02567 0 0 1    Fe1   STEEL   constant   7.83 
 #object extr      0 0 0.136        0 0 0      5 0.6025 -0.6336 -0.2475 0.8386 0.000 0.98150 0.85000 0.49075 0.85000 -0.49075   2 0 0 0 1 0.02578 0 0 1    Pb1   LEAD   constant   11.29
 #object extr      0 0 0.313        0 0 0      5 -0.6025 -0.6336 0.2475 0.8386 0.85000 0.49075 0.85000 -0.49075 0.000 -0.98150   2 0 0 0 1 0.02563 0 0 1    Fe2   STEEL   constant   7.83 
@@ -68,7 +22,7 @@ object extr			 0 0 2.0672       0 0 0      6 0.000 1.23553 1.07000 0.61776 1.070
 #object extr      0 0 0.534        0 0 0      4 0.000 0.000 0.85000 0.49075 0.85000 -0.49075 0.000 -0.98150   2 0 0 0 1 0.02573 0 0 1    Fe3   STEEL   constant   7.83 
 #object extr      0 0 0.534        0 0 0      3 0.000 0.000 0.000 -0.98150 -0.85000 -0.49075    2 0 0 0 1 0.02563 0 0 1    Pb3   LEAD   constant   11.29 
 #object extr      0 0 0.895        0 0 0      6 0.000 0.98150 0.85000 0.49075 0.85000 -0.49075 0.000 -0.98150 -0.85000 -0.49075 -0.85000 0.49075   2 0 0 0 1 0.1806 0 0 1   water_target  WATER    constant      1.0
-#object extr			 0 0 1.256        0 0 0      6 0.000 0.98150 0.85000 0.49075 0.85000 -0.49075 0.000 -0.98150 -0.85000 -0.49075 -0.85000 0.49075		2 0 0 0 1 0.00795 0 0 1   Pb4  LEAD    constant      11.29
+#object extr	  0 0 1.256        0 0 0      6 0.000 0.98150 0.85000 0.49075 0.85000 -0.49075 0.000 -0.98150 -0.85000 -0.49075 -0.85000 0.49075		2 0 0 0 1 0.00795 0 0 1   Pb4  LEAD    constant      11.29
 #object extr      0 0 1.389        0 0 0      5 -0.2475 0.8386 0.6025 -0.6336 0.000 -0.98150 -0.85000 -0.49075 -0.85000 0.49075    2 0 0 0 1 0.01289 0 0 1    Fe5   STEEL   constant   7.83 
 #object extr      0 0 1.389        0 0 0      5 0.6025 -0.6336 -0.2475 0.8386 0.000 0.98150 0.85000 0.49075 0.85000 -0.49075   2 0 0 0 1 0.01317 0 0 1    Pb5   LEAD   constant   11.29 
 
diff --git a/resources/earthparams/densities/MiniBooNE.dat b/resources/earthparams/densities/MiniBooNE.dat
index 4d2add8a..2d196dd9 100644
--- a/resources/earthparams/densities/MiniBooNE.dat
+++ b/resources/earthparams/densities/MiniBooNE.dat
@@ -1,49 +1,3 @@
-#
-# crust model file
-# PREM + mmc crust setup
-# crust data is consistent with MMC mediadef.
-#
-#   format as a context free grammar:
-#   [line] -> [entry] [comment]
-#   [comment] ->
-#   [comment] -> #(string comment)
-#   [entry] ->
-#   [entry] -> [object]
-#   [entry] -> [detector]
-#   [detector] -> detector (double x) (double y) (double z)
-#   [object] -> object [shape] [label] [material name] [density distribution]
-#   [shape] -> [sphere]
-#   [shape] -> [box]
-#   [sphere] -> sphere [shape coords] (double radius)
-#   [box] -> box [shape coords] (double dx) (double dy) (double dz)
-#   [shape coords] -> (double x) (double y) (double z)
-#   [label] -> (string label)
-#   [material name] -> (string material)
-#   [density distribution] -> [constant]
-#   [density distribution] -> [radial polynomial]
-#   [constant] -> constant (double density)
-#   [radial polynomial] -> radial_polynomial [axis coords] [polynomial]
-#   [axis coords] -> (double x) (double y) (double z)
-#   [polynomial] -> (int n_params) [polynomial {n_params}]
-#   [polynomial {n_params}] -> (double param_n) [polynomial {n_params-1}]
-#   [polynomial {1}] -> (double param_1)
-#
-#   format notes:
-#    - no space before first column
-#    - delimiter must be whitespaces
-#    - for overlapping objects the object defined later in the file is given precidence
-#    - MediumType must be consistent with EarthModelService::MediumType
-#    - number_of_parameters is the number of params for polynominal density function
-#    - parameters : ex.
-#       for 3 parameters p0, p1, p2 and distance from Earth's center x :
-#       density = p0 + p1*x + p2*x*x
-#
-#       CAUTION : the unit of density must be [g/cm3] (CGS unit)
-#
-#    - lines start from '#' or whitespace are treated
-#      as comments
-#
-
 # Dirt
 object box       0 0 -250            0 0 0      100 100 600    surr_earth           ROCK         constant    2.900 
 
diff --git a/resources/earthparams/densities/PREM_IceCube.dat b/resources/earthparams/densities/PREM_IceCube.dat
new file mode 100644
index 00000000..d201eb0a
--- /dev/null
+++ b/resources/earthparams/densities/PREM_IceCube.dat
@@ -0,0 +1,19 @@
+object sphere    0 0 0   0 0 0    6478000    atmo_radius           AIR         radial_polynomial 0 0 0    1    0.000811   # 0.673atm x 1.205e-3(g/cm3 for 1atm)
+object sphere    0 0 0   0 0 0    6374134    iceair_boundary       ICE         radial_polynomial 0 0 0    1    0.762944   # surface of ice, 0.832 x 0.917
+object sphere    0 0 0   0 0 0    6373934    clearice_boundary     ICE         radial_polynomial 0 0 0    1    0.921585   # 200m below ice surface, 1.005 x 0.917
+object sphere    0 0 0   0 0 0    6371324    rockice_boundary      ROCK        radial_polynomial 0 0 0    1    2.650      # surface of bed rock, 1.0 x 2.65
+object sphere    0 0 0   0 0 0    6356000    inner_crust           ROCK        radial_polynomial 0 0 0    1    2.900
+object sphere    0 0 0   0 0 0    6346600    moho_boundary         MANTLE      radial_polynomial 0 0 0    2    2.691     1.08679956050855438e-07   # surface of mantle
+object sphere    0 0 0   0 0 0    6151000    upper_transition      MANTLE      radial_polynomial 0 0 0    2    7.1089    -5.97159001726573544e-07
+object sphere    0 0 0   0 0 0    5971000    middle_transition     MANTLE      radial_polynomial 0 0 0    2    11.2494   -1.26036728927954783e-06
+object sphere    0 0 0   0 0 0    5771000    lower_transition      MANTLE      radial_polynomial 0 0 0    2    5.3197    -2.32867681682624407e-07
+object sphere    0 0 0   0 0 0    5701000    lowermantle_boundary  MANTLE      radial_polynomial 0 0 0    4    7.9565    -1.01649662533354259e-06   1.36199775701391389e-13    -1.19131495406828110e-20
+object sphere    0 0 0   0 0 0    3480000    coremantle_boundary   OUTERCORE   radial_polynomial 0 0 0    4    12.5815   -1.98367603202009108e-07   -8.97421093229181259e-14   -2.13773109929070169e-20
+object sphere    0 0 0   0 0 0    1221500    innercore_boundary    INNERCORE   radial_polynomial 0 0 0    3    13.0885   0   -2.17742748697875934e-13
+
+# IceCube detector: 1km^3 cylinder of ice
+object cylinder  0 0 6372184  0 0 0  564.19 0 1000  icecube        ICE         radial_polynomial 0 0 0    1    0.921585  # same as rockice
+
+
+# center of detector at IceCube
+detector      0 0 6372184
\ No newline at end of file
diff --git a/resources/earthparams/densities/PREM.dat b/resources/earthparams/densities/legacy/PREM.dat
similarity index 100%
rename from resources/earthparams/densities/PREM.dat
rename to resources/earthparams/densities/legacy/PREM.dat
diff --git a/resources/earthparams/densities/PREM_mmc.dat b/resources/earthparams/densities/legacy/PREM_mmc.dat
similarity index 100%
rename from resources/earthparams/densities/PREM_mmc.dat
rename to resources/earthparams/densities/legacy/PREM_mmc.dat
diff --git a/resources/earthparams/densities/mmc.dat b/resources/earthparams/densities/legacy/mmc.dat
similarity index 100%
rename from resources/earthparams/densities/mmc.dat
rename to resources/earthparams/densities/legacy/mmc.dat
diff --git a/resources/earthparams/materials/IceCube.dat b/resources/earthparams/materials/IceCube.dat
new file mode 100644
index 00000000..68ce603e
--- /dev/null
+++ b/resources/earthparams/materials/IceCube.dat
@@ -0,0 +1,24 @@
+ICE   2     # H20
+1000080160    0.8881016     # 0, 88.8% in weight  
+1000010010    0.1118984     # 2H, 11% in weight
+
+ROCK  2     # SiO2
+1000140280    0.4674349     # Si, 33% 
+1000080160    0.5325651     # 20, 66%
+
+INNERCORE 1 # Fe
+1000260560    1.0000000     # Fe56 100%
+
+OUTERCORE 1 # Fe
+1000260560    1.0000000     # Fe56 100%
+
+MANTLE  2   # SiO2
+1000140280    0.4674349     # Si, 33% 
+1000080160    0.5325651     # 20, 66%
+
+AIR 2       # N2 + O2
+1000070140    0.7562326     # N2 78% in volume
+1000080160    0.2437674     # O2 22% in volume
+
+VACUUM 1    # dummy, H
+1000010010    1.0000000     # N2 78% in volume
\ No newline at end of file