From e9b4514a620e7de66708bcfa15871ea1ee90e95d Mon Sep 17 00:00:00 2001 From: Nicholas Kamp Date: Mon, 8 Jan 2024 17:55:07 -0500 Subject: [PATCH 1/4] reorganizing examples and earthparams files --- 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 +++ 12 files changed, 390 insertions(+), 179 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/resources/Examples/Example1/DIS_IceCube.py b/resources/Examples/Example1/DIS_IceCube.py new file mode 100644 index 000000000..924e6d98a --- /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 000000000..b13f7fbf8 --- /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 a30669c0a..d175b5ab1 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 000000000..a3d86f07a --- /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": [ + "
" + ] + }, + "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": [ + "
" + ] + }, + "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": [ + "
" + ] + }, + "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": [ + "
" + ] + }, + "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 f1576cde6..e79103fc6 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 944f89514..cc515459d 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 4d2add8a5..2d196dd9a 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 000000000..d201eb0a9 --- /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 000000000..68ce603e5 --- /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 From 2efe1cd3c0252d9dcb6c2d8d8dca1c4f717ccb33 Mon Sep 17 00:00:00 2001 From: Nicholas Kamp Date: Mon, 8 Jan 2024 18:19:35 -0500 Subject: [PATCH 2/4] starting to add fill interpolation method --- python/LIDarkNews.py | 36 +++++++++++++++++++++++++++++++++++- 1 file changed, 35 insertions(+), 1 deletion(-) diff --git a/python/LIDarkNews.py b/python/LIDarkNews.py index 7f1c58f93..62f06ff54 100644 --- a/python/LIDarkNews.py +++ b/python/LIDarkNews.py @@ -27,6 +27,7 @@ def __init__(self, param_file=None, tolerance=1e-6, interp_tolerance=5e-2, + fill_tables_at_exit=True, **kwargs): # Defines a series of upscattering and decay objects # Each derive from the respective LeptonInjector classes @@ -99,7 +100,12 @@ def __init__(self, def SaveCrossSectionTables(self): for cross_section in self.cross_sections: + if fill_tables_at_exit: + cross_section.FillInterpolationTables() + else: + print('WARNING: Not filling PyDarkNewsCrossSection interpolation tables before exiting. Future updates to DarkNews can lead to inconsistent behavior if new entries are ever added to this table') cross_section.SaveInterpolationTables() + @@ -112,7 +118,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=True): DarkNewsCrossSection.__init__(self) # C++ constructor @@ -246,6 +252,34 @@ def _query_interpolation_table(self,inputs,mode): else: return 0 + # Fills the tables + def FillInterpolationTables(self,total=True,diff=True): + Emin = self.ups_case.Ethreshold + Emax = np.max(self.total_cross_section_table[:,0]) + if total: + E = Emin + while E Date: Tue, 9 Jan 2024 17:29:48 -0500 Subject: [PATCH 3/4] Flush out fill table methods for LIDarkNews --- python/LIDarkNews.py | 118 +++++++++++++++++++++++++++++-------------- 1 file changed, 80 insertions(+), 38 deletions(-) diff --git a/python/LIDarkNews.py b/python/LIDarkNews.py index 62f06ff54..10074e9b3 100644 --- a/python/LIDarkNews.py +++ b/python/LIDarkNews.py @@ -27,7 +27,6 @@ def __init__(self, param_file=None, tolerance=1e-6, interp_tolerance=5e-2, - fill_tables_at_exit=True, **kwargs): # Defines a series of upscattering and decay objects # Each derive from the respective LeptonInjector classes @@ -98,12 +97,14 @@ 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: - cross_section.FillInterpolationTables() - else: - print('WARNING: Not filling PyDarkNewsCrossSection interpolation tables before exiting. Future updates to DarkNews can lead to inconsistent behavior if new entries are ever added to this table') + print('Filling cross section table at %s'%cross_section.table_dir) + num = cross_section.FillInterpolationTables() + print('Added %d points'%num) cross_section.SaveInterpolationTables() @@ -118,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=True): + interpolate_differential=False): DarkNewsCrossSection.__init__(self) # C++ constructor @@ -176,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) @@ -203,20 +203,25 @@ 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 + #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: - # 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 + # # 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 + # Check if we are close enough to interpolate if np.abs(diff) >= self.interp_tolerance: Interpolate = False else: @@ -245,6 +250,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: @@ -252,30 +279,45 @@ def _query_interpolation_table(self,inputs,mode): else: return 0 - # Fills the tables + # Fills the total and differential cross section tables within interp_tolerance def FillInterpolationTables(self,total=True,diff=True): - Emin = self.ups_case.Ethreshold + 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 Date: Tue, 9 Jan 2024 17:48:27 -0500 Subject: [PATCH 4/4] Minor cleanup --- python/LIDarkNews.py | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/python/LIDarkNews.py b/python/LIDarkNews.py index 10074e9b3..4a4f4ac8e 100644 --- a/python/LIDarkNews.py +++ b/python/LIDarkNews.py @@ -203,7 +203,6 @@ def _interpolation_flags(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 @@ -212,15 +211,6 @@ def _interpolation_flags(self,inputs,mode): 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 # Check if we are close enough to interpolate if np.abs(diff) >= self.interp_tolerance: Interpolate = False