{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"WARNING: The TensorFlow contrib module will not be included in TensorFlow 2.0.\n",
"For more information, please see:\n",
" * https://github.com/tensorflow/community/blob/master/rfcs/20180907-contrib-sunset.md\n",
" * https://github.com/tensorflow/addons\n",
"If you depend on functionality not listed there, please file an issue.\n",
"\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/davide/miniconda3/envs/zfit_env/lib/python3.6/site-packages/zfit/util/execution.py:57: UserWarning: Not running on Linux. Determining available cpus for thread can failand be overestimated. Workaround (only if too many cpus are used):`zfit.run.set_n_cpu(your_cpu_number)`\n",
" warnings.warn(\"Not running on Linux. Determining available cpus for thread can fail\"\n"
]
}
],
"source": [
"import numpy as np\n",
"from math import pi\n",
"import matplotlib.pyplot as plt\n",
"from phasespace import Particle\n",
"#from phasespace.kinematics import lorentz_vector, lorentz_boost\n",
"import tensorflow as tf\n",
"import tensorflow_probability as tfp\n",
"import zfit\n",
"\n",
"\n",
"from utils.BtoKll_utils import *\n",
"#from phasespace.kinematics import lorentz_boost\n",
"ztf = zfit.ztf\n",
"ztyping = zfit.util.ztyping\n",
"ztypes = zfit.settings.ztypes"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"zfit.run.check_numerics = False\n",
"zfit.settings.set_verbosity(10)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"WARNING:tensorflow:From /Users/davide/miniconda3/envs/zfit_env/lib/python3.6/site-packages/tensorflow/python/ops/resource_variable_ops.py:435: colocate_with (from tensorflow.python.framework.ops) is deprecated and will be removed in a future version.\n",
"Instructions for updating:\n",
"Colocations handled automatically by placer.\n"
]
}
],
"source": [
"mmu_mass = 105.0\n",
"mB_mass = 5280.0\n",
"mKst_mass = 892.0\n",
"\n",
"minq=2*mmu_mass\n",
"maxq=(mB_mass-mKst_mass)\n",
"\n",
"el1 = Particle('l1', mmu_mass)\n",
"el2 = Particle('l2', mmu_mass)\n",
"\n",
"def modq(min_mass, max_mass, n_events):\n",
" \n",
" min_mass = tf.cast(min_mass, tf.float64)\n",
" max_mass = tf.cast(max_mass, tf.float64)\n",
"\n",
" min_mass = tf.broadcast_to(min_mass, shape=(n_events,))\n",
" \n",
" modq_mass = tfp.distributions.Uniform(low=min_mass, high=max_mass).sample()\n",
" \n",
" return modq_mass\n",
" \n",
"q=Particle('q', modq).set_children(el1, el2)\n",
"Kst = Particle('Kst', mKst_mass)\n",
"\n",
"B = Particle('B', mB_mass).set_children(Kst, q)\n",
"\n",
"#weights, particles = B.generate(n_events=10000)\n",
"weights_np, particles_np = B.generate(n_events=1000000)\n",
"weights, particles = B.generate_tensor(n_events=1000000)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"n_events=100000"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"105.0030608127684 5279.998960900106\n"
]
}
],
"source": [
"Kst=particles_np['Kst']\n",
"l1 =particles_np['l1']\n",
"l2 =particles_np['l2']\n",
"q_np=particles_np['q']\n",
"q2_array=scalar_product_4_np(q_np,q_np)\n",
"q2_range=np.arange(q2_array.min(), q2_array.max(),1000)\n",
"print(np.sqrt(q2_array.min()/4), np.sqrt(q2_array.max())+mKst_mass)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA3MAAAI/CAYAAADdpIDZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdeXzUZ7n//9c9M9n3yWTfN8K+hrCWpViKrW1da+ty1Ko9x+3oz6PnePQcj1/PqvXocWnV2tPTulZbrdZu0AUoZQ8Uyg4JSSAEyErYyTL3748ZQorQTpiQmUnez8eDR5PMB+YKw2Tm3fu678tYaxEREREREZHI4gh1ASIiIiIiIjJ4CnMiIiIiIiIRSGFOREREREQkAinMiYiIiIiIRCCFORERERERkQikMCciIiIiIhKBXKEu4HIej8cWFxeHugwREREREZGQ2LJlS5u1NuOtrgu7MFdcXExNTU2oyxAREREREQkJY0xjINepzVJERERERCQCKcyJiIiIiIhEIIU5ERERERGRCKQwJyIiIiIiEoEU5kRERERERCKQwpyIiIiIiEgEUpgTERERERGJQApzIiIiIiIiESjshoaPRud7+jjadZ6T53o4293HuZ5e+rzgdIDT4cDlMCTEuEiJiyIlLorkWBcup3K4iIiIiMhopjA3zLrO9bDmQCuvHTrBjqYu6lpP036me9B/TlKsi6zkWLKTY8lKjiUnJZaslFjyU+MoSo8nPy2eaJcCn4iIiIjISKUwNwz6vJYX9xznlxsaWV/XTq/XEuNyMCE3maUTsshLjSMnJY7U+Cjiop3ER7twOQy9Xkuf10tPn+XMhV5Onu+h62wPXed66TzbzfGT5znadZ66ujZaTl2gz2v779NhIC8tjiJ3AkXp8RSnJ1DiSWBMVhL5aXE4HCaEfyMiIiIiIhIshbnryFrLy3tb+Pdn9nCw7Qy5KbF84oZSbhqfyZT81CFtlezzWtpOX6Cp8ywNbWdpbD9DQ/tZGjvO8syOo5w429N/bWyUg/LMRMZkJlGRlURFZqJCnoiIiIhIhFGYu05OX+jla0/u4E/bminNSOCBD05n6fis67bXzekwZPlbLmcUuf/i9hNnu6lrPc2B46fZf/w0B1pOsa6unT+8dqT/mrgoJ2Oykxifk8z43GTG5yQzNjuJhBj9MxERERERCTfGWvvWVw2jqqoqW1NTE+oygtJy6jwffmgTB1pO8fklY/j04jKiwvTAkq5zPdS2nObA8VPsP36avcdOsqv5JF3nfCt5xkBxesIbAt743GQyk2IwRqt4IiIiIiJDzRizxVpb9VbXaclliHWc6eauBzdw9MR5Hr2nmhsqMkJd0ptKiYtiRlEaM4rS+r9mreVo13l2N59k99GT7G4+yY4jXTyz42j/NekJ0YzPTWZyfgqT8lKZnJ9CTkqsAp6IiIiIyDBRmBtCPX1ePvOrrTR1nuOXH59FdclftjtGAmMMualx5KbG8bbxWf1fP3m+h71HT7G7uYvdR0+y48hJfrL6YP/BK57EGH+4S/H9Nz+FzKTYUH0bIiIiIiIjmsLcEPrxqjrWH2znO++bErFB7s0kx0ZRXeJ+w/d2vqePPUd9K3evN3Wxo6mLVftauHiwZnZyLJPyU5ic5wt3k/JSSE+MCdF3ICIiIiIycijMDZG61tP86OVa3jE5h/fOyA91OcMmNsrJtMI0phVeatM8293L7uaTvnB3pIvXm07w4p7jXNyemZcax5SCFKYWpDK1II1JeSnERTtD9B2IiIiIiEQmhbkhct/z+4hxOfiX2yaEupSQi492UVXspqr40greqfM97Go+yY6mLrY3nWB70wme3XEM8J3EOTY7yR/uUplWmEqpJ1FjEkRERERE3oTC3BDYc/Qkz+86xt8uqSAjSS2EV5IUG8Xs0nRml6b3f63t9AW2Hz7BtsMneO3QCZ7a1syvNh7yX+9iSv6lcDe1IFXtmSIiIiIiAyjMDYGH1tSTEO3k4/NKQl1KRPEkxrBkXBZLxvkOWfF6LQfbTvPaoUsB78er6/oPWClwxzG1IK1/BW9CbjKxUWrPFBEREZHRSWEuSCfP9/DMjmbeNS2flPioUJcT0RwOQ3lmEuWZSbyvqgCAc9197DjSxbbDnWw7fIKahg7+vL0ZgCinYXxOMlMLUplelMb0wjTy0+I0HkFERERERgWFuSA9vf0o53u8vH9mQahLGZHiop1/cYLm8ZPn+1fvth3u5PEtTTy6vhGAjKQYphemMsMf7ibmpWj1TkRERERGJIW5ID2/6xglngSm5KeEupRRIys5lmUTs1k2MRuA3j4v+46fYmtjJ1sPnWBLYyfLdx0HfKt3E3JT+sPdjKI0slM0+05EREREIp/CXBBOX+hlQ107fzWnSK19IeRyOpiQm8KE3BQ+PMf3tdZTF3jtUCdbDnXyWuMJfrmhkf99tR6A3JTY/rbMGUVpjM9NJsrpCOF3ICIiIiIyeApzQXj1QBvdfV5uHJcZ6lLkMhlJMSydkM3SCb7Vu+5eL3uOnmRLYydbD3WytbGTp18/CkCMy8GU/FSmFaUyozCN6UVpeHRypoiIiIiEOYW5IKyrayM+2snMAfPUJDxFuxxMKUhlSkEq9+A7dfRo1zm2Np5g66FOtjR28vCr9fy07yAARenxTPcHuxmFaVRmJ+HU3DsRERERCSMKc0HY0tjJ1IJUtehFqJyUOG6dHMetk3MAON/Tx84jXf3hbs2BNp587QgACdFOphT4DlapKnYzrTCV5FidXioiIiIioaMwd43OXOhlz9GTfHZxeahLkSESG+WkqthNlX+l1VpLU+e5/nC3pbGT+1fW4rVgDFRmJTGz2E1VsW/vXV6qxiKIiIiIyPBRmLtG2w+fwGthelFaqEuR68QYQ4E7ngJ3PHdMzQN8Id43766TmsYOnnztCL/Y4BuLkJ0cS1VxGlX+1bux2Um4tGorIiIiIteJwtw12tV8EoDJ+akhrkSGU0KMi3nlHuaVewDo81r2HjvpD3edbGno6D9YJSHaydTCVKqKfKt30wrTSIzRU05EREREhobeWV6jfcdPkZEUgzshOtSlSAg5HaZ/LMJH5hYDcOTEOWoaOtjS2ElNQyc/fPkAXgsOA2Ozk5lZnMaMYjdVRWnkpsaF9hsQERERkYilMHeNDhw/xZisxFCXIWEoLzWOvKl5/a2Zp8738NqhE76Vu8YOHt/SxKPrfa2ZuSmx/n16vn13Y7OTdWqmiIiIiAREYe4aeL2W/cdPc1d1QahLkQiQFBvFgjEZLBiTAUBvn5c9R09R09hBTWMnG+vbeWp7MwCJMS6mDWjNnFqQSoJaM0VERETkCoJ6l2iMWQZ8H3ACD1lr/+uy2wuBR4FU/zVfsdY+G8x9hoMjJ85xrqePMVlJoS5FIpDL6WBSfgqT8lP42LyS/lMztzT6DlWpaejkf17aj7W+Ns7xOcn+kQhpVBW5yU6JDfW3ICIiIiJh4JrDnDHGCdwP3AQ0AZuNMU9Za3cPuOyfgN9Za39sjBkPPAsUB1FvWGhoPwNAiSchxJXISDDw1Mx3TvO1Znad6+E1/0iEzQ0dPLb5EI+sawB8bZwzi30nZs4sdlORmYhDrZkiIiIio04wK3PVQK219iCAMeYx4A5gYJizQLL/4xSgOYj7CxtNnecAKHDHh7gSGalS4qJYVJnJospMAHr6vOxuPklNYyc1DR2srWvnj9t8T6fU+Kj+cQgzi91Myksh2qWRCCIiIiIjXTBhLg84PODzJmDWZdd8A1hhjPkckAC8LYj7CxtNnWdxOQxZSTGhLkVGiSingykFqUwpSOXj832tmYc6zrKp3teWubmhgxf3tAAQ43IwtSCVmcVuZpa4mV6YSlJsVIi/AxEREREZasGEuSv1ddnLPr8beMRa+9/GmDnAL4wxE6213jf8QcbcC9wLUFhYGERJw+NwxzlyUmM1EFpCxhhDUXoCRekJvK/KdxBP66kLbGnsYFO9b+/dj1fX8aOVtTgMjMtJ9oW7YjczS9LITNK+OxEREZFIF0yYawIGHueYz1+2UX4cWAZgrV1vjIkFPEDLwIustQ8CDwJUVVVdHgjDTlPnWfJT1WIp4SUjKYZlE3NYNjEHgNMXetl26ASbGjrYXP/GfXfF6fFUFbup9o9FKPEkYIz23YmIiIhEkmDC3GagwhhTAhwB7gI+cNk1h4AlwCPGmHFALNAaxH2GhabOcyz0HzMvEq4SY1zMr/Awv8ID+Pbd7TzSRU1DJ5saOnhpz3Ge2NIEgCcxmqoiX1tmdbGbcTlJWnkWERERCXPXHOastb3GmM8Cy/GNHXjYWrvLGPNNoMZa+xTwd8DPjDH/H74WzI9aa8N+5e3N9Hktbacv6Hh4iThRTgfTCtOYVpjGJxeU4vVaDradZlO9b8/d5oYOnt91DICEaCfTi9KY6V+5m1aQRly0M8TfgYiIiIgMFNScOf/MuGcv+9rXB3y8G5gXzH2Em44z3Xitr6VNJJI5HIbyzCTKM5P4wCzfXtWjXefY3NDJ5npfuPvei755dy6HYWJeCtUlvn13VUVppCVEh/g7EBERERndggpzo1Hb6QsAeBIV5mTkyUmJ4/Ypcdw+JReArrM9bDnU0R/wHlnbwIOvHASgIjPRt++uxDfMPD8tTvvuRERERIaRwtwgtZ5SmJPRIyU+ihvHZnHj2CwAzvf08XpTV39b5tPbm/nNpkMA5KTE+k/MTGNmiZsxmUkaZi4iIiJyHSnMDdKllTm1mMnoExvlpLrETXWJG/DtId177GT/oSobDrbz1HbfobbJsa7+QeYzi9OYlJ9CjEv77kRERESGisLcIF0Mc9ozJwJOh2FCbgoTclP4yNxirLUc7jjHpoYOaho62NTQwct7Lw0zn1aYSnWxm+qSdKYXpRIfrR9BIiIiItdK76QGqe10NzEuB4kx+qsTuZwxhsL0eArT43nvjHwA2k9f8O25a+hgU30HP1pZi/fl2jccqlLtX8FLiY8K8XcgIiIiEjmUSAap7dQFPIkxOuhBJEDpiTEsm5jNsonZAJw638OWxk42+U/MvHioijFQmZXErJJL8+4ykzUCRERERORqFOYGqetcD6laPRC5ZkmxUSyqzGRRZSbgO1Rl2+ET/eHu8S1NPLq+EYAST4K/LdP3SydmioiIiFyiMDdIXed6SIlTmBMZKrFRTmaXpjO7NB2Anj4vu5pPsqm+nU31vkHmv605DPhOzLwY7GaVuCnLSFS4ExERkVFLYW6QTp7vodSTGOoyREasKKeDqQWpTC1I5d4FZXi9lv0tp9hU38HG+g7W1bXzp22+EzPdCdHMLE6juiSdWSVuxuUk49Q4BBERERklFOYG6eS5XpLj9NcmMlwcDsPY7GTGZifzV3N8J2Y2tp/tD3ebGtpZvus4AIkxLmYUpfWv3GkcgoiIiIxkSiWDdPK82ixFQskYQ7EngWJPAnfOLADgaNc5NtV39P+6b/k+QOMQREREZGTTu5pB6Onzcra7j+RYhTmRcJKTEscdU/O4Y2oeAB1nuvtHIWgcgoiIiIxUCnODcPJcDwDJWpkTCWvuhGhunpDNzRPeOA7hYsDTOAQREREZCRTmBuHk+V4A7ZkTiTBXG4ewub6DTRqHICIiIhFKqWQQLq7Mac+cSGS72jiEzf5DVa42DqG62E15psYhiIiISHhQmBuEk+d9YS5Je+ZERpSB4xA+uaD0DeMQNl1hHEJ1sZtZpW5mlaQzNjsJh8YhiIiISAgozA3CmQu+NssEnYYnMqK92TiEDfXtbDzoW70D30r9zGI3s/3hbnyuZt2JiIjI8FAqGYSz3X0AxEdrbpXIaHKlcQhNnWfZeLDDP++unRf3+GbdJcW4qCpOY1apb5D5xLwUopyOUJYvIiIiI5TC3CAozInIRflp8eTPiOc9M/IBONZ1no317Wys72DjwXZW7msFfD8vZhSlMavEzazSdCZrkLmIiIgMEYW5QTjnD3NxCnMicpnslNg3zLprPXWhf9Vu48EOvrNiP+AbZD69MK1/z920wlRio/QzRURERAZPYW4QLq3M6a9NRN5cRlIMt07O4dbJOYBvkPnAcPf9lw5g7QGi/YevXAx304tS9TNGREREAqJ3DINwtqeXaJdDhxuIyKC5E6JZNjGbZRN9g8y7zvZQ09jR35b5wKo6fvhyLS6HYXJ+CrNK06kucVNVlKYTdEVEROSKFOYG4Vx3n/bLiciQSImPYsm4LJaMywLg9IVeahouhbufvXKQH6+qw2FgYl6Kb89dSTozS9yadSkiIiKAwtygnO3uI157W0TkOkiMcbGoMpNFlZkAnO3uZWvjif62zEfXNfKzNfUYA+Oyk/vbMmeVuElLiA5x9SIiIhIKCnODcK67T4efiMiwiI92Mb/Cw/wKDwDne/p47dClcPfrjYf4v7UNAFRmJfWHu+oSNxlJMSGsXERERIaLwtwgnO3u1cEEIhISsVFO5pSlM6csHYALvX3saOpiY30HGw6288SWJn6+vhGA0owEZpWk9w8yz06JDWXpIiIicp0omQzCWa3MiUiYiHE5qSp2U1Xs5jOLy+np87LzSFf/nruntzfzm02HAChKj+/fczer1E1+WnyIqxcREZGhoDA3COd6+nBrb4qIhKEop4NphWlMK0zjbxaW0ee17G4+ycb6djYc7GD5ruP8rqYJgLzUOGaVupldks7s0nQK3HEYo1N6RUREIo3C3CCc7e4jP00rcyIS/pwOw6T8FCblp/CJG0rxei37jp9i48F2NtZ3sGpfK3/YegQYEO5K05lTmk5+msKdiIhIJFCYG4Rz3X3E6jRLEYlADodhXE4y43KS+ei8Eqy11LacZv3BdjYcbH/TcFfgVlumiIhIOFKYG4QLvV5iXApzIhL5jDFUZCVRkZXEX80pxlrLgZbTbLhKuJtd6jtQZbbCnYiISNhQmBuE7t4+YlyOUJchIjLkjDGMyUpijD/ceb2W2tbTrK/zhbuV+1r4/dZLe+4U7kREREJPYW4Quvu8RCvMicgo4HBcCncfmesLdwNX7l7ee7w/3OWnXQx3voCn0zJFRESGh8LcIHT3eol2KsyJyOjjcBgqs5OozL5yuHtpz3Ge2KJwJyIiMpwU5gLU2+fFa9HKnIgIVw53+1tOsaHONwphYLgrcMf1j0GYXZZOXmpciKsXEREZGRTmAtTd5wUU5kRErsThMIzNTmZstu+0zMvD3Qt7jvO4wp2IiMiQUpgLUHevP8ypzVJE5C1dKdztO36qvy3zSuFuTpkv4OUq3ImIiAREYS5A/WFOK3MiIoM2cM7dx94i3BW64/tPylS4ExERuTqFuQBdUJgTERkyVwt3F0chLN91nN/VvDHczSlLZ1aJwp2IiMhFCnMBurhnTnPmRESG3sBwd898X7jbe+zSyt3AcFeUHu/bc1fmW73LSVG4ExGR0SmoMGeMWQZ8H3ACD1lr/+sK19wJfAOwwHZr7QeCuc9QudCjPXMiIsPF4TCMz01mfO6Vw93zu47x25rDAJR4Ephdms5c/567jKSYEFcvIiIyPK45zBljnMD9wE1AE7DZGPOUtXb3gGsqgH8E5llrO40xmcEWHCo6zVJEJHSuFO72HDvZ35b59PZmfrPpEAAVmYnMLUvvb8tMS4gOcfUiIiLXRzArc9VArbX2IIAx5jHgDmD3gGs+Cdxvre0EsNa2BHF/IaUDUEREwofDYZiQm8KE3BQ+cUMpvX1edjWfZF1dO+sPtvO7miYeXd+IMTAuO5k5Zb6Vu5klbpJjo0JdvoiIyJAIJszlAYcHfN4EzLrsmjEAxpi1+Foxv2GtfT6I+wwZjSYQEQlfLqeDKQWpTClI5VOLyuju9fJ60wnW17Wzrq6dX2xo5H9frcdhYFJ+KnNKfSt3M4vTiI/W9nEREYlMwbyCmSt8zV7hz68AFgH5wBpjzERr7Yk3/EHG3AvcC1BYWBhESddPd18foJU5EZFIEO1yUFXspqrYzeeWVHC+p4/XDp1gfV0b6w+289Cag/xkdR1RTsOU/FTm+NsypxemERvlDHX5IiIiAQkmzDUBBQM+zwear3DNBmttD1BvjNmHL9xtHniRtfZB4EGAqqqqywNhWLi4Mhfj0ou8iEikiY1y9gc2gLPdvdQ0dLL+oG/l7v6Vtfzw5VqiXQ5mFKb1XzslP1X/E09ERMJWMGFuM1BhjCkBjgB3AZefVPlH4G7gEWOMB1/b5cEg7jNkNGdORGTkiI92sWBMBgvGZABw8nwPm+s7WO/fc/e9F/fz3RcgLspJVXEac8s8zClLZ2JuMi6124uISJi45jBnre01xnwWWI5vP9zD1tpdxphvAjXW2qf8ty01xuwG+oAvW2vbh6Lw4XZpZU4v4iIiI01ybBRLxmWxZFwWACfOdrPhYEd/W+a3nt8LQFKMi+oS3wDz2aXpjM9JxuG40q4DERGR6y+oXd/W2meBZy/72tcHfGyBL/p/RbSePl/3Z5T+j6yIyIiXGh/NsonZLJuYDUDrqQtsOOhbtVtf185Le1v810Uxq8Tdv3JXkZmIMQp3IiIyPHSEV4D6vL6VOZdTL9IiIqNNRlIMt03J5bYpuQAc7Trna8n0n5a5fNdxADyJ0cz2n5Q5t8xDcXq8wp2IiFw3CnMB6vX6VuZcaqcRERn1clLiePf0fN49PR+Awx1n/cHO15b59OtHAchOju0/TGVOaToF7vhQli0iIiOMwlyA+vxhzqkwJyIilylwx1PgjufOmQVYa6lvO9M/wPyV/a08+doR/3VxzC31MLfct3KXkRQT4spFRCSSKcwF6OKeOZdDe+ZEROTqjDGUZiRSmpHIh2YXYa1l//HTrK9rY11dO8/tPMpvaw4DUJmVxJyydOaVe5hV6iY5NirE1YuISCRRmAvQxT1zWpkTEZHBMMZQmZ1EZXYSH51XQp/XsvNIF+v8bZmPbT7EI+sacBiYnJ/KPP+q3YwiDTAXEZE3pzAXIO2ZExGRoeB0GKYUpDKlIJVPLSrjQm8fWxtPsK6ujbW1bfxk9UHuX1lHtMtBVVEa88o9zC1LZ1JeimbciYjIGyjMBajPa3EYNE9IRESGVIzL2X9Iyt8treTU+R42N3SwtradtbVt3Ld8H+CbcTerNJ25/rbMMVkagyAiMtopzAWo12u1X05ERK67pNgobhybxY1jfQPM205f6D8pc11dOy/uuTgGIcYf7HxtmTopU0Rk9FGYC1Bvn1f75UREZNh5Et84466p8yzrattZ6w93T21vBqDQHc/csnTm+tsyPYk6KVNEZKRTmAuQb2VOYU5EREIrPy2eO2deGoNQ23KatbVtrK1r55kdR3lss++kzLHZScwt8zCvPJ3qEjdJOilTRGTEUZgLUJ/X4nIqzImISPgwxlCRlURFlu+kzN4+LzubT7K2to31de38amMjD6+tx+kwTM5PYV6Zb8bd9EKdlCkiMhIozAWo12txas+ciIiEMZfTwdSCVKYWpPKZxeWc7+lj66HO/rbMH6+u40cra4lxOZhZ7O6fcTcpL0VbCUREIpDCXID6+tRmKSIikSU2ysncMg9zyzx8Cd9JmZvqfSdlrqvznZR53/J9JMW6mF2azjx/uCvP1EmZIiKRQGEuQD1eHYAiIiKRLSk2iiXjslgy7tJJmevq2llX6ztM5YXdvpMys5JjmFfu4YYKD/PKPGQmx4aybBERuQqFuQD1eS1R2jMnIiIjiCcxhtun5HK7/6TMwx1nWVvbxpraNlbubeEPW48AUJmV1B/uqkvcJMTo7YOISDjQT+MA+fbMKcyJiMjIVeCO567qQu6qLsTrtew+epJXa9tYW9vWf5iKy2GYXpjG/AoP88o9TMlPweXUnnIRkVBQmAuQb8+cXqxERGR0cDgME/NSmJiXwt8sLON8Tx9bGjtZc8AX7r734n6++8J+kmJczC5L97Vklnso9SRov52IyDBRmAtQr/bMiYjIKBYb5WReuS+wAXSc6WZ9XTuv1rbxam1r/3673JRY5pV7mF/hO3glI0nDy0VErheFuQD1as6ciIhIP3dCNLdOzuHWyTkAHGo/y5raVtbWtrFi93Ee39IE+IaXX1y1m1WSTly05tuJiAwVhbkA9Xk1mkBERORqCtPj+WB6ER+cVUSf17Krucu3anegjUfXNfKzNfVEOx1ML0plfrmH+RUZmm8nIhIkhbkA9WrPnIiISECcDsPk/FQm56fy6UXlnOvuY3NDh++kzANtfGfFfr6zYj/JsS7mlnmYV+FhfrmH4vR47bcTERkEhbkA9ek0SxERkWsSF+1kwZgMFozJ4B+Bdv98u1cPtPFqbRvP7zoGQF5qnH/VzsPcsnTSE7XfTkTkzSjMBajH6yUmSn9dIiIiwUpPjOG2KbncNiUXay0N7Wd9IxAOtPHczqP8tuYwABNyk5lf4WFBRQYzitKIjdJ+OxGRgZROAqSVORERkaFnjKHEk0CJJ4EPz/btt9txpItXD7Sy5kAbD79az09XHyQ2ysGsEt8IhIVjMijPTFRLpoiMegpzAdKeORERkevP6TBMLUhlakEqn72xgjMXetlY384r+9tYc6CVf3tmD//2zB6yk2O5ocLDDWMymF/uwZ0QHerSRUSGncJcgHwrc6GuQkREZHRJiHFx49gsbhybBcCRE+dYs9+3andxBIIxMCkvxRfuKjKYXphGtEsv2iIy8inMBajPqs1SREQk1PJS47irupC7qgv7WzLX7G/llQOt/GT1Qe5fWUdCtJM5ZencUJHBDRUeSjwJaskUkRFJYS5AXmtx6IVAREQkbAxsyfzckgpOne9hfV07aw608cqBVl7c0wL4AuCCMb6DVOaWeUiJjwpx5SIiQ0NhLkDWojAnIiISxpJio1g6IZulE7IBaGw/w5oDvr12T28/ym82HcZhYEpBKjdUZLCgwsPUglRc2kchIhFKYS5AvpW5UFchIiIigSpKT6AoPYEPzS6it8/LtsMneMUf7n708gF+8NIBkmJczC1P94e7DArT40NdtohIwBTmAqQ2SxERkcjlcjqoKnZTVezmizeNoetsD2vrfMHulf1tLN91HICi9HgW+PfazSlLJylWLZkiEr4U5gLk9aLN0yIiIiNESlLX6kAAACAASURBVHwUt0zK4ZZJOVhrqW87wyv+UzJ/v7WJX2xoxOkwTC9MZeGYDBaOyWRCbjIOtemISBhRmAuQ12o0gYiIyEhkjKE0I5HSjEQ+Oq+E7l4vWw91suZAK6v3t/KdFfv5zor9eBKjWVCRwcJK32y79MSYUJcuIqOcwlyA1GYpIiIyOkS7HMwuTWd2aTpfvnksbacvsOZAK6v2tbJyXwt/eO0IxsDkvBTfql1lJlMLUjXCSESGncJcgLxWbZYiIiKjkScxhndNy+dd0/Lp81p2Huli9X7fqt2PVtbyg5drSYmLYn6Fx9+SmUFWcmyoyxaRUUBhLkBer06zFBERGe2cDsOUglSmFKTyt0sq6Drbw5raVlbv84W7Z14/CsC4nOT+YDejKI1ol/ZqiMjQU5gLkG/PnNKciIiIXJISH8U7Jufyjsm5WGvZe+wUq/a1snp/Cw+tOchPVteRGONiblk6Cyt94S4/TeMPRGRoKMwFyKuh4SIiIvImjDGMy0lmXE4yn1pUxukLvayrbWP1ft9+uxW7feMPyjISWDgmk4WVGcwqcRMb5Qxx5SISqRTmAuS1FmU5ERERCVRijIulE7JZOiEbay11rWf699r9cmMjD6+tJzbKd9jKxZbMEk+C9uiLSMAU5gJktTInIiIi18gYQ3lmIuWZiXx8fgnnuvvYWN/Oqn2tvLK/lf/3590AFLrjWTgmg8VjM5hT6iEuWqt2InJ1QYU5Y8wy4PuAE3jIWvtfV7nuvcDjwExrbU0w9xkqfToARURERIZIXLSTRZWZLKrMBOBQ+1lWH/AdpHJxaHmMy8GcsnQWV2ayuDKTwnTttRORN7rmMGeMcQL3AzcBTcBmY8xT1trdl12XBPwtsDGYQkPNay0OpTkRERG5DgrT4/lwehEfnl3Ehd4+NtV3sHKvb67dvzy1i39hF2UZCSyuzOTGsZlUFbt1QqaIBLUyVw3UWmsPAhhjHgPuAHZfdt2/At8GvhTEfYWc2ixFRERkOMS4nNxQkcENFRl8/bbx1LedYeXeFlbua+Hn6xt56NV6EmNczC/3sHhsBosrM8nUXDuRUSmYMJcHHB7weRMwa+AFxphpQIG19mljTESHOa9Vm6WIiIgMvxJPAiXzS7hnfglnLvSyrq6dl/e2sGpfC8/vOgbAhNxkbhzra9ucWpCqcUoio0QwYe5KPyVs/43GOIDvAR99yz/ImHuBewEKCwuDKOn66bNWK3MiIiISUgkxLm4an8VN47P659qt3NfCyr0t3L+ylh++XEtafJT/EJVMFlRkkJYQHeqyReQ6CSbMNQEFAz7PB5oHfJ4ETARW+Y/YzQaeMsbcfvkhKNbaB4EHAaqqqixhxlqrNksREREJKwPn2n16UTldZ3tYfaCVVXtbWLW/lT9ua8ZhYFphmn/VLoPxOckafSAyggQT5jYDFcaYEuAIcBfwgYs3Wmu7AM/Fz40xq4AvReJpltYfLxXmREREJFylxEdx+5Rcbp+SS5/X8nrTCVbua2Xl3hbuW76P+5bvIys5hsX+UzTnV3hIjNGUKpFIds3PYGttrzHms8ByfKMJHrbW7jLGfBOosdY+NVRFhprXn+bUfi4iIiKRwOkwTCtMY1phGl+8aQwtp86zal8rq/a18MzrR3ls82GinIbZpeksGZvJknFZFLg1+kAk0hhrw6ursaqqytbUhNfiXXevlzH/9BxfvrmSzywuD3U5IiIiItesp89LTUMnL+89zkt7WjjYdgaAyqwklozzBTsdoiISWsaYLdbaqre6TmvrAbi4MqcuSxEREYl0UU7fMPI5Zel87dbxHGw9zUt7Wnhxz3F++spBHlhVR3pCNIvHZvK2cZncUJFBgtoxRcKSnpkBuBjmnEpzIiIiMsKUZiRSmpHIJxeUcuJsN6v3t/LinhaW7zrGE1uaiHY6mF2Wztv8q3Z5qXGhLllE/BTmAuDVASgiIiIyCqTGR3PH1DzumJpHT5+XzQ0dvLSnhZf2HOfrf9rF1/+0i7HZSbxtXBZLxmUyJT8Vh9oxRUJGYS4AarMUERGR0SbK6WBumYe5ZR7+6dZx1LWe4aU9vn12D6yq5Ucra/EkxnDj2AyWjMvihgoP8dF6aykynPSMC4DXe/E0S6U5ERERGX2MMZRnJlKemchfLyyj80w3q/a38OKeFp7bcYzf1TQR7XIwtyydJeOyWDI2k1y1Y4pcdwpzAbjYZqlTnUREREQgLSGad03L513T8n3tmPUdvLinhZf2Huef/7iTfwYm5iWzdHw2SydkUZmVpGHlIteBwlwANGdORERE5MqinA7mlnuYW+7hn98xjrrW07ywu4UXdh/juy/s57sv7KfQHc9N47NYOj6LqmK3/ge5yBBRmAvApT1z+sEjIiIicjW+dswkyjOT+NSiMlpOnufFPS2s2H2MX6xv5H9frcedEM2SsZksnZDNDRUeYqOcoS5bJGIpzAXA6jRLERERkUHLTI7lA7MK+cCsQk5f6GX1vlZW7D7G87uO8fiWJuKinNxQ4WHphGyWjM0kLSE61CWLRBSFuQD0edVmKSIiIhKMxBgXt07O4dbJOXT3etlY386KXcd5YfdxVuw+jtNhmFmcxk3js1k6PosCd3yoSxYJewpzAejfM6c0JyIiIhK0aJeDGyoyuKEig2/eMYEdR7pYses4K3Yf41+f3s2/Pr2bcTnJLB2fxdIJWYzPSdZ2F5ErUJgLgNosRURERK4PYwyT81OZnJ/Kl26upKHtjH+17hg/ePkA33/pAHmpcdw0Pou3T8zWASoiAyjMBUCnWYqIiIgMj2JPAp9cUMonF5TSdvoCL/sPUPn1pkM8sq4BT2IMSyf4gt3s0nSinI5QlywSMgpzAejT0HARERGRYedJjOHOmQXcObOA0xd6Wbm3hed3HuOPrx3h1xsPkRofxdvG+YLd/AoPMS6djCmji8JcAC4ODdeeOREREZHQSIxxcduUXG6bksv5nj5W729l+c5jLN91jCe2NJEY4+LGsZm8fWI2CysziI/W21wZ+fSvPABWbZYiIiIiYSM2ysnNE7K5eUI23b1e1tW18fzOY6zYfZyntjcTG+Vg0ZhM3j4pmxvHZpIUGxXqkkWuC4W5AHh1AIqIiIhIWIp2OVhUmcmiykz+7Z1eNjV08PzOY75fu44R7XQwrzydt0/M4abxWZplJyOKwlwALh6AoignIiIiEr5cTgdzyzzMLfPwjdsm8NrhTp7bcYzndh5j5b7XcT5pmF3q5pZJOSybkE16YkyoSxYJisJcAC6OJtB8ExEREZHI4HAYZhS5mVHk5mu3jmNX80me23mU53Yc42tP7uTrf9rF3LJ03jE5h5snZJMarxU7iTwKcwGw2FCXICIiIiLXyBjDxLwUJual8KWllew9doqnX2/m6deP8g+/38HXntzJ/AoPt07KYemEbFLitMdOIoPCXAAurcyFtg4RERERCY4xhnE5yYzLSeZLSyvZ1XySP7/ezDOvH+XLT7zOV5/cwYKKDN4xJYe3jcvS4SkS1hTmBkFZTkRERGTkGLhi95VlY9ne1MUz/mD30t4W3+EqYzK4dbIv2CXE6K2zhBf9ixwE7ZkTERERGZmMMUwtSGVqQSr/+PZxvHb4BE+/3syzO46yYvdxYlwObhybyTsm57J4rObYSXjQv8IAWG2ZExERERk1fIenpDGjKI1/vnU8NY2d/mDnOxkzPtrJ0vFZ3DE1j/kVHqKcjlCXLKOUwlwALh6AonU5ERERkdHF4TBUl7ipLnHzL7dNYFN9B09t963Y/XFbM+6EaG6dlMMdU3OZUZSmTi4ZVgpzg6DnpoiIiMjo5XQY5pSlM6csnf93+wRW72/lT9uO8Luaw/xiQyP5aXHcPiWXd07LY0xWUqjLlVFAYS4AarMUERERkYGiXQ5uGp/FTeOzOH2hlxW7jvHHbc389JWDPLCqjrHZSdwxNY/bp+aSlxoX6nJlhFKYC8DFLKeVORERERG5XGKMi3dPz+fd0/NpPXWBZ15v5k/bm/nW83v51vN7qS5xc8fUXG6ZmENagoaTy9BRmBsEo11zIiIiIvImMpJi+Oi8Ej46r4TG9jM8ta2ZP247wtee3Mk3ntrFwjGZvGd6HjeOyyTG5Qx1uRLhFOYCYNVnKSIiIiKDVJSewOeWVPDZG8vZ1XySP207wp+2NfPinuOkxkdx2+Rc3jMjnyn5KTo4Ra6JwlwA+qOcnmMiIiIiMkgDh5P/w7KxvFrbxu+3Xjo4pTwzkXdPz+Pd0/LJTokNdbkSQRTmBkFZTkRERESC4XI6WFSZyaLKTLrO9fDsjqP8fksT335+H/ct38f8cg/vmZ7PzROyiYtWG6a8OYW5AKjLUkRERESGWkpcFHdXF3J3dSENbWf4w9Ymfr/1CF/47TYSY1zcMimb90zPp7rErTZMuSKFuYD4h4brSSQiIiIi10GxJ4EvLq3kC28bw8b6Dn6/tYlnXj/K72qaKHTHc2dVPu+dUaA2THkDhbkAXFyZU5QTERERkevJMWAw+TfvmMDzO4/xeE0T31mxn+++sJ8bx2by/pmFLK7MwOV0hLpcCTGFuUHQwpyIiIiIDJf46Evz6xrazvC7msM8vqWJF/fUkJkUw3tn5PP+mQUUpSeEulQJEYW5AGjLnIiIiIiEUrEngb9fNpYv3jSGlfta+e3mQ/xkdR0PrKpjTmk6d1UXcPOEbGKjdGjKaKIwF4BLbZZamhMRERGR0HE5Hdw0PoubxmdxrOs8v9/axGObD/H5x7aREhfFu6blcVd1AWOzk0NdqgwDhblBUJuliIiIiISL7JRYPrO4nE8tLGPDwXYe23yYX288xCPrGphZnMaHZhexbGI2MS6t1o1UCnMBsJpNICIiIiJhyuEwzC33MLfcQ+eZbp7Y0sQvNzby+ce24UmM5v0zC7i7upD8tPhQlypDTGEuABejnBbmRERERCScpSVE88kFpXx8fglratv4xfpGfryqjh+vquPGsZl8aHYRCyoycDj0znYkCCrMGWOWAd8HnMBD1tr/uuz2LwKfAHqBVuAea21jMPcZUvo3LyIiIiIRwOEwLByTwcIxGRw5cY7fbDzEY5sP8eKeFgrd8XxwViHvqyrAnRAd6lIlCNc8nMIY4wTuB94OjAfuNsaMv+yy14Aqa+1k4Ang29d6f6GkLksRERERiVR5qXF86eZK1n1lCT+4exrZybH853N7mf2fL/H3T2xn77GToS5RrlEwK3PVQK219iCAMeYx4A5g98ULrLUrB1y/AfhQEPcXMtbfaKnTLEVEREQkUkW7HNw+JZfbp+Sy79gpHl3fwB+2NvG7mibmladzz7wSFldmqgUzggQzNj4PODzg8yb/167m48BzQdxfyOk0SxEREREZCSqzk/iPd01i/VeW8PfLKqlrOcPHH61hyXdX8/P1DZy50BvqEiUAwYS5K0WbKzYkGmM+BFQB913l9nuNMTXGmJrW1tYgSrpO1GYpIiIiIiNQWkI0n15Uzpp/WMwP7p5GclwUX//TLub850v857N7OHLiXKhLlDcRTJtlE1Aw4PN8oPnyi4wxbwO+Biy01l640h9krX0QeBCgqqoq7KKTTrMUERERkZEsynmpBXNLYycPr63noVd9v26ZlMPfLCxlQm5KqMuUywQT5jYDFcaYEuAIcBfwgYEXGGOmAT8FlllrW4K4r5C6eACKUZ+liIiIiIxwM4rSmFGUxpET53h0XQO/3niIP29vZsGYDD61sIzZpW69Lw4T19xmaa3tBT4LLAf2AL+z1u4yxnzTGHO7/7L7gETgcWPMNmPMU0FXHEL6NysiIiIio0VeahxfvWUca79yI1++uZLdzV3c/bMNvOuBdTy/8xheb9g11I06Qc2Zs9Y+Czx72de+PuDjtwXz54cLq01zIiIiIjJKpcRF8ZnF5Xx8fglPbGniwVcO8je/3EJZRgJ/vbCMd07NI9oVzFEccq30tx6A/jbL0JYhIiIiIhIysVFOPjS7iJf/biE/uHsaMS4nf//E6yy6byW/3NDIhd6+UJc46ijMDYLaLEVERERktHP5D0t55m/n88jHZpKdEss//XEni+9bpVA3zBTmAqAmSxERERGRNzLGsKgyk99/ai4/v6earAGh7lcbG+nu9Ya6xBFPYS4A1mo4gYiIiIjIlRhjWDAmgz8MCHVfe3Ini7/jC3U9fQp114vC3CCozVJERERE5MoGhrpH76kmMzmGrz25k6Xfe4VnXj86YIFEhorCXAD0z05EREREJDDGGBb6Q93/fqSKKKfhM7/eyjvvX8u6urZQlzeiKMwFQqdZioiIiIgMijGGJeOyeO7zC7jvvZNpOXWBD/xsIx95eBO7m0+GurwRQWFuEDTpXkRERERkcJwOw/uqClj5pUV89ZaxbDt8glt/uIYvP76dllPnQ11eRFOYC4CGhouIiIiIBCc2ysm9C8p45cuL+eQNpfxx2xFu/M5qHnylTidfXiOFuQBoaLiIiIiIyNBIiY/iq7eMY/kXFlBd4uY/nt3Lsv95hZV7W0JdWsRRmAtAf5hTmhMRERERGRKlGYk8/NGZ/N9HZwLwsUc2c88jmznccTbElUUOhblBMFqbExEREREZUovHZvL8Fxbw1VvGsvFgO0u/9wo/e+UgvZpP95YU5gKgHXMiIiIiItdPtMvBvQvKeOGLC5lbls6/P7uHdz6wlp1HukJdWlhTmAvAxQGHarMUEREREbl+clPjeOgjVdz/gekc67rAHfev5T+e3cP5nr5QlxaWFOZERERERCRsGGO4dXIOL31xIXdW5fPgKwe57YevapXuChTmAqA2SxERERGR4ZUSH8V/vnsyj95TTde5Ht71wFoeWFVLn1fvzi9SmAuATrMUEREREQmNhWMyWP6FBSwdn823n9/H+3+6nqZOnXgJCnODotMsRURERESGX1pCND/6wDT+5/1T2XfsFO/44aus3Ke5dApzAdFSroiIiIhIKBljeOe0PP78ufnkpMTxsf/bzH+v2Deq2y4V5gKgNksRERERkfBQ7EngyU/P5f1VBfzw5Vr+6uGNnDjbHeqyQkJhbhAU5kREREREQi82ysm33juZb79nMpvrO3n3A+toaDsT6rKGncJcAEbvwq2IiIiISPi6c2YBv/zELDrPdvOuB9ayuaEj1CUNK4W5APS3WeoAFBERERGRsFJd4ubJT88jLT6aD/5sIyt2HQt1ScNGYS4A1r82pzZLEREREZHwU+xJ4A+fnsv43GQ+9aut/Hl7c6hLGhYKc4OgLCciIiIiEp5S46P55SdmMaMojc8/9hqP1xwOdUnXncJcAKw2zYmIiIiIhL3EGBePfqyaeeUe/uH3r/PcjqOhLum6UpgLwMUspzZLEREREZHwFhft5MEPVzGtMI3PP7aNdbVtoS7pulGYGxSlORERERGRcBcX7eThj8ykxJPAJ39ew/7jp0Jd0nWhMBcAqz5LEREREZGIkhIfxSP3zCQu2sVf/2ILJ8/3hLqkIacwNwhqsxQRERERiRw5KXE88MHpHO44yxd/u33ELdIozA2CspyIiIiISGSpLnHz1VvG8eKe4zy2eWSdcKkwF4ARFuBFREREREaVj84tZm5ZOv/+zB6aOs+GupwhozAXgEtDw7U2JyIiIiISaRwOw7feMxmvtXzzz7tDXc6QUZgbBEU5EREREZHIVOCO59OLylix+zib6jtCXc6QUJgLgNosRUREREQi38fnl5KdHMt9y/eGupQhoTAXgIthTl2WIiIiIiKRKy7ayScXlLK5oZPXm06EupygKcwNglGjpYiIiIhIRHtfVT4J0U4eWdcQ6lKCpjAXAHVZioiIiIiMDMmxUdw6OYcVu45zobcv1OUERWEuABeHC6rNUkREREQk8r19Yg6nL/SyrrY91KUERWEuAFqZExEREREZOeaUpeNyGDY3RPaplkGFOWPMMmPMPmNMrTHmK1e4PcYY81v/7RuNMcXB3F+oaWVORERERCTyxUY5GZuTxLbDkX0IyjWHOWOME7gfeDswHrjbGDP+sss+DnRaa8uB7wHfutb7CyktzYmIiIiIjChjspJobD8b6jKCEszKXDVQa609aK3tBh4D7rjsmjuAR/0fPwEsMSby1rfSEqKZkp9CtEtdqSIiIiIiI0FmUiytpy70n48RiVxB/N484PCAz5uAWVe7xlrba4zpAtKBtiDud9jdND6Lm8ZnhboMEREREREZIqUZCYzNSeJsdx8JMcHEotAJpuorrbBdHmsDuQZjzL3AvQCFhYVBlCQiIiIiIvLW7qwq4M6qglCXEZRg+gabgIHffT7QfLVrjDEuIAX4iyNjrLUPWmurrLVVGRkZQZQkIiIiIiIyOgQT5jYDFcaYEmNMNHAX8NRl1zwFfMT/8XuBl20kN6WKiIiIiIiEiWtus/TvgfsssBxwAg9ba3cZY74J1FhrnwL+F/iFMaYW34rcXUNRtIiIiIiIyGgX1E4/a+2zwLOXfe3rAz4+D7wvmPsQERERERGRv6Sz9kVERERERCKQwpyIiIiIiEgEUpgTERERERGJQApzIiIiIiIiEUhhTkREREREJAKZcBv7ZoxpBRpDXcdlPEBbqIsQQI9FONFjET70WIQPPRbhQ49F+NBjET70WISPt3osiqy1GW/1h4RdmAtHxpgaa21VqOsQPRbhRI9F+NBjET70WIQPPRbhQ49F+NBjET6G6rFQm6WIiIiIiEgEUpgTERERERGJQApzgXkw1AVIPz0W4UOPRfjQYxE+9FiEDz0W4UOPRfjQYxE+huSx0J45ERERERGRCKSVORERERERkQikMOdnjFlmjNlnjKk1xnzlCrfHGGN+6799ozGmePirHPmMMQXGmJXGmD3GmF3GmM9f4ZpFxpguY8w2/6+vh6LW0cIY02CM2eH/u665wu3GGPMD/3PjdWPM9FDUOdIZYyoH/JvfZow5aYz5wmXX6LlxnRhjHjbGtBhjdg74mtsY84Ix5oD/v2lX+b0f8V9zwBjzkeGremS6ymNxnzFmr/9n0JPGmNSr/N43/Xkmg3OVx+IbxpgjA34O3XKV3/um77tkcK7yWPx2wOPQYIzZdpXfq+fFELrae9nr9ZqhNkvAGOME9gM3AU3AZuBua+3uAdd8Gphsrf0bY8xdwLuste8PScEjmDEmB8ix1m41xiQBW4B3XvZYLAK+ZK19R4jKHFWMMQ1AlbX2irNQ/C/UnwNuAWYB37fWzhq+Ckcf/8+sI8Asa23jgK8vQs+N68IYswA4DfzcWjvR/7VvAx3W2v/yvxlNs9b+w2W/zw3UAFWAxfczbYa1tnNYv4ER5CqPxVLgZWttrzHmWwCXPxb+6xp4k59nMjhXeSy+AZy21n7nTX7fW77vksG50mNx2e3/DXRZa795hdsa0PNiyFztvSzwUa7Da4ZW5nyqgVpr7UFrbTfwGHDHZdfcATzq//gJYIkxxgxjjaOCtfaotXar/+NTwB4gL7RVyVu4A9+Lh7XWbgBS/T/I5PpZAtQNDHJyfVlrXwE6LvvywNeFR/G9WF/uZuAFa22H/8X4BWDZdSt0FLjSY2GtXWGt7fV/ugHIH/bCRqGrPC8CEcj7LhmEN3ss/O9X7wR+M6xFjVJv8l72urxmKMz55AGHB3zexF8GiP5r/C8YXUD6sFQ3ShlfK+s0YOMVbp5jjNlujHnOGDNhWAsbfSywwhizxRhz7xVuD+T5I0PrLq7+oqznxvDJstYeBd+LN5B5hWv0/Bh+9wDPXeW2t/p5JkPjs/6W14ev0kqm58XwugE4bq09cJXb9by4Ti57L3tdXjMU5nyutMJ2ef9pINfIEDHGJAK/B75grT152c1bgSJr7RTgh8Afh7u+UWaetXY68HbgM/5WjoH03BhGxpho4Hbg8SvcrOdG+NHzYxgZY74G9AK/usolb/XzTIL3Y6AMmAocBf77CtfoeTG87ubNV+X0vLgO3uK97FV/2xW+9qbPDYU5nyagYMDn+UDz1a4xxriAFK6ttUDegjEmCt8//l9Za/9w+e3W2pPW2tP+j58FoowxnmEuc9Sw1jb7/9sCPImvPWagQJ4/MnTeDmy11h6//AY9N4bd8Ystxf7/tlzhGj0/hon/oIB3AB+0VzkQIICfZxIka+1xa22ftdYL/Oz/b+/+gy2t7/qAvz/dNdjElkSy/gLWpQWdErVpskI6asoUgxBt1lgoSxylFmdFpVPHUUNsxRTHGbAqTRtapYEOQQ2k2OidcSOmYlvrJJQlRpMFqRu6hhsy+cEypEgQN/n0j3NwTk7u3XuW+/O55/WaubPPj+9z7/fc7z7nPJ/7+f7I0r9j98UGGT+zfleSu5Yr475Ye8s8y67LZ4ZgbuT+JOdU1Vnjv3rvT7IwVWYhyXMzylya0UBrf0VaY+N+3bcmeai7f3GZMl/x3HjFqjovo//Hj29cLedHVb1oPHg3VfWiJBcl+dBUsYUk31sjr8pogPXHNriq82TZv7C6Nzbc5OfClUl+c4ky9yS5qKpeMu5udtH4GGuoqi5O8sYkr+vup5cpM8v7Gas0NWb69Vn6dzzLcxdr41uT/El3Ly510n2x9k7wLLsunxk7V1/l4RvPfnVNRr+sHUlu6+7DVXV9kkPdvZBRo9xRVUcyysjt37wab2vflOR7knxwYgrdn0yyO0m6+5cyCqZ/sKqOJ/lMkv0C63Xz5UneNY4Pdib5te7+7aq6Ovmr9jiY0UyWR5I8neT7Nqmu215VvTCj2d9+YOLYZFu4N9ZJVb0jyQVJXlpVi0l+OskNSd5ZVVcl+UiSy8Zl9ya5uru/v7uPVdXPZPTwmiTXd7deHauwTFu8KckpSd4zfr9633j26a9K8rbufm2WeT/bhJewbSzTFhdU1csz6hp2NOP3q8m2WO65axNewraxVFt0961ZYoy1+2LdLfcsuy6fGZYmAAAAGCDdLAEAAAZIMAcAzwUluAAAHB1JREFUADBAgjkAAIABEswBAAAMkGAOAABggARzAAAAAySYAwAAGCDBHAAAwAAJ5gAAAAZIMAcAADBAgjkAAIABEswBAAAMkGAOAABggARzAAAAAySYAwAAGCDBHAAAwAAJ5gAAAAZIMAcAADBAgjkAAIABEswBAAAMkGAOAABggARzAAAAAySYAwAAGCDBHAAAwAAJ5gAAAAZIMAcAADBAgjkAAIABEswBAAAMkGAOAABggARzAAAAAySYAwAAGCDBHAAAwAAJ5gAAAAZIMAcAADBAOze7AtNe+tKX9p49eza7GgAAAJvigQce+FR371qp3JYL5vbs2ZNDhw5tdjUAAAA2RVX92SzldLMEAAAYIMEcAADAAAnmAAAABkgwBwAAMECCOQAAgAESzAEAAAyQYA4AAGCABHMAAAADJJgDAAAYoJ2bXQGSPdf+1gnPH73h2zeoJgAAwFDIzAEAAAyQzNwGWSn7tlk/V9YPAACGSTA3ACcKyARjAAAwnwRzA7fajJ9AEQAAhkkwx7J00QQAgK1LMLdGNmtMHAAAMJ9mms2yqi6uqoer6khVXbvE+VOq6q7x+fuqas/4+BdV1e1V9cGqeqiq3rS21QcAAJhPK2bmqmpHkpuTvCbJYpL7q2qhux+cKHZVkie6++yq2p/kxiSXJ7ksySnd/fVV9cIkD1bVO7r76Fq/EDae8XYAALB5ZsnMnZfkSHc/0t3PJrkzyb6pMvuS3D7evjvJhVVVSTrJi6pqZ5K/nuTZJJ9ek5oDAADMsVmCudOTPDqxvzg+tmSZ7j6e5Mkkp2UU2P15ko8l+UiSn+/uY6usMwAAwNybZQKUWuJYz1jmvCSfTfJVSV6S5Per6r919yOfd3HVgSQHkmT37t0zVImtbrUTwuimCQAAJzZLMLeY5MyJ/TOSPLZMmcVxl8pTkxxL8oYkv93df5nkE1X1B0n2Jvm8YK67b0lyS5Ls3bt3OlBkDhmPBwAAJzZLMHd/knOq6qwkH02yP6MgbdJCkiuTvDfJpUnu7e6uqo8k+YdV9StJXpjkVUn+7VpVnvlk/TsAAJhhzNx4DNw1Se5J8lCSd3b34aq6vqpeNy52a5LTqupIkh9N8tzyBTcn+ZIkH8ooKPzP3f3Ha/waAAAA5s5Mi4Z398EkB6eOXTex/UxGyxBMX/fUUscBAABYnZkWDQcAAGBrmSkzB/PCeDwAAIZCMMe2s9plEQAAYAh0swQAABggwRwAAMAACeYAAAAGyJg5OAknGo9nchQAADaSYA7WiJkwAQDYSLpZAgAADJDMHGwQmTsAANaSzBwAAMAACeYAAAAGSDdL2CJW6oZ5IrpoAgDMH8EcbAOWTAAAmD+COdjmTLwCALA9GTMHAAAwQII5AACAAdLNEuacbpgAAMMkmAOeN4EgAMDmmSmYq6qLk7wlyY4kb+vuG6bOn5Lk7UlemeTxJJd399Gq+u4kPz5R9BuSvKK7P7AWlQfW32qWTAAAYP2sOGauqnYkuTnJJUnOTXJFVZ07VeyqJE9099lJbkpyY5J0969298u7++VJvifJUYEcAADA6s0yAcp5SY509yPd/WySO5PsmyqzL8nt4+27k1xYVTVV5ook71hNZQEAABiZpZvl6UkendhfTHL+cmW6+3hVPZnktCSfmihzeb4wCATmlPF2AACrM0tmbjrDliR9MmWq6vwkT3f3h5b8AVUHqupQVR365Cc/OUOVAAAA5tsswdxikjMn9s9I8thyZapqZ5JTkxybOL8/J+hi2d23dPfe7t67a9euWeoNAAAw12bpZnl/knOq6qwkH80oMHvDVJmFJFcmeW+SS5Pc292dJFX115JcluTVa1VpYBjMhAkAsH5WDObGY+CuSXJPRksT3Nbdh6vq+iSHunshya1J7qiqIxll5PZPfItXJ1ns7kfWvvoAAADzaaZ15rr7YJKDU8eum9h+JqPs21LX/vckr3r+VdwaZBgAAICtZJYxcwAAAGwxgjkAAIABmqmbJcBGsw4dAMCJCeaAQTpRsCfQAwDmgWAOmDsCQQBgOxDMAduO2WcBgHlgAhQAAIABEswBAAAMkG6WABPMogkADIXMHAAAwADJzAGsITNlAgAbRTAHcBLMlAkAbBW6WQIAAAyQzBzABjG5CgCwlmTmAAAABkgwBwAAMEC6WQJsEevZDdMsmwCw/cjMAQAADJBgDgAAYIBm6mZZVRcneUuSHUne1t03TJ0/Jcnbk7wyyeNJLu/uo+Nz35Dkl5P8zSSfS/KN3f3MWr0AgHlhjTsAYNKKmbmq2pHk5iSXJDk3yRVVde5UsauSPNHdZye5KcmN42t3JvmVJFd398uSXJDkL9es9gAAAHNqlszceUmOdPcjSVJVdybZl+TBiTL7krx5vH13krdWVSW5KMkfd/cfJUl3P75G9QZgjaw242cCFQDYHLOMmTs9yaMT+4vjY0uW6e7jSZ5MclqSr0nSVXVPVb2/qn5i9VUGAABglsxcLXGsZyyzM8k3J/nGJE8n+d2qeqC7f/fzLq46kORAkuzevXuGKgEAAMy3WTJzi0nOnNg/I8ljy5UZj5M7Ncmx8fH/0d2f6u6nkxxM8orpH9Ddt3T33u7eu2vXrpN/FQAAAHNmlszc/UnOqaqzknw0yf4kb5gqs5DkyiTvTXJpknu7u6vqniQ/UVUvTPJskn+Q0QQpAGwTFiQHgM2xYjDX3cer6pok92S0NMFt3X24qq5Pcqi7F5LcmuSOqjqSUUZu//jaJ6rqFzMKCDvJwe42tzYAAMAqzbTOXHcfzKiL5OSx6ya2n0ly2TLX/kpGyxMAAACwRmYZMwcAAMAWM1NmDgDWg/F2APD8CeYAWDerXZAcAFiebpYAAAADJDMHwJa0UlZPN0wA5p3MHAAAwAAJ5gAAAAZIMAcAADBAxswBMEjG1AEw72TmAAAABkhmDoBtyYLkAGx3gjkA5o4umgBsB7pZAgAADJDMHACcBFk9ALYKmTkAAIABkpkDgCkrZd8AYCuQmQMAABggmTkAWEPG1AGwUWTmAAAABmimYK6qLq6qh6vqSFVdu8T5U6rqrvH5+6pqz/j4nqr6TFV9YPz1S2tbfQAAgPm0YjfLqtqR5OYkr0mymOT+qlro7gcnil2V5InuPruq9ie5Mcnl43Mf7u6Xr3G9AQAA5tosY+bOS3Kkux9Jkqq6M8m+JJPB3L4kbx5v353krVVVa1hPANgWVjNTpvF2AEyaJZg7PcmjE/uLSc5frkx3H6+qJ5OcNj53VlX9YZJPJ/lX3f37q6syAMynEwWCAj2A+TNLMLdUhq1nLPOxJLu7+/GqemWS36iql3X3pz/v4qoDSQ4kye7du2eoEgAAwHybJZhbTHLmxP4ZSR5bpsxiVe1McmqSY93dSf4iSbr7gar6cJKvSXJo8uLuviXJLUmyd+/e6UARAFiBJREA5s8ss1nen+Scqjqrql6QZH+ShakyC0muHG9fmuTe7u6q2jWeQCVV9beSnJPkkbWpOgAAwPxaMTM3HgN3TZJ7kuxIclt3H66q65Mc6u6FJLcmuaOqjiQ5llHAlySvTnJ9VR1P8tkkV3f3sfV4IQAAAPNklm6W6e6DSQ5OHbtuYvuZJJctcd2vJ/n1VdYRAACAKTMFcwDAsBlTB7D9zDJmDgAAgC1GZg4AOCFZPYCtSWYOAABggGTmAIAVs28AbD0ycwAAAAMkMwcArIoxdQCbQ2YOAABggGTmAIB1daLMnawdwPMnMwcAADBAMnMAwJYlqwewPMEcALBpLIkA8PzpZgkAADBAgjkAAIAB0s0SABgk69sB804wBwBsS4I9YLvTzRIAAGCABHMAAAADJJgDAAAYoJnGzFXVxUnekmRHkrd19w1T509J8vYkr0zyeJLLu/voxPndSR5M8ubu/vm1qToAwPNnQXJg6FYM5qpqR5Kbk7wmyWKS+6tqobsfnCh2VZInuvvsqtqf5MYkl0+cvynJu9eu2gAA68fkKcAQzNLN8rwkR7r7ke5+NsmdSfZNldmX5Pbx9t1JLqyqSpKq+s4kjyQ5vDZVBgAAYJZulqcneXRifzHJ+cuV6e7jVfVkktOq6jNJ3phRVu/HVl9dAIDNp4smsBXMkpmrJY71jGX+dZKbuvupE/6AqgNVdaiqDn3yk5+coUoAAADzbZbM3GKSMyf2z0jy2DJlFqtqZ5JTkxzLKIN3aVX9XJIXJ/lcVT3T3W+dvLi7b0lyS5Ls3bt3OlAEABgM4+2AjTJLMHd/knOq6qwkH02yP8kbpsosJLkyyXuTXJrk3u7uJN/yXIGqenOSp6YDOQAAAE7eisHceAzcNUnuyWhpgtu6+3BVXZ/kUHcvJLk1yR1VdSSjjNz+9aw0AMBQydwBa2Wmdea6+2CSg1PHrpvYfibJZSt8jzc/j/oBAACwhFkmQAEAAGCLmSkzBwDAxrDsATArmTkAAIABEswBAAAMkGAOAABggARzAAAAAySYAwAAGCCzWQIADMRKC46fiJkwYfuRmQMAABggmTkAgDmwUlZP5g6GRzAHAIBgDwZIN0sAAIABkpkDAGBFJ8rcydrB5pCZAwAAGCDBHAAAwAAJ5gAAAAbImDkAAFbFTJiwOWTmAAAABkhmDgCAdWUmTFgfM2Xmquriqnq4qo5U1bVLnD+lqu4an7+vqvaMj59XVR8Yf/1RVb1+basPAAAwn1YM5qpqR5Kbk1yS5NwkV1TVuVPFrkryRHefneSmJDeOj38oyd7ufnmSi5P8clXJBgIAAKzSLJm585Ic6e5HuvvZJHcm2TdVZl+S28fbdye5sKqqu5/u7uPj41+cpNei0gAAAPNulmDu9CSPTuwvjo8tWWYcvD2Z5LQkqarzq+pwkg8muXoiuAMAAOB5mqXLYy1xbDrDtmyZ7r4vycuq6u8kub2q3t3dz3zexVUHkhxIkt27d89QJQAAtoOVljVYiQlUmGezZOYWk5w5sX9GkseWKzMeE3dqkmOTBbr7oSR/nuTrpn9Ad9/S3Xu7e++uXbtmrz0AAMCcmiUzd3+Sc6rqrCQfTbI/yRumyiwkuTLJe5NcmuTe7u7xNY929/Gq+uokX5vk6FpVHgCA+WbZA+bZisHcOBC7Jsk9SXYkua27D1fV9UkOdfdCkluT3FFVRzLKyO0fX/7NSa6tqr9M8rkkP9Tdn1qPFwIAADBPZlomoLsPJjk4dey6ie1nkly2xHV3JLljlXUEAABgijXfAADYllaaXEU3TIZulglQAAAA2GJk5gAAmEsydwydzBwAAMAACeYAAAAGSDdLAABYgjXs2Opk5gAAAAZIMAcAADBAgjkAAIABEswBAAAMkAlQAADgJFmjjq1AMAcAAGvMTJhsBN0sAQAABkgwBwAAMECCOQAAgAEyZg4AALYQk6swK5k5AACAAZKZAwCADbRS5g1mJTMHAAAwQII5AACAAZqpm2VVXZzkLUl2JHlbd98wdf6UJG9P8sokjye5vLuPVtVrktyQ5AVJnk3y49197xrWHwAA5ooFyXnOipm5qtqR5OYklyQ5N8kVVXXuVLGrkjzR3WcnuSnJjePjn0ryj7r765NcmeSOtao4AADAPJulm+V5SY509yPd/WySO5PsmyqzL8nt4+27k1xYVdXdf9jdj42PH07yxeMsHgAAAKswSzfL05M8OrG/mOT85cp09/GqejLJaRll5p7zj5P8YXf/xfOvLgAAsBxr1M2XWYK5WuJYn0yZqnpZRl0vL1ryB1QdSHIgSXbv3j1DlQAAAObbLN0sF5OcObF/RpLHlitTVTuTnJrk2Hj/jCTvSvK93f3hpX5Ad9/S3Xu7e++uXbtO7hUAAADMoVmCufuTnFNVZ1XVC5LsT7IwVWYhowlOkuTSJPd2d1fVi5P8VpI3dfcfrFWlAQAA5t2KwVx3H09yTZJ7kjyU5J3dfbiqrq+q142L3ZrktKo6kuRHk1w7Pn5NkrOT/FRVfWD89WVr/ioAAADmTHVPD3/bXHv37u1Dhw5tdjW+wEqDSQEAYMhMjrJ1VNUD3b13pXKzdLMEAABgi5llNksAAGCbs6zB8MjMAQAADJBgDgAAYIB0swQAAFakG+bWIzMHAAAwQII5AACAARLMAQAADJBgDgAAYIBMgAIAAKzaiSZIMTnK+pCZAwAAGCDBHAAAwAAJ5gAAAAbImDkAAGBdWXB8fcjMAQAADJBgDgAAYIAEcwAAAANkzBwAALCprFH3/MjMAQAADNBMwVxVXVxVD1fVkaq6donzp1TVXePz91XVnvHx06rq96rqqap669pWHQAAYH6tGMxV1Y4kNye5JMm5Sa6oqnOnil2V5InuPjvJTUluHB9/JslPJfmxNasxAAAAM2XmzktypLsf6e5nk9yZZN9UmX1Jbh9v353kwqqq7v7z7v5fGQV1AAAArJFZgrnTkzw6sb84PrZkme4+nuTJJKetRQUBAAD4QrPMZllLHOvnUWb5H1B1IMmBJNm9e/eslwEAANvciWa6TOZ7tstZMnOLSc6c2D8jyWPLlamqnUlOTXJs1kp09y3dvbe79+7atWvWywAAAObWLMHc/UnOqaqzquoFSfYnWZgqs5DkyvH2pUnu7e6ZM3MAAACcnBW7WXb38aq6Jsk9SXYkua27D1fV9UkOdfdCkluT3FFVRzLKyO1/7vqqOprkbyZ5QVV9Z5KLuvvBtX8pAAAA82OWMXPp7oNJDk4du25i+5kkly1z7Z5V1A8AAIAlzLRoOAAAAFvLTJk5AACArWieZ7uUmQMAABggwRwAAMAACeYAAAAGSDAHAAAwQII5AACAARLMAQAADJBgDgAAYIAEcwAAAAMkmAMAABignZtdAQAAgPWy59rfWvbc0Ru+fQNrsvZk5gAAAAZIMAcAADBAgjkAAIABEswBAAAMkGAOAABggARzAAAAAySYAwAAGKCZgrmquriqHq6qI1V17RLnT6mqu8bn76uqPRPn3jQ+/nBVfdvaVR0AAGB+rRjMVdWOJDcnuSTJuUmuqKpzp4pdleSJ7j47yU1Jbhxfe26S/UleluTiJP9h/P0AAABYhVkyc+clOdLdj3T3s0nuTLJvqsy+JLePt+9OcmFV1fj4nd39F939f5McGX8/AAAAVmGWYO70JI9O7C+Ojy1ZpruPJ3kyyWkzXgsAAMBJ2jlDmVriWM9YZpZrU1UHkhwY7z5VVQ/PUK+N9NIkn9rsSpBEW2wl2mLr0BZbh7bYOrTF1qEttg5tMaVu3LQfvVJbfPUs32SWYG4xyZkT+2ckeWyZMotVtTPJqUmOzXhtuvuWJLfMUuHNUFWHunvvZtcDbbGVaIutQ1tsHdpi69AWW4e22Dq0xdaxVm0xSzfL+5OcU1VnVdULMprQZGGqzEKSK8fblya5t7t7fHz/eLbLs5Kck+R/r7bSAAAA827FzFx3H6+qa5Lck2RHktu6+3BVXZ/kUHcvJLk1yR1VdSSjjNz+8bWHq+qdSR5McjzJD3f3Z9fptQAAAMyNWbpZprsPJjk4dey6ie1nkly2zLU/m+RnV1HHrWDLdgGdQ9pi69AWW4e22Dq0xdahLbYObbF1aIutY03aoka9IQEAABiSWcbMAQAAsMUI5saq6uKqeriqjlTVtUucP6Wq7hqfv6+q9mx8Lbe/qjqzqn6vqh6qqsNV9S+WKHNBVT1ZVR8Yf1231PdibVTV0ar64Ph3fWiJ81VV/258b/xxVb1iM+q53VXV1078n/9AVX26qn5kqox7Y51U1W1V9Ymq+tDEsS+tqvdU1Z+O/33JMtdeOS7zp1V15VJlmN0ybfFvqupPxu9B76qqFy9z7Qnfzzg5y7TFm6vqoxPvQ69d5toTPndxcpZpi7sm2uFoVX1gmWvdF2tkuefY9fy80M0ySVXtSPJ/krwmo+UU7k9yRXc/OFHmh5J8Q3dfXVX7k7y+uy/flApvY1X1lUm+srvfX1V/I8kDSb5zqi0uSPJj3f0dm1TNuVJVR5Ps7e4l10IZf1D/8ySvTXJ+krd09/kbV8P5M37P+miS87v7zyaOXxD3xrqoqlcneSrJ27v768bHfi7Jse6+Yfww+pLufuPUdV+a5FCSvRmts/pAkld29xMb+gK2kWXa4qKMZtI+XjVaNWq6LcbljuYE72ecnGXa4s1Jnurunz/BdSs+d3FylmqLqfO/kOTJ7r5+iXNH475YE8s9xyb5p1mnzwuZuZHzkhzp7ke6+9kkdybZN1VmX5Lbx9t3J7mwqpZaFJ1V6O6Pdff7x9v/L8lDSU7f3Fqxgn0ZfXh0d78vyYvHb2asnwuTfHgykGN9dff/zGi25kmTnwu3Z/SBPe3bkrynu4+NP5Dfk+TidavoHFiqLbr7d7r7+Hj3fRmta8s6W+a+mMUsz12chBO1xfh59Z8keceGVmoOneA5dt0+LwRzI6cneXRifzFfGED8VZnxB8aTSU7bkNrNqRp1Zf17Se5b4vTfr6o/qqp3V9XLNrRi86eT/E5VPVBVB5Y4P8v9w9ran+U/lN0bG+fLu/tjyegDPMmXLVHG/bHx/lmSdy9zbqX3M9bGNeMur7ct053MfbGxviXJx7v7T5c5775YB1PPsev2eSGYG1kqwzbd/3SWMqyRqvqSJL+e5Ee6+9NTp9+f5Ku7++8m+fdJfmOj6zdnvqm7X5HkkiQ/PO7KMcm9sYGq6gVJXpfkvyxx2r2x9bg/NlBV/cuM1rX91WWKrPR+xur9xyR/O8nLk3wsyS8sUcZ9sbGuyImzcu6LNbbCc+yyly1xbMX7QjA3spjkzIn9M5I8tlyZqtqZ5NQ8v64FrKCqviijG+BXu/u/Tp/v7k9391Pj7YNJvqiqXrrB1Zwb3f3Y+N9PJHlXRt1jJs1y/7B2Lkny/u7++PQJ98aG+/hzXYrH/35iiTLujw0ynizgO5J8dy8zIcAM72esUnd/vLs/292fS/KfsvTv2H2xQcbPrN+V5K7lyrgv1tYyz7Hr9nkhmBu5P8k5VXXW+K/e+5MsTJVZSPLcrDKXZjTQ2l+R1ti4X/etSR7q7l9cpsxXPDdesarOy+j/8eMbV8v5UVUvGg/gTVW9KMlFST40VWwhyffWyKsyGmD9sQ2u6jxZ9i+s7o0NN/m5cGWS31yizD1JLqqql4y7m100PsYaqqqLk7wxyeu6++llyszyfsYqTY2Zfn2W/h3P8tzF2vjWJH/S3YtLnXRfrK0TPMeu2+fFztVVeXsYz351TUa/sB1Jbuvuw1V1fZJD3b2QUcPcUVVHMsrI7d+8Gm9r35Tke5J8cGIK3Z9MsjtJuvuXMgqmf7Cqjif5TJL9Aut18+VJ3jWOD3Ym+bXu/u2qujr5q/Y4mNFMlkeSPJ3k+zaprtteVb0wo9nffmDi2GRbuDfWSVW9I8kFSV5aVYtJfjrJDUneWVVXJflIksvGZfcmubq7v7+7j1XVz2T08Jok13e3Xh2rsExbvCnJKUneM36/et949umvSvK27n5tlnk/24SXsG0s0xYXVNXLM+oedjTj96vJtljuuWsTXsK2sVRbdPetWWKMtftiXS33HLtunxeWJgAAABgg3SwBAAAGSDAHAAAwQII5AACAARLMAQAADJBgDgAAYIAEcwAAAAMkmAMAABggwRwAAMAA/X9kg84FMo7+RwAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 1080x720 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"val=dphi2_np(np.sqrt(q2_range),mmu_mass,mmu_mass)*dphi2_np(mB_mass,mKst_mass,np.sqrt(q2_range))\n",
"plt.subplot(2,1,1)\n",
"plt.plot(q2_range/1e6, val);\n",
"fig=plt.gcf()\n",
"plt.subplot(2,1,2)\n",
"plt.hist(q2_array/1e6, weights=weights_np, bins=100,density=True);\n",
"fig=plt.gcf()\n",
"fig.set_size_inches(15,10);"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"class GaussianSampleAndWeights():\n",
"\n",
" def __call__(self, n_to_produce, limits, dtype):\n",
" \n",
" n_to_produce = tf.cast(n_to_produce, dtype=tf.int32)\n",
"\n",
" weights, phase_space_sample = B.generate_tensor(n_events=n_to_produce) #check order \n",
" \n",
" weights_out = 1/weights\n",
" phase_space_sample_tensor = tf.concat([phase_space_sample[\"Kst\"],phase_space_sample[\"q\"],phase_space_sample[\"l1\"],phase_space_sample[\"l2\"],], axis=1)\n",
" \n",
" weights_max = None\n",
" \n",
" thresholds = tf.random_uniform(shape=(n_to_produce,), dtype=dtype)\n",
" \n",
" return phase_space_sample_tensor, thresholds, weights_out, weights_max, n_to_produce\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"class dGamma(zfit.pdf.ZPDF):\n",
" \n",
" _PARAMS = ['mB','mKst','ml']\n",
" _N_OBS = 16\n",
"\n",
" def _unnormalized_pdf(self, x):\n",
" \n",
" ml = self.params['ml']\n",
" mB = self.params['mB']\n",
" mKst = self.params['mKst']\n",
" \n",
" list_ps= x.unstack_x()\n",
" new_list_ps = [tf.expand_dims(ele, axis=1) for ele in list_ps] \n",
" pKstx, pKsty, pKstz, pKstE, qx, qy, qz, qE, p1x, p1y, p1z, p1E, p2x, p2y, p2z, p2E = new_list_ps\n",
" \n",
" pKst = lorentz_vector(tf.concat([pKstx, pKsty, pKstz],axis=-1), pKstE)\n",
" p1 = lorentz_vector(tf.concat([p1x, p1y, p1z],axis=-1), p1E)\n",
" p2 = lorentz_vector(tf.concat([p2x, p2y, p2z],axis=-1), p2E)\n",
" qvec_B = tf.concat([qx, qy, qz],axis=-1)\n",
" q = lorentz_vector(qvec_B, qE)\n",
" \n",
" z = tf.cast(tf.expand_dims(tf.stack([0., 0., 1., 0.], axis=0), axis=0),dtype=tf.float64)\n",
" \n",
" #assert q==p1+p2\n",
" q2=scalar_product_4(q,q)\n",
" modq = ztf.sqrt(q2)\n",
" \n",
" #q_in_q_rf = lorentz_vector(tf.zeros_like(qvec_B), modq)\n",
" p1q = lorentz_boost(p1, qvec_B/qE)\n",
" \n",
" costheta_l= tf.expand_dims(get_costheta_l(p1q, z), axis=1)\n",
" \n",
" \n",
" pdf = matrix_elt(q2, ml, mB, mKst, costheta_l)\n",
" #pdf = tf.ones_like(q2, dtype=tf.float64)\n",
" print(pdf)\n",
" return pdf[:, 0]\n"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"lepton_min_E = 105\n",
"lepton_max_E = 2700\n",
"\n",
"lepton_min_p = -2700\n",
"lepton_max_p = 2700\n",
"\n",
"Kst_min_E = 892.\n",
"Kst_max_E = 2700\n",
"\n",
"Kst_min_p = -2700\n",
"Kst_max_p = 2700\n",
"\n",
"q_min_E = 2700\n",
"q_max_E = 4400\n",
"\n",
"q_min_p = -2700\n",
"q_max_p = 2700"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"lower =((Kst_min_p, Kst_min_p, Kst_min_p, mKst_mass, q_min_p, q_min_p, q_min_p, q_min_E, lepton_min_p, lepton_min_p, lepton_min_p, lepton_min_E, lepton_min_p, lepton_min_p, lepton_min_p, lepton_min_E,),)\n",
"upper =((Kst_max_p , Kst_max_p , Kst_max_p , Kst_max_E, q_max_p, q_max_p, q_max_p, q_max_E, lepton_max_p, lepton_max_p , lepton_max_p , lepton_max_E, lepton_max_p , lepton_max_p , lepton_max_p , lepton_max_E,),)\n",
"\n",
"obs = zfit.Space([\"pKstx\",\"pKsty\",\"pKstz\",\"pKstE\", \n",
" \"qx\", \"qy\", \"qz\", \"qE\",\n",
" \"p1x\" ,\"p1y\" ,\"p1z\" ,\"p1E\" , \n",
" \"p2x\" ,\"p2y\" ,\"p2z\" ,\"p2E\" ],\n",
" limits=(lower,upper))"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"mlepton_par = zfit.Parameter('ml', mmu_mass)\n",
"mB_par = zfit.Parameter('mB', mB_mass)\n",
"mKst_par = zfit.Parameter('mKst', mKst_mass)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"pdf = dGamma(obs=obs, ml = mlepton_par, mB = mB_par, mKst = mKst_par)\n",
"#pdf = dGamma(obs=obs, ml = m_lepton)\n",
"pdf._sample_and_weights=GaussianSampleAndWeights"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"WARNING:tensorflow:From /Users/davide/miniconda3/envs/zfit_env/lib/python3.6/site-packages/zfit/core/sample.py:163: to_int32 (from tensorflow.python.ops.math_ops) is deprecated and will be removed in a future version.\n",
"Instructions for updating:\n",
"Use tf.cast instead.\n",
"Tensor(\"create_sampler/while/unnormalized_pdf/add_21:0\", shape=(?, 1), dtype=float64)\n"
]
}
],
"source": [
"sampler = pdf.create_sampler(n=n_events)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"zfit.settings.set_verbosity(10)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"for i in range(1):\n",
" sampler.resample()\n"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA3kAAAEvCAYAAAD4uAgWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3de7Bd5Xnf8e+vgIlTO0GYSxVdKtIoF+KZ2OQUaN1JqUmEwJ6IzIQEu2NkwozSDG7sOm0tkszg+tLBbWLHnrh4lKJYZBwwwc6gceQQGcN4PBMwghBsLLsomMIxKlIijO3SkOI+/WO/J96S9j7onLPPZa/9/czs2Ws/613rrPXqHL3z7PeyUlVIkiRJkrrhHyz3BUiSJEmSRsckT5IkSZI6xCRPkiRJkjrEJE+SJEmSOsQkT5IkSZI6xCRPkiRJkjrk5OW+gPk644wzasOGDct9GZKkRfbAAw/8dVWdudzXMS5sHyVpcgxrI8c2yduwYQP79u1b7suQJC2yJP9zua9hnNg+StLkGNZGOlxTkiRJkjrEJE+SJEmSOsQkT5IkSZI6xCRPkiRJkjrEJE+SJEmSOsQkT5IkSZI6xCRPkiRJkjrEJE+SJEmSOsQkT5IkSZI6xCRPkiRJkjrEJE+SJEmSOuTk5b4ASSduw/Y/GRh//IbXLfGVSJK0PAa1hbaD0tFM8qQVaFgyJ0mSJL0YkzxpGZnMSZIkadSckydJkiRJHfKiPXlJdgKvBw5V1Stb7HTg48AG4HHgF6rqmSQBPghcBjwHvLmqHmzHbAV+s532PVW1q8V/Evgo8FJgD/DWqqoR3Z8kSZI6zjnr0tFOZLjmR4HfBW7ui20H7qqqG5Jsb5/fAVwKbGyvC4AbgQtaUng9MAUU8ECS3VX1TCuzDbiXXpK3Gfj0wm9NmhxOQpckSdKMFx2uWVWfA44cE94C7Grbu4DL++I3V8+9wGlJVgOXAHur6khL7PYCm9u+76uqP2+9dzf3nUuSJEmSNEfznZN3dlUdBGjvZ7X4GuDJvnLTLTZbfHpAXJIkSZI0D6NeeCUDYjWP+OCTJ9uS7Euy7/Dhw/O8REmSJEnqrvk+QuHpJKur6mAbcnmoxaeBdX3l1gJPtfhFx8TvafG1A8oPVFU7gB0AU1NTLs4iSZLUUT5mSJq/+fbk7Qa2tu2twB198avScyHwbBvOeSewKcmqJKuATcCdbd+3klzYVua8qu9ckiRJkqQ5OpFHKNxCrxfujCTT9FbJvAG4Lck1wBPAFa34HnqPTzhA7xEKVwNU1ZEk7wbub+XeVVUzi7n8Ct99hMKncWVNSZIkSZq3F03yquoNQ3ZdPKBsAdcOOc9OYOeA+D7glS92HZLmxmcGSZIkTab5zsmTNEfOLZAmR5LvAT4HnEqvrb29qq5P8lHgXwLPtqJvrqqH2pSFD9IbDfNciz/YzrUV+M1W/j1VtQtJkmZhkidJ0ug9D7y2qr6d5BTg80lmpiP8h6q6/ZjylwIb2+sC4EbggiSn05smMUVv9ekHkuxuz5yVJGmgUT9CQZKkiVc9324fT2mv2VaF3gLc3I67FzitrV59CbC3qo60xG4vsHkxr12SNP5M8iRJWgRJTkryEL3HDO2tqvvarvcmeTjJB5Kc2mJrgCf7Dp9usWFxSZKGMsmTJGkRVNV3qupV9J4Be36SVwLXAT8K/FPgdOAdrXgGnWKW+FGSbEuyL8m+w4cPj+T6JUnjyzl5kiQtoqr6RpJ7gM1V9Vst/HyS3wf+ffs8DazrO2wt8FSLX3RM/J4BP2MHsANgampqtmGh0kRxpWlNKnvyJEkasSRnJjmtbb8U+GngK22eHW01zcuBL7VDdgNXpedC4NmqOgjcCWxKsirJKmBTi0mSNJQ9edKI+agEScBqYFeSk+h9oXpbVX0qyWeTnElvGOZDwL9p5ffQe3zCAXqPULgaoKqOJHk3cH8r966qOrKE9yFJGkMmeZIkjVhVPQy8ekD8tUPKF3DtkH07gZ0jvUBJUqc5XFOSJEmSOsSePEmSJC0bpzlIo2dPniRJkiR1iEmeJEmSJHWIwzWlCeMzgyRJk862UF1nT54kSZIkdYhJniRJkiR1iMM1pXlyNTBJkiStRPbkSZIkSVKH2JMnCXASuiRJUlfYkydJkiRJHWKSJ0mSJEkdYpInSZIkSR3inDxJkiQtOlellpaOPXmSJEmS1CEmeZIkSZLUISZ5kiRJktQhzsmTJEmSGDxv0OfFahzZkydJkiRJHWJPnvQiJn01ML/VlCRJGi/25EmSJElSh5jkSZIkSVKHmORJkiRJUoeY5EmSJElSh5jkSZIkSVKHmORJkjRiSb4nyReS/GWSR5L8pxY/J8l9SR5N8vEkL2nxU9vnA23/hr5zXdfiX01yyfLckSRpnPgIBamZ9EclSBqp54HXVtW3k5wCfD7Jp4G3Ax+oqluTfAS4BrixvT9TVT+U5ErgfcAvJjkXuBL4ceAHgM8k+eGq+s5y3JR0orrUpg67Fx8npJXMJE/SnNngSbOrqgK+3T6e0l4FvBZ4Y4vvAt5JL8nb0rYBbgd+N0la/Naqeh74WpIDwPnAny/+XUiSxpXDNSVJWgRJTkryEHAI2Av8FfCNqnqhFZkG1rTtNcCTAG3/s8Ar+uMDjpEkaSCTPEmSFkFVfaeqXgWspdf79mODirX3DNk3LH6UJNuS7Euy7/Dhw/O9ZElSR5jkSZK0iKrqG8A9wIXAaUlmpkqsBZ5q29PAOoC2//uBI/3xAcf0/4wdVTVVVVNnnnnmYtyGJGmMmORJkjRiSc5Mclrbfinw08B+4G7g51uxrcAdbXt3+0zb/9k2r283cGVbffMcYCPwhaW5C0nSuFpQkpfk37Wlob+U5Ja2ZLTLQ0uSJt1q4O4kDwP3A3ur6lPAO4C3twVUXgHc1MrfBLyixd8ObAeoqkeA24AvA38KXOvKmpKkFzPv1TWTrAF+FTi3qv5PktvoLfN8GS4PLUmaYFX1MPDqAfHH6M3POzb+t8AVQ871XuC9o75GSVJ3LXS45snAS9v8ge8FDtJbHvr2tn8XcHnb3tI+0/ZffOzy0FX1NWBmeWhJkiRJ0hzNO8mrqq8DvwU8QS+5exZ4AJeHliRJkqRls5Dhmqvo9cKdA3wD+CPg0gFFR7I8dPuZ24BtAOvXr5/jFUvfNexh3pIkSdK4m3eSR2+lsK9V1WGAJJ8E/jlteejWWzdoeejp+SwPDb0looEdAFNTUwMTQUnLZ1jy/PgNr1viK5EkaXHZ5mklW8icvCeAC5N8b5tbdzG91b9cHlqSJEmSlsm8e/Kq6r4ktwMPAi8Af0Gvl+1PgFuTvKfF+peH/oO2PPQReitqUlWPtJU5v9zO4/LQkiRJY8DpD9LKtJDhmlTV9cD1x4RdHlqSJEmSlslCH6EgSZIkSVpBTPIkSZIkqUNM8iRJkiSpQxY0J0+SJEnd5wIr0nixJ0+SJEmSOsSePEmSJGlEfEi6VgJ78iRJkiSpQ+zJU6c5h0CSJEmTxiRP0qJz6IokSdLScbimJEmSJHWISZ4kSZIkdYhJniRJkiR1iEmeJEmSJHWISZ4kSZIkdYhJniRJkiR1iEmeJEmSJHWIz8lTJ/jQ8/E06N/NZ+dJ0vKxPZW6wSRPkiRJWmR+saml5HBNSZIkSeoQkzxJkkYsybokdyfZn+SRJG9t8Xcm+XqSh9rrsr5jrktyIMlXk1zSF9/cYgeSbF+O+5EkjReHa0qSNHovAL9WVQ8meTnwQJK9bd8Hquq3+gsnORe4Evhx4AeAzyT54bb7w8DPANPA/Ul2V9WXl+QuJEljySRPkqQRq6qDwMG2/a0k+4E1sxyyBbi1qp4HvpbkAHB+23egqh4DSHJrK2uSJ0kayuGakiQtoiQbgFcD97XQW5I8nGRnklUttgZ4su+w6RYbFpckaSiTPEmSFkmSlwGfAN5WVd8EbgT+CfAqej19vz1TdMDhNUv82J+zLcm+JPsOHz48kmuXJI0vh2tKkrQIkpxCL8H7WFV9EqCqnu7b/3vAp9rHaWBd3+Frgafa9rD436uqHcAOgKmpqeOSQEkr07DnEvpoBS2USZ6kFcUGT12QJMBNwP6qen9ffHWbrwfwc8CX2vZu4A+TvJ/ewisbgS/Q68nbmOQc4Ov0Fmd549LchSRpXJnkSZI0eq8B3gR8MclDLfbrwBuSvIrekMvHgV8GqKpHktxGb0GVF4Brq+o7AEneAtwJnATsrKpHlvJGJEnjxyRPY2VYL48krSRV9XkGz6fbM8sx7wXeOyC+Z7bjJEk6lguvSJIkSVKHmORJkiRJUoeY5EmSJElShzgnT5IkacI4x13qNnvyJEmSJKlD7MmTNBZ8fp4kSdKJsSdPkiRJkjrEJE+SJEmSOsThmlqxnBQuSZIkzZ09eZIkSZLUISZ5kiRJktQhDteUJEmSVhBXlNZC2ZMnSZIkSR2yoCQvyWlJbk/ylST7k/yzJKcn2Zvk0fa+qpVNkg8lOZDk4STn9Z1nayv/aJKtC70pSZIkSZpUC+3J+yDwp1X1o8BPAPuB7cBdVbURuKt9BrgU2Nhe24AbAZKcDlwPXACcD1w/kxhKkiRJkuZm3kleku8Dfgq4CaCq/q6qvgFsAXa1YruAy9v2FuDm6rkXOC3JauASYG9VHamqZ4C9wOb5XpckSZIkTbKFLLzyg8Bh4PeT/ATwAPBW4OyqOghQVQeTnNXKrwGe7Dt+usWGxSXpRTk5XZJ6fL6spBkLGa55MnAecGNVvRr433x3aOYgGRCrWeLHnyDZlmRfkn2HDx+e6/VKkiRJUuctpCdvGpiuqvva59vpJXlPJ1ndevFWA4f6yq/rO34t8FSLX3RM/J5BP7CqdgA7AKampgYmgho/fvMoSZIkjc68e/Kq6n8BTyb5kRa6GPgysBuYWSFzK3BH294NXNVW2bwQeLYN67wT2JRkVVtwZVOLSZIkSZLmaKEPQ/+3wMeSvAR4DLiaXuJ4W5JrgCeAK1rZPcBlwAHguVaWqjqS5N3A/a3cu6rqyAKvS5IkSZIm0oKSvKp6CJgasOviAWULuHbIeXYCOxdyLZIkSZKkhT8nT5IkSZK0gpjkSZIkSVKHLHROniStSD4/T5IkTSqTPEmSJGkM+AWmTpTDNSVJGrEk65LcnWR/kkeSvLXFT0+yN8mj7X1ViyfJh5IcSPJwkvP6zrW1lX80ydZhP1OSpBn25GnJ+NBzSRPkBeDXqurBJC8HHkiyF3gzcFdV3ZBkO7AdeAdwKbCxvS4AbgQuSHI6cD29layrnWd3VT2z5HckSRobJnmSJI1YVR0EDrbtbyXZD6wBtgAXtWK7gHvoJXlbgJvb44buTXJaktWt7N6Z58e2RHEzcMuS3YykFc9hnDqWSZ4kSYsoyQbg1cB9wNktAaSqDiY5qxVbAzzZd9h0iw2LH/sztgHbANavXz/aG9CK48gYSS/GOXmSJC2SJC8DPgG8raq+OVvRAbGaJX50oGpHVU1V1dSZZ545v4uVJHWGSZ4kSYsgySn0EryPVdUnW/jpNgyT9n6oxaeBdX2HrwWemiUuSdJQJnmSJI1YkgA3Afur6v19u3YDMytkbgXu6Itf1VbZvBB4tg3rvBPYlGRVW4lzU4tJkjSUc/IkSRq91wBvAr6Y5KEW+3XgBuC2JNcATwBXtH17gMuAA8BzwNUAVXUkybuB+1u5d80swiJJ0jAmeZIkjVhVfZ7B8+kALh5QvoBrh5xrJ7BzdFenceECK5Lmy+GakiRJktQh9uRp5PzmUSvZoN9PnyMkSZK6xJ48SZIkSeoQe/IkSZKkDnL0yuSyJ0+SJEmSOsQkT5IkSZI6xCRPkiRJkjrEJE+SJEmSOsQkT5IkSZI6xNU1tSA+E0+SJElaWUzyJE28YV9WuMy0JEkaRw7XlCRJkqQOMcmTJEmSpA5xuKYkSZI0IZyiMBnsyZMkSZKkDrEnT5IkaZm5WrWkUbInT5IkSZI6xCRPkiRJkjrEJE+SJEmSOsQ5eTohzhWQJGnhbE8lLQWTPEkawmWmJUnSOHK4piRJkiR1iEmeJEmSJHWIwzUlSZKkCecUhW6xJ0+SJEmSOsQkT5KkEUuyM8mhJF/qi70zydeTPNRel/Xtuy7JgSRfTXJJX3xzix1Isn2p70OSNJ5M8iRJGr2PApsHxD9QVa9qrz0ASc4FrgR+vB3z35KclOQk4MPApcC5wBtaWUmSZuWcPEmSRqyqPpdkwwkW3wLcWlXPA19LcgA4v+07UFWPASS5tZX98ogvV5LUMfbkSZK0dN6S5OE2nHNVi60BnuwrM91iw+KSJM1qwT15bTjJPuDrVfX6JOcAtwKnAw8Cb6qqv0tyKnAz8JPA3wC/WFWPt3NcB1wDfAf41aq6c6HXpfkZtrKSJGnBbgTeDVR7/23gl4AMKFsM/iK2Bp04yTZgG8D69etHca2SBLjq5rgaxXDNtwL7ge9rn99Hb87BrUk+Qi95u7G9P1NVP5TkylbuF4+Zi/ADwGeS/HBVfWcE1yZJI2eDp/moqqdntpP8HvCp9nEaWNdXdC3wVNseFj/23DuAHQBTU1MDE0FJ0uRYUJKXZC3wOuC9wNuTBHgt8MZWZBfwTnpJ3pa2DXA78Lut/LC5CH++kGuTJGklSbK6qg62jz8HzKy8uRv4wyTvp/dl50bgC/R6+Da2ETJfp/eF6BvRWHBkjKTltNCevN8B/iPw8vb5FcA3quqF9rl//sDfzy2oqheSPNvKrwHu7Tvn0DkHDkeRJI2DJLcAFwFnJJkGrgcuSvIqekMuHwd+GaCqHklyG70FVV4Arp0ZzZLkLcCdwEnAzqp6ZIlvRZI0huad5CV5PXCoqh5IctFMeEDRepF9sx1zdNDhKJKkMVBVbxgQvmmW8u+lNyrm2PgeYM8IL02SNAEW0pP3GuBn28Ncv4fenLzfAU5LcnLrzeufPzAz52A6ycnA9wNHmH0ugiRJkiRpDub9CIWquq6q1lbVBnrzBD5bVf8auBv4+VZsK3BH297dPtP2f7aqqsWvTHJqm3cwMxdBkiRJkjRHi/Ew9HcAtyZ5D/AXfHd4yk3AH7SFVY7QSwxnnYsgSZIkSZqbkSR5VXUPcE/bfoze6pjHlvlb4Iohxw+ciyBJkiRJmpvF6MnTGHBpZ0mSJKmbTPIkaUQGfXniA9IlSdJSM8mTJEmSNCd+sbmyzXt1TUmSJEnSymNPniRJ0jw5x13SSmRPniRJkiR1iEmeJEmSJHWISZ4kSZIkdYhz8iaA8wWk5TPs788VyCRJXWObt3LYkydJkiRJHWJPniRJ0glwZIykcWFPniRJkiR1iEmeJEmSJHWISZ4kSZIkdYhJniRJkiR1iEmeJEmSJHWISZ4kSZIkdYhJniRJkiR1iM/J6xCf3yONj2F/r4/f8LolvhJJkhaXbd7SsydPkiRJkjrEJE+SJEmSOsQkT5IkSZI6xDl5kiSNWJKdwOuBQ1X1yhY7Hfg4sAF4HPiFqnomSYAPApcBzwFvrqoH2zFbgd9sp31PVe1ayvuYVM5xlzTu7MmTJGn0PgpsPia2HbirqjYCd7XPAJcCG9trG3Aj/H1SeD1wAXA+cH2SVYt+5ZKksWdPniRJI1ZVn0uy4ZjwFuCitr0LuAd4R4vfXFUF3JvktCSrW9m9VXUEIMleeonjLYt8+ZK0JFx1c/HYkydJ0tI4u6oOArT3s1p8DfBkX7npFhsWlyRpViZ5kiQtrwyI1Szx40+QbEuyL8m+w4cPj/TiJEnjx+GakiQtjaeTrK6qg2045qEWnwbW9ZVbCzzV4hcdE79n0ImragewA2BqampgIqjjucCKpK4yyRtDNkqSNJZ2A1uBG9r7HX3xtyS5ld4iK8+2RPBO4D/3LbayCbhuia9ZkjSGTPIkSRqxJLfQ64U7I8k0vVUybwBuS3IN8ARwRSu+h97jEw7Qe4TC1QBVdSTJu4H7W7l3zSzCIknSbEzyJGkFcaWxbqiqNwzZdfGAsgVcO+Q8O4GdI7w0SdIEcOEVSZIkSeoQe/IkSVKnOZdd0qQxyZMkSZK0Ygz6YsZpC3PjcE1JkiRJ6hCTPEmSJEnqEJM8SZIkSeoQkzxJkiRJ6hAXXlnBXA1M0gwnoUuSpBNlT54kSZIkdYg9eZIkqTMcBSN107C/bUe1DDbvnrwk65LcnWR/kkeSvLXFT0+yN8mj7X1ViyfJh5IcSPJwkvP6zrW1lX80ydaF35YkSZIkTaaFDNd8Afi1qvox4ELg2iTnAtuBu6pqI3BX+wxwKbCxvbYBN0IvKQSuBy4Azgeun0kMJUmSJElzM+8kr6oOVtWDbftbwH5gDbAF2NWK7QIub9tbgJur517gtCSrgUuAvVV1pKqeAfYCm+d7XZIkSZI0yUay8EqSDcCrgfuAs6vqIPQSQeCsVmwN8GTfYdMtNiwuSZIkSZqjBS+8kuRlwCeAt1XVN5MMLTogVrPEB/2sbfSGerJ+/fq5X6wkSZKkznBBlsEWlOQlOYVegvexqvpkCz+dZHVVHWzDMQ+1+DSwru/wtcBTLX7RMfF7Bv28qtoB7ACYmpoamAiOI1cCkzQfNmySJGmQhayuGeAmYH9Vvb9v125gZoXMrcAdffGr2iqbFwLPtuGcdwKbkqxqC65sajFJkiRJ0hwtpCfvNcCbgC8meajFfh24AbgtyTXAE8AVbd8e4DLgAPAccDVAVR1J8m7g/lbuXVV1ZAHXJUmSOs5RMJI03LyTvKr6PIPn0wFcPKB8AdcOOddOYOd8r0WSJEmS1DOS1TUlSZIkSSvDglfXlCRJWiwOy5Q0H5O+OJk9eZIkSZLUISZ5kiRJktQhJnmSJEmS1CHOyZOkjpn0eQiSJE06k7wl5gRySZpsSR4HvgV8B3ihqqaSnA58HNgAPA78QlU9kyTAB+k9Z/Y54M1V9eByXLckaXw4XFOSpKX3r6rqVVU11T5vB+6qqo3AXe0zwKXAxvbaBty45FcqSRo7JnmSJC2/LcCutr0LuLwvfnP13AuclmT1clygJGl8mORJkrS0CvizJA8k2dZiZ1fVQYD2flaLrwGe7Dt2usUkSRrKOXmSJC2t11TVU0nOAvYm+cosZTMgVscV6iWL2wDWr18/mqtcYs5Zl6TRMcmTJGkJVdVT7f1Qkj8GzgeeTrK6qg624ZiHWvFpYF3f4WuBpwaccwewA2Bqauq4JFCS1DPoC6Uurj7tcE1JkpZIkn+Y5OUz28Am4EvAbmBrK7YVuKNt7wauSs+FwLMzwzolSRrGnjxJkpbO2cAf956MwMnAH1bVnya5H7gtyTXAE8AVrfweeo9POEDvEQpXL/0lS5LGjUneInFugaSVxoekL7+qegz4iQHxvwEuHhAv4NoluDRJUoeY5EmSJEmaWF38EtQ5eZIkSZLUISZ5kiRJktQhDteUJElLxjnrkrT47MmTJEmSpA6xJ0+SJEmSjjHOC7LYkydJkiRJHWJP3gI5t0DSuBv0/9g4fEspSZIGsydPkiRJkjrEnjxJkrQoHO0iScvDnjxJkiRJ6hB78iRJ0oLYYydJK4s9eZIkSZLUISZ5kiRJktQhJnmSJEmS1CHOyTtBzjeQNEmG/Z/n8/MkSZNuHNpIe/IkSZIkqUNM8iRJkiSpQ0zyJEmSJKlDTPIkSZIkqUNM8iRJkiSpQ1xdU5IknRBXmpak4VbSqpsmeQPYiEnSYCupAZMkSYM5XFOSJEmSOsQkT5IkSZI6xCRPkiRJkjpkxczJS7IZ+CBwEvDfq+qGZb4kSZKW3XK0j85Nl6TxtiKSvCQnAR8GfgaYBu5Psruqvry8VyZJOhEuyLI4bB8lafwNaiMXu31cEUkecD5woKoeA0hyK7AFWNRGzG8qJUkr3LK0j5Kk8bZS5uStAZ7s+zzdYpIkTTLbR0nSnK2UnrwMiNVxhZJtwLb28dtJvrqoV7X0zgD+erkvYoWxTo5mfRzN+jjaiquPvG8kp/nHIznLeLJ97Flxv9srgHVyNOvjeNbJ0VZUfYyofYQhbeRKSfKmgXV9n9cCTx1bqKp2ADuW6qKWWpJ9VTW13NexklgnR7M+jmZ9HM366CTbR/zdHsQ6OZr1cTzr5GiTVh8rZbjm/cDGJOckeQlwJbB7ma9JkqTlZvsoSZqzFdGTV1UvJHkLcCe9JaJ3VtUjy3xZkiQtK9tHSdJ8rIgkD6Cq9gB7lvs6lllnh9osgHVyNOvjaNbH0ayPDrJ9BPzdHsQ6OZr1cTzr5GgTVR+pOm7+tiRJkiRpTK2UOXmSJEmSpBEwyVtCSf5rkq8keTjJHyc5rW/fdUkOJPlqkkv64ptb7ECS7X3xc5Lcl+TRJB9vE/LHSpIrkjyS5P8lmTpm38TVx2yG3XcXJdmZ5FCSL/XFTk+yt/377k2yqsWT5EOtXh5Ocl7fMVtb+UeTbF2Oe1moJOuS3J1kf/tbeWuLT2R9qJtsG49n+zg3k9JG2j4ezTbyRVSVryV6AZuAk9v2+4D3te1zgb8ETgXOAf6K3gT7k9r2DwIvaWXObcfcBlzZtj8C/Mpy39886uPHgB8B7gGm+uITWR+z1NPQ++7iC/gp4DzgS32x/wJsb9vb+/52LgM+Te9ZYhcC97X46cBj7X1V21613Pc2j7pYDZzXtl8O/I/29zGR9eGrmy/bxoF1Yvt44nU1MW2k7eNx9WEbOcvLnrwlVFV/VlUvtI/30nveEcAW4Naqer6qvgYcAM5vrwNV9VhV/R1wK7AlSYDXAre343cBly/VfYxKVe2vqkEP7J3I+pjFwPte5mtaNFX1OeDIMeEt9P5d4eh/3y3AzdVzL3BaktXAJcDeqjpSVc8Ae4HNi3/1o1VVB6vqwbb9LWA/sIYJrQ91k23j8Wwf52Ri2kjbx6PZRs7OJG/5/BK9bxOg9wv5ZN++6RYbFn8F8I2+RnEm3hXWx9GG3fckObuqDkLvP3XgrBaf6+/K2EqyAXg1cB/Wh7rLtnF21snxJv3/N9sDbCMHWTGPUOiKJJ8B/tGAXb9RVXe0Mr8BvAB8bOawAeWLwUl4zVJ+xTmR+hh02IBYJ+pjnrp+fwsxrG46VWdJXgZ8AnhbVX2z9+X84KIDYp2rD40f28bj2T6OzCTc43xMTHtgGzmYSSULL5kAAAHiSURBVN6IVdVPz7a/TeZ8PXBxVc38Ak0D6/qKrQWeatuD4n9Nr4v55PbtXH/5FeXF6mOIztbHPM1WH5Pi6SSrq+pgG1pxqMWH1c00cNEx8XuW4DpHLskp9Bqvj1XVJ1t4YutD48m28Xi2jyMz6W3kRLcHtpHDOVxzCSXZDLwD+Nmqeq5v127gyiSnJjkH2Ah8Abgf2NhWxnoJcCWwuzWAdwM/347fCgz71m8cWR9HG3jfy3xNS203vX9XOPrfdzdwVVsx60Lg2TY0405gU5JVbVWtTS02Vtp8mpuA/VX1/r5dE1kf6ibbxjmxTo436W3kxLYHtpEvYq4rtfia/4veBOkngYfa6yN9+36D3upQXwUu7YtfRm+1oL+iN4RjJv6D9P5jPwD8EXDqct/fPOrj5+h9e/I88DRw5yTXx4vU1cD77uILuAU4CPzf9vtxDb15JXcBj7b301vZAB9u9fJFjl6F7pfa78MB4Orlvq951sW/oDdk5OG+/zcum9T68NXNl23jwDqxfZxbfU1EG2n7eFx92EbO8kq7MUmSJElSBzhcU5IkSZI6xCRPkiRJkjrEJE+SJEmSOsQkT5IkSZI6xCRPkiRJkjrEJE+SJEmSOsQkT5IkSZI6xCRPkiRJkjrk/wPrle+U22ZogAAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 1080x360 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"a=sampler.to_pandas() \n",
"plt.subplot(1,2,1)\n",
"plt.hist(l1[:,0],weights=weights_np,bins=60);\n",
"plt.subplot(1,2,2)\n",
"plt.hist(a['p1x'].values[:],bins=60);\n",
"fig=plt.gcf()\n",
"fig.set_size_inches(15,5)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"p1x=a['p1x'].values[:].reshape(-1,1)\n",
"p1y=a['p1y'].values[:].reshape(-1,1)\n",
"p1z=a['p1z'].values[:].reshape(-1,1)\n",
"p1E=a['p1E'].values[:].reshape(-1,1)\n",
"\n",
"p2x=a['p2x'].values[:].reshape(-1,1)\n",
"p2y=a['p2y'].values[:].reshape(-1,1)\n",
"p2z=a['p2z'].values[:].reshape(-1,1)\n",
"p2E=a['p2E'].values[:].reshape(-1,1)\n",
"\n",
"qx =a['qx'].values[:].reshape(-1,1)\n",
"qy =a['qy'].values[:].reshape(-1,1)\n",
"qz =a['qz'].values[:].reshape(-1,1)\n",
"qE =a['qE'].values[:].reshape(-1,1)\n",
"\n",
"pKstx =a['pKstx'].values[:].reshape(-1,1)\n",
"pKsty =a['pKsty'].values[:].reshape(-1,1)\n",
"pKstz =a['pKstz'].values[:].reshape(-1,1)\n",
"pKstE =a['pKstE'].values[:].reshape(-1,1)\n",
"\n",
"p1=np.concatenate([p1x,p1y,p1z,p1E], axis=1)\n",
"p2=np.concatenate([p2x,p2y,p2z,p2E], axis=1)\n",
"q = np.concatenate([qx,qy,qz,qE], axis=1)\n",
"pKst = np.concatenate([pKstx,pKsty,pKstz,pKstE], axis=1)\n",
"#q = p1 + p2"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"q2_from_sample=scalar_product_4_np(q,q)"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA3AAAAEvCAYAAAAErSPcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAatElEQVR4nO3df5Bd51kf8O+DjE0JJTGOoES2kacxnXEopaDIMJQ0A7VxCETQ2kQJUxzqjpqC2zIMBac/jDEw4/ArQxtPi4tdQgLYIW1azUQgAunQDpOkkk1IkI2JcEW8diZxsOvUhBCUPP3jXoeb5e7ukXe1d4/u5zOj0Tnnfc/d9/r43L1fve953+ruAAAAsPN91qIbAAAAwDACHAAAwEgIcAAAACMhwAEAAIyEAAcAADASAhwAAMBInLfoBqz23Oc+t/fu3bvoZgAAACzEvffe+5Hu3j2vbMcFuL179+b48eOLbgYAAMBCVNUfrVVmCCUAAMBICHAAAAAjIcABAACMhAAHAAAwEgIcAADASAhwAAAAIyHAAQAAjIQABwAAMBICHAAAwEgIcAAAACMhwAEAAIzEeYtuAMnem962bvmp2166TS0BAAB2Mj1wAAAAIyHAAQAAjIQABwAAMBKegdsinmMDAADONj1wAAAAIyHAAQAAjIQABwAAMBICHAAAwEiYxGQENpogZT0mTwEAgHOHHjgAAICREOAAAABGQoADAAAYCQEOAABgJAQ4AACAkTAL5TbZzEySAAAAiR44AACA0RgU4Krqmqp6sKpOVtVNc8pfVFX3VdXpqrp2TvnnV9UjVfX6rWg0AADAMtowwFXVriS3J3lJkiuSvKKqrlhV7QNJXpXkl9Z4mR9J8lvPvJkAAAAMeQZuf5KT3f1QklTV3UkOJLn/6QrdfWpa9qnVJ1fVVyX5oiS/lmTf5pvMmdjo2btTt710m1oCAABs1pAhlHuSPDyzvzI9tqGq+qwkP5XkX25Q71BVHa+q44899tiQlwYAAFg6QwJczTnWA1//u5Mc6e6H16vU3Xd0977u3rd79+6BLw0AALBchgyhXElyycz+xUkeHfj6X5Pk66rqu5N8XpLzq+qp7v5LE6EAAACwviEB7liSy6vqsiSPJDmY5JVDXry7v+Pp7ap6VZJ9wtvOst4zcp6PAwCAnWXDIZTdfTrJjUmOJnkgyZu7+0RV3VpVL0uSqnphVa0kuS7Jz1bVibPZaAAAgGU0pAcu3X0kyZFVx26e2T6WydDK9V7j55P8/Bm3EAAAgCQDF/IGAABg8QQ4AACAkRDgAAAARmLQM3Asp/VmqEzMUgkAANtNDxwAAMBICHAAAAAjIcABAACMhAAHAAAwEgIcAADASAhwAAAAIyHAAQAAjIQABwAAMBICHAAAwEict+gGMF57b3rbuuWnbnvpNrUEAACWgx44AACAkdADN9BGvU0AAABnmx44AACAkRDgAAAARkKAAwAAGAkBDgAAYCQEOAAAgJEwCyULYQ05AAA4cwIcZ42lFwAAYGsZQgkAADASAhwAAMBIDApwVXVNVT1YVSer6qY55S+qqvuq6nRVXTtz/Cuq6p1VdaKq3ltVL9/KxgMAACyTDZ+Bq6pdSW5PclWSlSTHqupwd98/U+0DSV6V5PtXnf6xJN/Z3e+vquclubeqjnb3/92S1rOUTIACAMCyGjKJyf4kJ7v7oSSpqruTHEjy6QDX3aemZZ+aPbG7/2Bm+9Gq+nCS3UkEONZlAhQAAPjLhgyh3JPk4Zn9lemxM1JV+5Ocn+QP55QdqqrjVXX8scceO9OXBgAAWApDAlzNOdZn8kOq6ouTvDHJd3X3p1aXd/cd3b2vu/ft3r37TF4aAABgaQwJcCtJLpnZvzjJo0N/QFV9fpK3Jfk33f2uM2seAAAATxsS4I4lubyqLquq85McTHJ4yItP6781yS90968882YCAACw4SQm3X26qm5McjTJriR3dfeJqro1yfHuPlxVL8wkqF2Y5Fuq6oe7+wVJvj3Ji5JcVFWvmr7kq7r7PWfjzUBilkoAAM5dQ2ahTHcfSXJk1bGbZ7aPZTK0cvV5b0rypk22EQAAgAxcyBsAAIDFE+AAAABGQoADAAAYCQEOAABgJAQ4AACAkRg0CyWcS9ZbZsASAwAA7GR64AAAAEZCgAMAABgJAQ4AAGAkBDgAAICRMIkJnAEToAAAsEh64AAAAEZCgAMAABgJAQ4AAGAkBDgAAICRMIkJbJH1JjhJTHICAMDmCXAwY6MQBgAAi2QIJQAAwEgIcAAAACMhwAEAAIyEAAcAADASJjGBbWKWSgAANkuAgx1ivYAn3AEAkBhCCQAAMBoCHAAAwEgMCnBVdU1VPVhVJ6vqpjnlL6qq+6rqdFVdu6rs+qp6//TP9VvVcAAAgGWzYYCrql1Jbk/ykiRXJHlFVV2xqtoHkrwqyS+tOvcLkvxQkiuT7E/yQ1V14eabDQAAsHyG9MDtT3Kyux/q7k8kuTvJgdkK3X2qu9+b5FOrzv3GJG/v7se7+4kkb09yzRa0GwAAYOkMmYVyT5KHZ/ZXMulRG2LeuXsGngtMWYIAAIBkWA9czTnWA19/0LlVdaiqjlfV8ccee2zgSwMAACyXIQFuJcklM/sXJ3l04OsPOre77+jufd29b/fu3QNfGgAAYLkMCXDHklxeVZdV1flJDiY5PPD1jya5uqounE5ecvX0GAAAAGdowwDX3aeT3JhJ8HogyZu7+0RV3VpVL0uSqnphVa0kuS7Jz1bViem5jyf5kUxC4LEkt06PAQAAcIaGTGKS7j6S5MiqYzfPbB/LZHjkvHPvSnLXJtoIAABABi7kDQAAwOIJcAAAACMxaAglsLNZJw4AYDnogQMAABgJAQ4AAGAkBDgAAICREOAAAABGQoADAAAYCQEOAABgJAQ4AACAkbAOHCw5a8gBAIyHHjgAAICREOAAAABGQoADAAAYCQEOAABgJExiAqxrvUlOTHACALC99MABAACMhB44WAIbLRUAAMA46IEDAAAYCQEOAABgJAQ4AACAkRDgAAAARkKAAwAAGAkBDgAAYCQEOAAAgJEYFOCq6pqqerCqTlbVTXPKL6iqe6bl766qvdPjn11Vb6iq91XVA1X1mq1tPgAAwPLYcCHvqtqV5PYkVyVZSXKsqg539/0z1W5I8kR3P7+qDiZ5bZKXJ7kuyQXd/Ter6nOT3F9Vv9zdp7b6jQDbb6MFwk/d9tJtagkAwHLYMMAl2Z/kZHc/lCRVdXeSA0lmA9yBJLdMt9+S5PVVVUk6ybOq6rwkfyXJJ5J8dGuaDux0Ah4AwNYaMoRyT5KHZ/ZXpsfm1unu00meTHJRJmHuT5J8MMkHkvxkdz++yTYDAAAspSE9cDXnWA+ssz/JJ5M8L8mFSf5XVf3G0715nz656lCSQ0ly6aWXDmgScC7YqIduPXrvAIBlNKQHbiXJJTP7Fyd5dK060+GSz07yeJJXJvm17v7z7v5wkt9Osm/1D+juO7p7X3fv271795m/CwAAgCUwJMAdS3J5VV1WVecnOZjk8Ko6h5NcP92+Nsk7urszGTb59TXxrCRfneT3t6bpAAAAy2XDADd9pu3GJEeTPJDkzd19oqpuraqXTavdmeSiqjqZ5PuSPL3UwO1JPi/J72USBP9zd793i98DAADAUhjyDFy6+0iSI6uO3Tyz/fFMlgxYfd5T844DAABw5gYt5A0AAMDiCXAAAAAjIcABAACMhAAHAAAwEgIcAADASAhwAAAAIyHAAQAAjMSgdeAAziV7b3rbuuWnbnvpNrUEAODMCHDAKK0XwjYbwM7mawMAbIYAB5xzNuphAwAYK8/AAQAAjIQABwAAMBICHAAAwEgIcAAAACMhwAEAAIyEWSgBtpAlCACAs0kPHAAAwEjogQM4A9aYAwAWSQ8cAADASAhwAAAAI2EIJcA22ezwS5OgAAB64AAAAEZCgAMAABgJAQ4AAGAkBDgAAICRGBTgquqaqnqwqk5W1U1zyi+oqnum5e+uqr0zZV9eVe+sqhNV9b6q+pytaz4AAMDy2DDAVdWuJLcneUmSK5K8oqquWFXthiRPdPfzk7wuyWun556X5E1JXt3dL0jy4iR/vmWtBwAAWCJDeuD2JznZ3Q919yeS3J3kwKo6B5K8Ybr9liTfUFWV5Ook7+3u302S7v7j7v7k1jQdAABguQwJcHuSPDyzvzI9NrdOd59O8mSSi5J8aZKuqqNVdV9V/cDmmwwAALCchizkXXOO9cA65yX5O0lemORjSX6zqu7t7t/8jJOrDiU5lCSXXnrpgCYBAAAsnyEBbiXJJTP7Fyd5dI06K9Pn3p6d5PHp8d/q7o8kSVUdSfKVST4jwHX3HUnuSJJ9+/atDocAbGDvTW9bt/zUbS/dppYAAGfTkAB3LMnlVXVZkkeSHEzyylV1Die5Psk7k1yb5B3d3VV1NMkPVNXnJvlEkr+bySQnAJyhjUIaAHDu2zDAdffpqroxydEku5Lc1d0nqurWJMe7+3CSO5O8sapOZtLzdnB67hNV9dOZhMBOcqS7fQMBAAB4Bob0wKW7jyQ5surYzTPbH09y3RrnvimTpQQAAADYhEELeQMAALB4AhwAAMBIDBpCCcC4rTcBihkqAWA89MABAACMhAAHAAAwEgIcAADASAhwAAAAI2ESEwDWtd4EKIlJUABgOwlwAEtuo4C2mfOFOwDYWoZQAgAAjIQABwAAMBICHAAAwEgIcAAAACMhwAEAAIyEWSgBOGssQQAAW0sPHAAAwEgIcAAAACNhCCUAC2OIJQCcGT1wAAAAIyHAAQAAjIQhlADsWBsNsVyP4ZcAnIv0wAEAAIyEAAcAADASAhwAAMBICHAAAAAjMSjAVdU1VfVgVZ2sqpvmlF9QVfdMy99dVXtXlV9aVU9V1fdvTbMBAACWz4YBrqp2Jbk9yUuSXJHkFVV1xapqNyR5orufn+R1SV67qvx1SX51880FAABYXkOWEdif5GR3P5QkVXV3kgNJ7p+pcyDJLdPttyR5fVVVd3dVfWuSh5L8yZa1GgA2sN4SBJYYAGCshgS4PUkentlfSXLlWnW6+3RVPZnkoqr60yQ/mOSqJIZPArAjbLS+nIAHwE415Bm4mnOsB9b54SSv6+6n1v0BVYeq6nhVHX/ssccGNAkAAGD5DOmBW0lyycz+xUkeXaPOSlWdl+TZSR7PpKfu2qr68STPSfKpqvp4d79+9uTuviPJHUmyb9++1eEQAACADAtwx5JcXlWXJXkkycEkr1xV53CS65O8M8m1Sd7R3Z3k656uUFW3JHlqdXgDAABgmA0D3PSZthuTHE2yK8ld3X2iqm5Ncry7Dye5M8kbq+pkJj1vB89mowFgkUyQAsCiDOmBS3cfSXJk1bGbZ7Y/nuS6DV7jlmfQPgAAAKYGBTgAWCYbzVIJAIsyZBZKAAAAdgA9cACwhawxB8DZpAcOAABgJPTAAcA2MoMlAJuhBw4AAGAkBDgAAICRMIQSAEbCBCkA6IEDAAAYCQEOAABgJAyhBIAdYqMhkgAgwAHAOcISBQDnPkMoAQAARkKAAwAAGAkBDgAAYCQEOAAAgJEQ4AAAAEbCLJQAwLo2Wt7ADJcA20eAA4AlIIQBnBsEOABgU4RDgO0jwAEAG4YwAHYGk5gAAACMhB44AOCs2kzvnuGXAJ9JDxwAAMBICHAAAAAjMSjAVdU1VfVgVZ2sqpvmlF9QVfdMy99dVXunx6+qqnur6n3Tv79+a5sPAACwPDZ8Bq6qdiW5PclVSVaSHKuqw919/0y1G5I80d3Pr6qDSV6b5OVJPpLkW7r70ar6siRHk+zZ6jcBAJyb1nt+zvNxwDIaMonJ/iQnu/uhJKmqu5McSDIb4A4kuWW6/ZYkr6+q6u7fmalzIsnnVNUF3f1nm245ALDUrD8HLKMhAW5Pkodn9leSXLlWne4+XVVPJrkokx64p/2DJL8jvAEA20HAA85FQwJczTnWZ1Knql6QybDKq+f+gKpDSQ4lyaWXXjqgSQAAAMtnyCQmK0kumdm/OMmja9WpqvOSPDvJ49P9i5O8Ncl3dvcfzvsB3X1Hd+/r7n27d+8+s3cAAACwJIb0wB1LcnlVXZbkkSQHk7xyVZ3DSa5P8s4k1yZ5R3d3VT0nyduSvKa7f3vrmg0AcPYYfgnsVBv2wHX36SQ3ZjKD5ANJ3tzdJ6rq1qp62bTanUkuqqqTSb4vydNLDdyY5PlJ/m1VvWf65wu3/F0AAAAsgSE9cOnuI0mOrDp288z2x5NcN+e8H03yo5tsIwAAABkY4AAAzjUbDZME2IkEOACAM+QZOWBRBDgAgJEQHAEBDgBgi60XtIQsYDMEOACAbaQXDdgMAQ4AYAcxuQqwHgEOAOAcoXcPzn0bLuQNAADAziDAAQAAjIQABwAAMBICHAAAwEiYxAQAYEksaoZLk6fA1tEDBwAAMBICHAAAwEgYQgkAwFl1NoduGp7JshHgAAAYrfXCoXDHucgQSgAAgJEQ4AAAAEbCEEoAAM5JGz17Z4glYyTAAQCwlAQ8xkiAAwCAMyT8sSgCHAAAzLGZ5Q/MjsnZIsABAMA20nvHZghwAACwgwh4rEeAAwCAEdnM0E7hb/wGBbiquibJzyTZleTnuvu2VeUXJPmFJF+V5I+TvLy7T03LXpPkhiSfTPLPu/volrUeAAAYTPgbvw0DXFXtSnJ7kquSrCQ5VlWHu/v+mWo3JHmiu59fVQeTvDbJy6vqiiQHk7wgyfOS/EZVfWl3f3Kr3wgAAHD2GNq5Mwzpgduf5GR3P5QkVXV3kgNJZgPcgSS3TLffkuT1VVXT43d3958l+T9VdXL6eu/cmuYDAAA7wWZ69zYiHP6FIQFuT5KHZ/ZXkly5Vp3uPl1VTya5aHr8XavO3fOMWwsAACydsxUOxxgMhwS4mnOsB9YZcm6q6lCSQ9Pdp6rqwQHt2k7PTfKRRTeCz+Ca7Dyuyc7jmuw8rsnO45rsPK7JznPOXpN67aJbsKYvWatgSIBbSXLJzP7FSR5do85KVZ2X5NlJHh94brr7jiR3DGjLQlTV8e7et+h28Bdck53HNdl5XJOdxzXZeVyTncc12Xlck53lswbUOZbk8qq6rKrOz2RSksOr6hxOcv10+9ok7+junh4/WFUXVNVlSS5P8r+3pukAAADLZcMeuOkzbTcmOZrJMgJ3dfeJqro1yfHuPpzkziRvnE5S8ngmIS/Tem/OZMKT00m+xwyUAAAAz8ygdeC6+0iSI6uO3Tyz/fEk161x7o8l+bFNtHEn2LHDO5eYa7LzuCY7j2uy87gmO49rsvO4JjuPa7KD1GSkIwAAADvdkGfgAAAA2AEEuKmquqaqHqyqk1V105zyC6rqnmn5u6tq7/a3cnlU1SVV9T+q6oGqOlFV/2JOnRdX1ZNV9Z7pn5vnvRZbq6pOVdX7pv/Nj88pr6r6d9N75b1V9ZWLaOcyqKq/MfP//3uq6qNV9b2r6rhPtkFV3VVVH66q35s59gVV9faqev/07wvXOPf6aZ33V9X18+pw5ta4Jj9RVb8//Wx6a1U9Z41z1/2c45lZ45rcUlWPzHxGfdMa5677PY1nZo1rcs/M9ThVVe9Z41z3yYIYQpmkqnYl+YMkV2Wy9MGxJK/o7vtn6nx3ki/v7ldX1cEk39bdL19Ig5dAVX1xki/u7vuq6q8muTfJt666Ji9O8v3d/c0LauZSqqpTSfZ199z1YKa/fP9Zkm9KcmWSn+nuK7evhctp+jn2SJIru/uPZo6/OO6Ts66qXpTkqSS/0N1fNj3240ke7+7bpl84L+zuH1x13hckOZ5kXybrpN6b5Ku6+4ltfQPnoDWuydWZzJR9umqy+tPqazKtdyrrfM7xzKxxTW5J8lR3/+Q65234PY1nZt41WVX+U0me7O5b55SdivtkIfTATexPcrK7H+ruTyS5O8mBVXUOJHnDdPstSb6hquYtVM4W6O4Pdvd90+3/l+SBJHsW2yoGOpDJL4Lu7nclec40kHN2fUOSP5wNb2yf7v6fmczCPGv298YbknzrnFO/Mcnbu/vxaWh7e5JrzlpDl8i8a9Ldv97dp6e778pkfVq2yRr3yRBDvqfxDKx3Tabfc789yS9va6PYkAA3sSfJwzP7K/nLYeHTdaYf/k8muWhbWrfkpsNV/3aSd88p/pqq+t2q+tWqesG2Nmx5dZJfr6p7q+rQnPIh9xNb72DW/iXrPlmML+ruDyaTf5RK8oVz6rhfFucfJfnVNco2+pxja904HdZ61xpDjd0ni/F1ST7U3e9fo9x9siAC3MS8nrTVY0uH1GGLVdXnJfkvSb63uz+6qvi+JF/S3X8ryb9P8t+2u31L6mu7+yuTvCTJ90yHX8xyr2yzqjo/ycuS/MqcYvfJzuZ+WYCq+teZrE/7i2tU2ehzjq3zH5L89SRfkeSDSX5qTh33yWK8Iuv3vrlPFkSAm1hJcsnM/sVJHl2rTlWdl+TZeWbDABioqj47k/D2i939X1eXd/dHu/up6faRJJ9dVc/d5mYune5+dPr3h5O8NZOhLbOG3E9srZckua+7P7S6wH2yUB96evjw9O8Pz6njftlm04livjnJd/QaEwEM+Jxji3T3h7r7k939qST/KfP/W7tPttn0u+7fT3LPWnXcJ4sjwE0cS3J5VV02/Zfsg0kOr6pzOMnTs4Ndm8lD0P715yyZjru+M8kD3f3Ta9T5a08/h1hV+zP5//mPt6+Vy6eqnjWdVCZV9awkVyf5vVXVDif5zpr46kwefv7gNjd12az5r6Tuk4Wa/b1xfZL/PqfO0SRXV9WF06FjV0+PcRZU1TVJfjDJy7r7Y2vUGfI5xxZZ9Yz0t2X+f+sh39PYWn8vye9398q8QvfJYp236AbsBNPZqG7M5JfmriR3dfeJqro1yfHuPpxJmHhjVZ3MpOft4OJavBS+Nsk/TPK+melr/1WSS5Oku/9jJkH6n1bV6SR/muSgUH3WfVGSt07zwHlJfqm7f62qXp18+rocyWQGypNJPpbkuxbU1qVQVZ+bycxs/2Tm2Oz1cJ9sg6r65SQvTvLcqlpJ8kNJbkvy5qq6IckHklw3rbsvyau7+x939+NV9SOZfEFNklu72+iOLbDGNXlNkguSvH36Ofau6ezSz0vyc939TVnjc24Bb+Gcs8Y1eXFVfUUmQyJPZfpZNntN1vqetoC3cM6Zd026+87Mea7afbJzWEYAAABgJAyhBAAAGAkBDgAAYCQEOAAAgJEQ4AAAAEZCgAMAABgJAQ4AAGAkBDgAAICREOAAAABG4v8DIBVxR6i8BYYAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 1080x360 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.hist(q2_from_sample/1e6,bins=100, density=True);\n",
"fig=plt.gcf()\n",
"fig.set_size_inches(15,5);"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA3MAAAEvCAYAAADvmpjfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAWxklEQVR4nO3df4ylV3kf8O8TbyAN+WGDDTJru+s221YUtYGuDCpSiuKW+EfapRKO7FZhcSxtK5mENJXKkv7hKGmqpU2hREqRHOzGVMTGIkS2aivg8kOof9hhbRBgLOKt2dqLXbzI4ITSBJk8/eO+C8N6ZnZn7uydOTOfj7Sae8977r1n1mffma+f8563ujsAAACM5Qc2ewAAAACsnTAHAAAwIGEOAABgQMIcAADAgIQ5AACAAQlzAAAAA9q12QNYzfnnn9979uzZ7GEAAABsigcffPBr3X3Bcse2dJjbs2dPjhw5stnDAAAA2BRV9b9XOmaZJQAAwICEOQAAgAEJcwAAAAMS5gAAAAYkzAEAAAxImAMAABiQMAcAADAgYQ4AAGBAwhwAAMCAhDkAAIABCXMAAAAD2rXZAwAAgJP2HLpnxWPHDl+9wJHA1ifMAQCcRcIJozJ3tz5hDgAATkOwYSsS5gAAWBPBZuvw32JnswEKAADAgFTmAAAYwmpVqEQlip1HmAMA2CTzhBPBBhDmAIAhuDZoDEImiX+viyLMAcA25BcpWJvThdARCdbbnzAHALADCfwbZ54g6O+aeQhzAACwSVTPmIcwBwBsGL+YbqztuPQP2DjCHADANiQIwvYnzAEAw1MRhMWa538W+Pe6cYQ5AADYgVRvx/cDmz0AAAAA1k5lDgA2ia3htwfVjTH478R2JMwBADueYM1WJYSyGmEOAOYgBKyNvy+AjXPaMFdVtyb52SRPd/crp7b/mOQfJ/l2kv+V5Pru/sZ07B1JbkjynSS/1N0fmdqvSPKeJOckeV93H974bwcAYGPtxMrIqN/zqOOG9TqTDVB+L8kVp7Tdl+SV3f13kvxJknckSVW9Ism1Sf729Jr/UlXnVNU5SX4nyZVJXpHkuqkvAAAA63Daylx3f6qq9pzS9tElT+9P8qbp8f4kd3T3XyT5clUdTXLZdOxodz+WJFV1x9T3i3ONHgBYKJUP5mUOwcbZiGvmfiHJB6fHuzMLdycdn9qS5IlT2l+zAZ8NAFvWqL+0bta4R/37Atgsc4W5qvq3SZ5L8oGTTct06yy/nLNXeM+DSQ4mySWXXDLP8AAAkgiKwPa07jBXVQcy2xjl8u4+GcyOJ7l4SbeLkjw5PV6p/ft0981Jbk6Sffv2LRv4AGCp0/2ibpdEALajdYW5aWfKtyf5B939rSWH7k7y+1X1riQvT7I3yR9nVrHbW1WXJvlKZpuk/LN5Bg4AZ8p2+ABsR2dya4Lbk7w+yflVdTzJTZntXvnCJPdVVZLc393/srsfrqo7M9vY5LkkN3b3d6b3eWuSj2R2a4Jbu/vhs/D9AAAA7Ahnspvldcs037JK/99M8pvLtN+b5N41jQ4AWNY814C5fgwYmdUW33Mm95kDAABgi9mIWxMAwLZlm34AtiqVOQAAgAGpzAGwEG4fAAAbS2UOAABgQMIcAADAgCyzBGBHs9EIAKNSmQMAABiQyhzAgNwwFQAQ5gA4Y3akBICtwzJLAACAAanMATAEG5UAwPdTmQMAABiQMAcAADAgyywBAIAtw7L6MyfMAewwbmsAANuDZZYAAAADEuYAAAAGJMwBAAAMyDVzAFuQi78BgNMR5gDYEgRYAFgbYQ6A7xKoAGAcrpkDAAAYkMocwFmyE+/nprIHAIujMgcAADAgYQ4AAGBAwhwAAMCAhDkAAIABCXMAAAADEuYAAAAG5NYEwI52uq30t+stBACA8Z22MldVt1bV01X1hSVtL66q+6rq0enreVN7VdVvV9XRqvpcVb16yWsOTP0fraoDZ+fbAQAA2BnOZJnl7yW54pS2Q0k+1t17k3xsep4kVybZO/05mOS9ySz8JbkpyWuSXJbkppMBEAAAgLU77TLL7v5UVe05pXl/ktdPj29L8skkb5/a39/dneT+qjq3qi6c+t7X3c8kSVXdl1lAvH3u7wBgQKdb3gkAcDrrvWbuZd39VJJ091NV9dKpfXeSJ5b0Oz61rdT+PFV1MLOqXi655JJ1Dg/g7BPIAIDNtNG7WdYybb1K+/Mbu2/u7n3dve+CCy7Y0MEBAABsF+utzH21qi6cqnIXJnl6aj+e5OIl/S5K8uTU/vpT2j+5zs8GAAB4np22S/V6w9zdSQ4kOTx9vWtJ+1ur6o7MNjt5dgp8H0ny75dsevKGJO9Y/7ABWInlnwCwM5w2zFXV7ZlV1c6vquOZ7Up5OMmdVXVDkseTXDN1vzfJVUmOJvlWkuuTpLufqarfSPLpqd+vn9wMBQAAgLU7k90sr1vh0OXL9O0kN67wPrcmuXVNowMAAGBZG70BCgAAAAsgzAEAAAxImAMAABiQMAcAADCg9d6aAGBHsM0/ALBVqcwBAAAMSGUOGN5q1bNjh69e4EgAABZHZQ4AAGBAKnPAlue6NQCA5xPmgC1BYAMAWBvLLAEAAAYkzAEAAAxImAMAABiQMAcAADAgYQ4AAGBAdrMEtjW7ZAIA25UwByyEUAUAsLEsswQAABiQMAcAADAgYQ4AAGBAwhwAAMCAhDkAAIABCXMAAAADEuYAAAAG5D5zwBlb7V5xxw5fvcCRAACgMgcAADAglTngu1arvAEAsLWozAEAAAxIZQ7YEKp6AACLpTIHAAAwoLnCXFX9q6p6uKq+UFW3V9UPVdWlVfVAVT1aVR+sqhdMfV84PT86Hd+zEd8AAADATrTuMFdVu5P8UpJ93f3KJOckuTbJO5O8u7v3Jvl6khuml9yQ5Ovd/RNJ3j31AwAAYB3mXWa5K8lfqapdSX44yVNJfjrJh6bjtyV54/R4//Q80/HLq6rm/HwAAIAdad0boHT3V6rqt5I8nuT/JflokgeTfKO7n5u6HU+ye3q8O8kT02ufq6pnk7wkydfWOwbg+U63EYmbewMAbA/zLLM8L7Nq26VJXp7kRUmuXKZrn3zJKseWvu/BqjpSVUdOnDix3uEBAABsa/PcmuAfJvlyd59Ikqr6cJK/n+Tcqto1VecuSvLk1P94kouTHJ+WZf54kmdOfdPuvjnJzUmyb9++54U92OncAgAAgGS+a+YeT/Laqvrh6dq3y5N8Mcknkrxp6nMgyV3T47un55mOf7y7hTUAAIB1mOeauQeq6kNJHkryXJLPZFZRuyfJHVX176a2W6aX3JLkv1XV0cwqctfOM3BgfVT2AAC2h3mWWaa7b0py0ynNjyW5bJm+f57kmnk+DwAAgJl5b00AAADAJhDmAAAABiTMAQAADEiYAwAAGJAwBwAAMCBhDgAAYEDCHAAAwICEOQAAgAEJcwAAAAMS5gAAAAYkzAEAAAxImAMAABjQrs0eAGxHew7ds+rxY4evnuv1AACgMgcAADAglTnYBCpvAADMS2UOAABgQMIcAADAgIQ5AACAAQlzAAAAA7IBCqyTTUwAANhMKnMAAAADEuYAAAAGJMwBAAAMSJgDAAAYkDAHAAAwIGEOAABgQG5NACtw6wEAgO1ltd/vjh2+eoEj2RgqcwAAAAMS5gAAAAYkzAEAAAxormvmqurcJO9L8sokneQXknwpyQeT7ElyLMnPdffXq6qSvCfJVUm+leQt3f3QPJ8P83JdHAAAo5q3MveeJH/U3X8ryd9N8kiSQ0k+1t17k3xsep4kVybZO/05mOS9c342AADAjrXuMFdVP5bkp5LckiTd/e3u/kaS/Ulum7rdluSN0+P9Sd7fM/cnObeqLlz3yAEAAHaweSpzfy3JiST/tao+U1Xvq6oXJXlZdz+VJNPXl079dyd5Ysnrj09tAAAArNE8YW5XklcneW93vyrJ/833llQup5Zp6+d1qjpYVUeq6siJEyfmGB4AAMD2NU+YO57keHc/MD3/UGbh7qsnl09OX59e0v/iJa+/KMmTp75pd9/c3fu6e98FF1wwx/AAAAC2r3WHue7+P0meqKq/OTVdnuSLSe5OcmBqO5Dkrunx3UneXDOvTfLsyeWYAAAArM1ctyZI8otJPlBVL0jyWJLrMwuId1bVDUkeT3LN1PfezG5LcDSzWxNcP+dnAwAA7Fhzhbnu/mySfcscunyZvp3kxnk+DwAAgJl5K3Nw1q12Y+9jh69e4EgAAGDrEOYY2mpBDwAAtrN5drMEAABgkwhzAAAAAxLmAAAABiTMAQAADEiYAwAAGJAwBwAAMCC3JmDTub0AAACsncocAADAgIQ5AACAAQlzAAAAAxLmAAAABiTMAQAADEiYAwAAGJAwBwAAMCBhDgAAYEBuGs6GON2Nv48dvnpBIwEAgJ1BZQ4AAGBAKnMsxOkqdwAAwNqozAEAAAxImAMAABiQMAcAADAgYQ4AAGBAwhwAAMCAhDkAAIABCXMAAAADEuYAAAAG5KbhfNdqN/Y+dvjqBY4EAAA4HZU5AACAAc0d5qrqnKr6TFX99+n5pVX1QFU9WlUfrKoXTO0vnJ4fnY7vmfezAQAAdqqNqMy9LckjS56/M8m7u3tvkq8nuWFqvyHJ17v7J5K8e+oHAADAOswV5qrqoiRXJ3nf9LyS/HSSD01dbkvyxunx/ul5puOXT/0BAABYo3krc/85yb9J8pfT85ck+UZ3Pzc9P55k9/R4d5InkmQ6/uzUHwAAgDVad5irqp9N8nR3P7i0eZmufQbHlr7vwao6UlVHTpw4sd7hAQAAbGvz3JrgdUn+SVVdleSHkvxYZpW6c6tq11R9uyjJk1P/40kuTnK8qnYl+fEkz5z6pt19c5Kbk2Tfvn3PC3tsjtVuWwAAACzeuitz3f2O7r6ou/ckuTbJx7v7nyf5RJI3Td0OJLlrenz39DzT8Y93t7AGAACwDmfjPnNvT/IrVXU0s2vibpnab0nykqn9V5IcOgufDQAAsCPMs8zyu7r7k0k+OT1+LMlly/T58yTXbMTnAQAA7HRnozIHAADAWSbMAQAADEiYAwAAGJAwBwAAMCBhDgAAYEDCHAAAwICEOQAAgAEJcwAAAAMS5gAAAAa0a7MHwOLsOXTPZg8BAADYICpzAAAAAxLmAAAABiTMAQAADEiYAwAAGJAwBwAAMCBhDgAAYEBuTbDNuP0AAADsDCpzAAAAA1KZG4zKGwAAkKjMAQAADEmYAwAAGJAwBwAAMCBhDgAAYEDCHAAAwICEOQAAgAEJcwAAAAMS5gAAAAYkzAEAAAxImAMAABjQrs0eAN9vz6F7NnsIAADAAFTmAAAABrTuMFdVF1fVJ6rqkap6uKreNrW/uKruq6pHp6/nTe1VVb9dVUer6nNV9eqN+iYAAAB2mnmWWT6X5F9390NV9aNJHqyq+5K8JcnHuvtwVR1KcijJ25NcmWTv9Oc1Sd47fd1xLKUEAADmte7KXHc/1d0PTY//LMkjSXYn2Z/ktqnbbUneOD3en+T9PXN/knOr6sJ1jxwAAGAH25Br5qpqT5JXJXkgycu6+6lkFviSvHTqtjvJE0tednxqAwAAYI3mDnNV9SNJ/iDJL3f3n67WdZm2Xub9DlbVkao6cuLEiXmHBwAAsC3NFeaq6gczC3If6O4PT81fPbl8cvr69NR+PMnFS15+UZInT33P7r65u/d1974LLrhgnuEBAABsW/PsZllJbknySHe/a8mhu5McmB4fSHLXkvY3T7tavjbJsyeXYwIAALA28+xm+bokP5/k81X12antV5McTnJnVd2Q5PEk10zH7k1yVZKjSb6V5Po5PhsAAGBHW3eY6+7/meWvg0uSy5fp30luXO/nAQAA8D0bspslAAAAiyXMAQAADEiYAwAAGJAwBwAAMCBhDgAAYEDCHAAAwICEOQAAgAHNc9NwVrDn0D2bPQQAAGCbU5kDAAAYkDAHAAAwIGEOAABgQMIcAADAgIQ5AACAAQlzAAAAAxLmAAAABiTMAQAADEiYAwAAGJAwBwAAMCBhDgAAYEC7NnsAI9pz6J7NHgIAALDDqcwBAAAMSJgDAAAYkDAHAAAwIGEOAABgQMIcAADAgIQ5AACAAQlzAAAAAxLmAAAABiTMAQAADEiYAwAAGNDCw1xVXVFVX6qqo1V1aNGfDwAAsB0sNMxV1TlJfifJlUlekeS6qnrFIscAAACwHSy6MndZkqPd/Vh3fzvJHUn2L3gMAAAAw1t0mNud5Iklz49PbQAAAKzBrgV/Xi3T1t/XoepgkoPT029W1ZdO6X9+kq+dhbHBcsw3Fsl8Y5HMNxbJfGOR1jXf6p1nYSQb46+udGDRYe54kouXPL8oyZNLO3T3zUluXukNqupId+87O8OD72e+sUjmG4tkvrFI5huLtJPm26KXWX46yd6qurSqXpDk2iR3L3gMAAAAw1toZa67n6uqtyb5SJJzktza3Q8vcgwAAADbwaKXWaa7701y7xxvseISTDgLzDcWyXxjkcw3Fsl8Y5F2zHyr7j59LwAAALaURV8zBwAAwAYYKsxV1RVV9aWqOlpVhzZ7PGwPVXWsqj5fVZ+tqiNT24ur6r6qenT6et7UXlX129Mc/FxVvXpzR89WV1W3VtXTVfWFJW1rnl9VdWDq/2hVHdiM74Wtb4X59mtV9ZXpHPfZqrpqybF3TPPtS1X1M0va/bxlVVV1cVV9oqoeqaqHq+ptU7vzGxtulfm2489vwyyzrKpzkvxJkn+U2S0OPp3kuu7+4qYOjOFV1bEk+7r7a0va/kOSZ7r78PQP/bzufvt0kvjFJFcleU2S93T3azZj3Iyhqn4qyTeTvL+7Xzm1rWl+VdWLkxxJsi+ze3M+mOTvdffXN+FbYgtbYb79WpJvdvdvndL3FUluT3JZkpcn+R9J/sZ02M9bVlVVFya5sLsfqqofzey89MYkb4nzGxtslfn2c9nh57eRKnOXJTna3Y9197eT3JFk/yaPie1rf5Lbpse3ZXbCONn+/p65P8m50wkGltXdn0ryzCnNa51fP5Pkvu5+ZvoF574kV5z90TOaFebbSvYnuaO7/6K7v5zkaGY/a/285bS6+6nufmh6/GdJHkmyO85vnAWrzLeV7Jjz20hhbneSJ5Y8P57V/yPCmeokH62qB6vq4NT2su5+KpmdQJK8dGo3D9kIa51f5h3zeuu0tO3Wk8veYr6xQapqT5JXJXkgzm+cZafMt2SHn99GCnO1TNsYa0TZ6l7X3a9OcmWSG6dlSisxDzmbVppf5h3zeG+Sv57kJ5M8leQ/Te3mG3Orqh9J8gdJfrm7/3S1rsu0mW+syTLzbcef30YKc8eTXLzk+UVJntyksbCNdPeT09enk/xhZiX4r55cPjl9fXrqbh6yEdY6v8w71q27v9rd3+nuv0zyu5md4xLzjTlV1Q9m9ov1B7r7w1Oz8xtnxXLzzfltrDD36SR7q+rSqnpBkmuT3L3JY2JwVfWi6ULaVNWLkrwhyRcym1snd9Q6kOSu6fHdSd487cr12iTPnlxOAmuw1vn1kSRvqKrzpiUkb5ja4LROua73n2Z2jktm8+3aqnphVV2aZG+SP46ft5yBqqoktyR5pLvfteSQ8xsbbqX55vyW7NrsAZyp7n6uqt6a2T/wc5Lc2t0Pb/KwGN/Lkvzh7ByRXUl+v7v/qKo+neTOqrohyeNJrpn635vZTlxHk3wryfWLHzIjqarbk7w+yflVdTzJTUkOZw3zq7ufqarfyOyHUJL8enef6SYX7CArzLfXV9VPZraU6FiSf5Ek3f1wVd2Z5ItJnktyY3d/Z3ofP285ndcl+fkkn6+qz05tvxrnN86OlebbdTv9/DbMrQkAAAD4npGWWQIAADAR5gAAAAYkzAEAAAxImAMAABiQMAcAADAgYQ4AAGBAwhwAAMCAhDkAAIAB/X+qOUodAFnpegAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 1080x360 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.hist(p1[:,3],bins=100);\n",
"fig = plt.gcf()\n",
"fig.set_size_inches(15,5)"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA3MAAAEvCAYAAADvmpjfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAUaklEQVR4nO3dXail13kf8P9TKfZFE5BUyYrQR0ctc1EntE4YZEN64datPkPHKVWwofXYMagFGRrIRcZ2QCWJy4TSpE7riKp4aglcyyKJ66FWqkxFgtsLORq7xh9RVA/O1JpqkORIcVIELqqfXux33KOZ8/2xz15n/35wOHs/77v3XhvW7DP/vda7VnV3AAAAGMtf2O8GAAAAsHXCHAAAwICEOQAAgAEJcwAAAAMS5gAAAAYkzAEAAAzoyv1uwHquvfbaPnTo0H43AwAAYF988Ytf/HZ3X7fasYUOc4cOHcqZM2f2uxkAAAD7oqr+51rHTLMEAAAYkDAHAAAwIGEOAABgQMIcAADAgIQ5AACAAQlzAAAAAxLmAAAABiTMAQAADEiYAwAAGJAwBwAAMCBhDgAAYEBX7ncDAGCRHTr+uXWPnztxz5xaAgCvJ8wBcOCtF8iEMQBGZZolAADAgIzMAbDUNppGCQCLysgcAADAgIzMATA8o2sALCMjcwAAAAMS5gAAAAYkzAEAAAxImAMAABiQBVAAYI9stDCLDcsB2AkjcwAAAAMS5gAAAAYkzAEAAAxImAMAABiQMAcAADAgq1kCwA5stGIlAOwVI3MAAAADEuYAAAAGJMwBAAAMSJgDAAAYkDAHAAAwIGEOAABgQMIcAADAgOwzB8AQ7OcGAK+3YZirqpuTPJLkh5N8L8lD3f3RqromyaeTHEpyLslPd/crVVVJPprk7iSvJnlvd39peq5jSX5heupf7u6Hd/ftADAqYQ0AtmYzI3OvJfm57v5SVf1Qki9W1ekk703yZHefqKrjSY4n+fkkdyU5PP28NcmDSd46hb8HkhxJ0tPznOruV3b7TQHA6DYKt+dO3DOnlgCwqDa8Zq67L1wcWevuP0/yTJIbkxxNcnFk7eEk75xuH03ySM88leSqqrohyR1JTnf3y1OAO53kzl19NwAAAEtiSwugVNWhJD+W5AtJru/uC8ks8CV503TajUmeW/Gw81NtrToAAABbtOkwV1U/mOS3kvxsd//ZeqeuUut16pe+zn1Vdaaqzrz00kubbR4AAMBS2VSYq6ofyCzIfbK7f3sqvzBNn8z0+8Wpfj7JzSseflOS59epv053P9TdR7r7yHXXXbeV9wIAALA0NrOaZSX5eJJnuvtXVxw6leRYkhPT78+uqH+gqh7NbAGU73T3hap6Isk/r6qrp/NuT/LB3XkbADAeK3gCsBObWc3yJ5L8oyRfraovT7UPZRbiHquq9yf5VpJ7p2OPZ7YtwdnMtiZ4X5J098tV9UtJnp7O+8XufnlX3gUAC09wAYDdtWGY6+7/ltWvd0uSd6xyfie5f43nOpnk5FYaCAAAwOW2tJolAAAAi0GYAwAAGJAwBwAAMCBhDgAAYECbWc0SAFgwG60Oeu7EPXNqCQD7xcgcAADAgIQ5AACAAQlzAAAAAxLmAAAABiTMAQAADMhqlgDsmo1WWAQAdo8wBwBLxrYGAAeDMAfAphl5A4DFIcwBwAEkeAMcfBZAAQAAGJAwBwAAMCBhDgAAYEDCHAAAwICEOQAAgAEJcwAAAAMS5gAAAAYkzAEAAAxImAMAABiQMAcAADAgYQ4AAGBAV+53AwCYr0PHP7fmsXMn7pljS1hU6/WRjehDAPNjZA4AAGBAwhwAAMCATLMEOGB2MkVuJ48FAObLyBwAAMCAhDkAAIABmWYJAOyajabqWu0SYPcYmQMAABiQMAcAADAgYQ4AAGBAwhwAAMCAhDkAAIABCXMAAAADEuYAAAAGJMwBAAAMSJgDAAAYkDAHAAAwoCv3uwEAbN2h45/b7ybAtqzXd8+duGeOLQEYn5E5AACAARmZA1hARt7gckb1AF7PyBwAAMCAhDkAAIABCXMAAAADcs0cALAQXCsKsDUbjsxV1cmqerGqvrai9s+q6n9V1Zenn7tXHPtgVZ2tqmer6o4V9Tun2tmqOr77bwUAAGB5bGaa5SeS3LlK/de6+y3Tz+NJUlVvTvKuJD8yPeY3quqKqroiyceS3JXkzUnePZ0LAADANmw4zbK7P19Vhzb5fEeTPNrd303yx1V1Nslt07Gz3f3NJKmqR6dz/3DLLQYAAGBHC6B8oKq+Mk3DvHqq3ZjkuRXnnJ9qa9UvU1X3VdWZqjrz0ksv7aB5AAAAB9d2w9yDSf5qkrckuZDkX071WuXcXqd+ebH7oe4+0t1Hrrvuum02DwAA4GDb1mqW3f3CxdtV9e+S/Kfp7vkkN6849aYkz0+316oDLB2r9gEAO7WtkbmqumHF3Z9KcnGly1NJ3lVVb6yqW5McTvIHSZ5Ocriqbq2qN2S2SMqp7TcbAABguW04MldVn0ry9iTXVtX5JA8keXtVvSWzqZLnkvzjJOnur1fVY5ktbPJakvu7+/9Oz/OBJE8kuSLJye7++q6/GwAAgCWxmdUs371K+ePrnP+RJB9Zpf54kse31DoAAABWtZPVLAEAANgnwhwAAMCAhDkAAIABCXMAAAAD2tY+cwBszF5yMD8b/Xs7d+KeObUEYH6MzAEAAAxImAMAABiQMAcAADAgYQ4AAGBAFkAB2CYLnAAA+8nIHAAAwICMzAEAB56tC4CDyMgcAADAgIQ5AACAAZlmCQAsvfWmYZqCCSwqYQ5Yaq6jAQBGZZolAADAgIQ5AACAAZlmCRxoNvYGAA4qI3MAAAADEuYAAAAGJMwBAAAMSJgDAAAYkDAHAAAwIKtZAqzDapgAwKIS5gAA1rHRlzrnTtwzp5YAvJ4wBwCwA8IesF9cMwcAADAgI3MAAHtovZE7o3bAThiZAwAAGJAwBwAAMCBhDgAAYEDCHAAAwIAsgAIsPBt3AwBczsgcAADAgIQ5AACAAZlmCSwEUykBALZGmAMA2Cc7+SLLhuOAaZYAAAADEuYAAAAGJMwBAAAMSJgDAAAYkDAHAAAwIGEOAABgQMIcAADAgIQ5AACAAQlzAAAAAxLmAAAABnTlRidU1ckkP5nkxe7+0al2TZJPJzmU5FySn+7uV6qqknw0yd1JXk3y3u7+0vSYY0l+YXraX+7uh3f3rQAALI9Dxz+37vFzJ+6ZU0uA/bKZkblPJLnzktrxJE929+EkT073k+SuJIenn/uSPJh8P/w9kOStSW5L8kBVXb3TxgMAACyrDcNcd38+ycuXlI8muTiy9nCSd66oP9IzTyW5qqpuSHJHktPd/XJ3v5LkdC4PiAAAAGzSdq+Zu767LyTJ9PtNU/3GJM+tOO/8VFurDgAAwDZseM3cFtUqtV6nfvkTVN2X2RTN3HLLLbvXMmBPuXYDAGC+thvmXqiqG7r7wjSN8sWpfj7JzSvOuynJ81P97ZfUf3+1J+7uh5I8lCRHjhxZNfAB+2OjwLZXjwUA4HLbnWZ5Ksmx6faxJJ9dUX9PzbwtyXemaZhPJLm9qq6eFj65faoBAACwDZvZmuBTmY2qXVtV5zNblfJEkseq6v1JvpXk3un0xzPbluBsZlsTvC9JuvvlqvqlJE9P5/1id1+6qAoAAACbtGGY6+53r3HoHauc20nuX+N5TiY5uaXWAQAAsKrdXgAFAIABrHcts0WrYAzbvWYOAACAfWRkDgDgALKKMBx8whwAAK9j71AYg2mWAAAAAxLmAAAABiTMAQAADEiYAwAAGJAwBwAAMCCrWcKSsUksAMDBIMwB32dPIgCAcQhzAABsiVkesBhcMwcAADAgYQ4AAGBAwhwAAMCAhDkAAIABCXMAAAADEuYAAAAGZGsCOGDsFQfAftro75CtC2D3GJkDAAAYkDAHAAAwINMsAQCYG9MwYfcYmQMAABiQMAcAADAg0yxhAa03BcX0EwAAEmEOhmPrAQAAEmEOAIAFYnYKbJ5r5gAAAAYkzAEAAAxImAMAABiQMAcAADAgYQ4AAGBAwhwAAMCAbE0AAMCBYFsDlo0wBwDAENYLa7CMhDnYB/4YAQCwU66ZAwAAGJAwBwAAMCDTLGEPmEYJAItlo7/NFkhhREbmAAAABiTMAQAADEiYAwAAGJAwBwAAMCALoMA2WeQEAID9JMwBALD0rHbJiEyzBAAAGJCROZbaet/C+QYOAIBFZmQOAABgQEbmYA0WOAEAYJHtaGSuqs5V1Ver6stVdWaqXVNVp6vqG9Pvq6d6VdWvV9XZqvpKVf34brwBAACAZbQb0yz/Vne/pbuPTPePJ3myuw8neXK6nyR3JTk8/dyX5MFdeG0AAICltBfXzB1N8vB0++Ek71xRf6RnnkpyVVXdsAevDwAAcODt9Jq5TvK7VdVJ/m13P5Tk+u6+kCTdfaGq3jSde2OS51Y89vxUu7DyCavqvsxG7nLLLbfssHksO9e9AQC7wQrYLKKdhrmf6O7np8B2uqr+aJ1za5VaX1aYBcKHkuTIkSOXHQcAAGCHYa67n59+v1hVn0lyW5IXquqGaVTuhiQvTqefT3LziofflOT5nbw+AADst41mAhm5Y69s+5q5qvqLVfVDF28nuT3J15KcSnJsOu1Yks9Ot08lec+0quXbknzn4nRMAAAAtmYnI3PXJ/lMVV18nv/Q3f+5qp5O8lhVvT/Jt5LcO53/eJK7k5xN8mqS9+3gtQEAAJbatsNcd38zyd9Ypf4nSd6xSr2T3L/d1wMAgBFZPIW9shdbEwAAALDHhDkAAIAB7XRrAthX9pEDAGBZCXMAALBPbGvATghz7DsfYgAAsHXCHAvPVEoAALicMAcAAAvKDCbWYzVLAACAARmZAwCAQdmQfLkZmQMAABiQMAcAADAg0yyZCytSAgDA7jIyBwAAMCBhDgAAYEDCHAAAwICEOQAAgAFZAIVdYYETAACYL2EOAAAOoI2+bLep+PhMswQAABiQMAcAADAg0yzZFNfEAQDAYhHmAABgCbmmbnymWQIAAAzIyBwAAHCZ9UbujNotBmGO73NdHAAAm2GK5mIwzRIAAGBAwhwAAMCAhDkAAIABCXMAAAADsgDKErHACQAAHBzCHAAAsKtsazAfwtxgLAMLAAAkrpkDAAAYkpG5A8Z1cQAAsByMzAEAAAxImAMAABiQMAcAADAg18wtGNe8AQBwkFmdffcIcwAAwMIQ9jbPNEsAAIABCXMAAAADEuYAAAAGJMwBAAAMyAIoAADAMNZbIGXZFkcR5vaB7QcAAICdMs0SAABgQMIcAADAgEyz3AOmUQIAwPwt24bjcx+Zq6o7q+rZqjpbVcfn/foAAAAHwVxH5qrqiiQfS/J3k5xP8nRVneruP5xnO3bKyBsAALDf5j3N8rYkZ7v7m0lSVY8mOZpkqDAHAACM56BtazDvaZY3Jnluxf3zUw0AAIAtmPfIXK1S69edUHVfkvumu/+7qp7d81YxT9cm+fZ+N4IDSd9iL+hX7BV9i72ib21T/cp+t2BNf3mtA/MOc+eT3Lzi/k1Jnl95Qnc/lOSheTaK+amqM919ZL/bwcGjb7EX9Cv2ir7FXtG3lsu8p1k+neRwVd1aVW9I8q4kp+bcBgAAgOHNdWSuu1+rqg8keSLJFUlOdvfX59kGAACAg2Dum4Z39+NJHp/367IwTKFlr+hb7AX9ir2ib7FX9K0lUt298VkAAAAslHlfMwcAAMAuEObYNVX1L6rqj6rqK1X1maq6asWxD1bV2ap6tqruWFG/c6qdrarjK+q3VtUXquobVfXpacEcllRV3VtVX6+q71XVkUuO6VvsibX6EKylqk5W1YtV9bUVtWuq6vT0mXO6qq6e6lVVvz71r69U1Y+veMyx6fxvVNWx/XgvLI6qurmqfq+qnpn+Fv7Tqa5vIcyxq04n+dHu/utJ/keSDyZJVb05s5VLfyTJnUl+o6quqKorknwsyV1J3pzk3dO5SfIrSX6tuw8neSXJ++f6Tlg0X0vy95N8fmVR32KvbNCHYC2fyOyzaKXjSZ6cPnOenO4ns751ePq5L8mDyew/6EkeSPLWJLcleeDif9JZWq8l+bnu/mtJ3pbk/unzSN9CmGP3dPfvdvdr092nMttHMEmOJnm0u7/b3X+c5GxmHyK3JTnb3d/s7v+T5NEkR6uqkvztJL85Pf7hJO+c1/tg8XT3M9397CqH9C32yqp9aJ/bxILr7s8nefmS8tHMPmuS13/mHE3ySM88leSqqrohyR1JTnf3y939SmZflF4aEFki3X2hu7803f7zJM8kuTH6FhHm2Ds/k+R3pts3JnluxbHzU22t+l9K8qcrguHFOlxK32KvrNWHYKuu7+4Lyew/5UneNNW3+vkFqapDSX4syReib5F92JqAsVXVf0nyw6sc+nB3f3Y658OZTQn45MWHrXJ+Z/UvE3qd8znANtO3VnvYKjV9i92gr7DX1upj+h6rqqofTPJbSX62u/9sNtlk9VNXqelbB5Qwx5Z0999Z7/h0Me1PJnlH//99L84nuXnFaTcleX66vVr925lNCbhyGkFZeT4H1EZ9aw36Fntlvb4FW/FCVd3Q3RemqW4vTvW1+tj5JG+/pP77c2gnC6yqfiCzIPfJ7v7tqaxvYZolu6eq7kzy80n+Xne/uuLQqSTvqqo3VtWtmV2Q+wdJnk5yeFpd8A2ZLWRxagqBv5fkH0yPP5ZkrZEZlpu+xV5ZtQ/tc5sY06nMPmuS13/mnErynmnlwbcl+c40Ve6JJLdX1dXT4hS3TzWW1HS998eTPNPdv7rikL6FTcPZPVV1Nskbk/zJVHqqu//JdOzDmV1H91pm0wN+Z6rfneRfJbkiycnu/shU/yuZLThwTZL/nuQfdvd35/h2WCBV9VNJ/nWS65L8aZIvd/cd0zF9iz2xVh+CtVTVpzIb+bg2yQuZrRz4H5M8luSWJN9Kcm93vzz9B/3fZLYAxatJ3tfdZ6bn+ZkkH5qe9iPd/e/n+T5YLFX1N5P81yRfTfK9qfyhzK6b07eWnDAHAAAwINMsAQAABiTMAQAADEiYAwAAGJAwBwAAMCBhDgAAYEDCHAAAwICEOQAAgAEJcwAAAAP6f5pdPylm7h5fAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 1080x360 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.hist(p1[:,0],bins=100);\n",
"fig = plt.gcf()\n",
"fig.set_size_inches(15,5)"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA3MAAAEvCAYAAADvmpjfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAW8klEQVR4nO3db4ylV30f8O8v3kCL0sYGL5Su7a6jbKs4qBVoZdymahFOwcY061ZYMv3DCixZkZyGNJHqNbyw1BZpUSogNCmSFTuYiGIoofWqdkocA6J9YYc1IMA41Fvj2hu7eKn/JK0VqJNfX9xnw+x6dmZ37t2ZOTOfj7Sa5znn3HvP9R7fvd855zlPdXcAAAAYyw9tdAcAAAA4c8IcAADAgIQ5AACAAQlzAAAAAxLmAAAABiTMAQAADGjHRndgJeeff37v3r17o7sBAACwIR544IHvdvfO5eo2dZjbvXt3Dh8+vNHdAAAA2BBV9T9PVWeZJQAAwICEOQAAgAEJcwAAAAMS5gAAAAYkzAEAAAxImAMAABiQMAcAADAgYQ4AAGBAwhwAAMCAhDkAAIABCXMAAAAD2rHRHQAA2K52H7hrxfpHD161Tj0BRiTMAaxgpS9avmQBABvJMksAAIABCXMAAAADsswSYBNyHc2Z8d+LrcpSb2AlwhwAf0YoAoBxCHPA8FYLIBv1uoIPbA3+Xwc2K2EOYJvZqPC7mrP5hdlSNQC2ImEO4CwRIACAs0mYA9gAm3V2jBezxG5r8PcIbEXCHAAwhM062y0oAhtFmAMWwpcZAID1JcwBm54libB5+MXNmfMZBpwtwhwAQ/CFmJUYH8B2JMwBwBY0YrgZsc8AG0mYAzYFX+LOzGbdCIITzTOu/T0CsJpVw1xV3ZbkrUme6u7XTGW/nOTvJ/l+kv+R5J3d/exUd1OS65L8SZKf7+7PTuVXJPmVJOck+fXuPrj4twMA68svItgo816/6JdCML7TmZn7aJJfTfKxJWX3JLmpu1+oqvcnuSnJjVV1SZJrk/xkkr+c5Her6q9Oj/m1JH8vydEkX6qqQ939zcW8DQBYm1HD2Gbttw1SNo/NOkaAxfmh1Rp09xeTPH1S2e909wvT6X1JLpiO9yW5o7u/193fTnIkyaXTnyPd/Uh3fz/JHVNbAAAA1mAR18y9K8knp+NdmYW7445OZUny+Enlr1/AawOwiVi2xajMYgEjmivMVdV7k7yQ5OPHi5Zp1ll+BrBP8ZzXJ7k+SS666KJ5ugecIcujtgZfSrcGf48ArGbNYa6q9me2Mcrl3X08mB1NcuGSZhckeWI6PlX5Cbr7liS3JMnevXuXDXzAeHwxBQBYrDWFuWlnyhuT/N3ufn5J1aEk/76qPpDZBih7kvxeZjN2e6rq4iR/kNkmKf9ono4DMBaBHgAW63RuTfCJJG9Icn5VHU1yc2a7V740yT1VlST3dffPdveDVfWpJN/MbPnlDd39J9Pz/FySz2Z2a4LbuvvBs/B+AACYk2X3MIZVw1x3v32Z4ltXaP++JO9bpvzuJHefUe8AAABY1qq3JgAAAGDzWcStCYBNxNIYYDNz7STA4ghzAACcEfeUhM1BmAPgtJlVAYDNwzVzAAAAAxLmAAAABmSZJQAAC2MjLlg/ZuYAAAAGZGYOAIB1Y+YOFkeYA06bnQwBADYPYQ42gN9KAgAwL9fMAQAADMjMHGwzlkoCAGwNwhwAAENwmQKcyDJLAACAAZmZg8FYJgkAQCLMAQCwDViiyVZkmSUAAMCAhDkAAIABCXMAAAADEuYAAAAGZAMU2ITsWAkAwGqEOTgLhDEAWJt5/g317y/bjWWWAAAAAxLmAAAABmSZJayRpRwAAGwkM3MAAAADMjMHAMC2t9qKm0cPXrVOPYHTZ2YOAABgQMIcAADAgIQ5AACAAa0a5qrqtqp6qqq+saTs5VV1T1U9PP08byqvqvpwVR2pqq9V1euWPGb/1P7hqtp/dt4OAADA9nA6M3MfTXLFSWUHktzb3XuS3DudJ8mVSfZMf65P8pFkFv6S3Jzk9UkuTXLz8QAIAADAmVs1zHX3F5M8fVLxviS3T8e3J7l6SfnHeua+JOdW1auTvDnJPd39dHc/k+SevDggAgAAcJrWes3cq7r7ySSZfr5yKt+V5PEl7Y5OZacqBwAAYA0WfZ+5WqasVyh/8RNUXZ/ZEs1cdNFFi+sZLGOle8q4nwwAAJvZWmfmvjMtn8z086mp/GiSC5e0uyDJEyuUv0h339Lde7t7786dO9fYPQAAgK1trWHuUJLjO1LuT3LnkvJ3TLtaXpbkuWkZ5meTvKmqzps2PnnTVAYAAMAarLrMsqo+keQNSc6vqqOZ7Up5MMmnquq6JI8luWZqfneStyQ5kuT5JO9Mku5+uqr+VZIvTe3+ZXefvKkKAAAAp2nVMNfdbz9F1eXLtO0kN5zieW5LctsZ9Q4AAIBlrXWZJQAAABtImAMAABiQMAcAADCgRd9nDraMle5BBwAAG83MHAAAwICEOQAAgAEJcwAAAAMS5gAAAAZkAxQAAFjFShujPXrwqnXsCfyAMAcAAHNYbQdsYY+zxTJLAACAAQlzAAAAA7LMkqFZ1gAAwHZlZg4AAGBAZubY0labuQMAgFGZmQMAABiQMAcAADAgYQ4AAGBAwhwAAMCAhDkAAIAB2c0SAADOopV213ZPXOZhZg4AAGBAZuYAAGCDrHZPXDN3rMTMHAAAwICEOQAAgAEJcwAAAAMS5gAAAAYkzAEAAAzIbpZseqvt8gQAANuRmTkAAIABmZljw5l5AwCAMzfXzFxV/fOqerCqvlFVn6iqP1dVF1fV/VX1cFV9sqpeMrV96XR+ZKrfvYg3AAAAsB2tOcxV1a4kP59kb3e/Jsk5Sa5N8v4kH+zuPUmeSXLd9JDrkjzT3T+e5INTOwAAANZg3mvmdiT581W1I8nLkjyZ5I1JPj3V357k6ul433Seqf7yqqo5Xx8AAGBbWnOY6+4/SPJvkjyWWYh7LskDSZ7t7hemZkeT7JqOdyV5fHrsC1P7V5z8vFV1fVUdrqrDx44dW2v3AAAAtrQ1b4BSVedlNtt2cZJnk/yHJFcu07SPP2SFuh8UdN+S5JYk2bt374vqAQBgu1hto7hHD161Tj1hM5pnmeVPJ/l2dx/r7v+X5DNJ/laSc6dll0lyQZInpuOjSS5Mkqn+R5M8PcfrAwAAbFvzhLnHklxWVS+brn27PMk3k3w+ydumNvuT3DkdH5rOM9V/rrvNvAEAAKzBPNfM3Z/ZRiZfTvL16bluSXJjkl+sqiOZXRN36/SQW5O8Yir/xSQH5ug3AADAtjbXTcO7++YkN59U/EiSS5dp+8dJrpnn9RiXG4MDAMBizXtrAgAAADbAXDNzAADAxllp9ZOdLrc+M3MAAAADEuYAAAAGJMwBAAAMSJgDAAAYkDAHAAAwIGEOAABgQMIcAADAgIQ5AACAAQlzAAAAAxLmAAAABrRjozsAAAAs3u4Dd61Y/+jBq9apJ5wtZuYAAAAGJMwBAAAMyDJLTotpegAA2FzMzAEAAAzIzBwLsdrMHQAAsFhm5gAAAAYkzAEAAAxImAMAABiQMAcAADAgYQ4AAGBAwhwAAMCAhDkAAIABCXMAAAADEuYAAAAGJMwBAAAMSJgDAAAYkDAHAAAwoLnCXFWdW1Wfrqrfr6qHqupvVtXLq+qeqnp4+nne1Laq6sNVdaSqvlZVr1vMWwAAANh+dsz5+F9J8l+6+21V9ZIkL0vyniT3dvfBqjqQ5ECSG5NcmWTP9Of1ST4y/QQAANbZ7gN3rVj/6MGr1qknrNWaZ+aq6i8m+TtJbk2S7v5+dz+bZF+S26dmtye5ejrel+RjPXNfknOr6tVr7jkAAMA2Ns/M3I8lOZbkN6rqbyR5IMm7k7yqu59Mku5+sqpeObXfleTxJY8/OpU9OUcfWKDVfjsDAABsHvNcM7cjyeuSfKS7X5vk/2a2pPJUapmyflGjquur6nBVHT527Ngc3QMAANi65glzR5Mc7e77p/NPZxbuvnN8+eT086kl7S9c8vgLkjxx8pN29y3dvbe79+7cuXOO7gEAAGxdaw5z3f2/kjxeVX9tKro8yTeTHEqyfyrbn+TO6fhQkndMu1peluS548sxAQAAODPz7mb5z5J8fNrJ8pEk78wsIH6qqq5L8liSa6a2dyd5S5IjSZ6f2gIAALAGc4W57v5qkr3LVF2+TNtOcsM8rwcAAKyPlTbHc9uCzWGum4YDAACwMYQ5AACAAQlzAAAAAxLmAAAABiTMAQAADEiYAwAAGJAwBwAAMCBhDgAAYEDCHAAAwICEOQAAgAEJcwAAAAPasdEdAAAAxrL7wF0r1j968Kp16sn2ZmYOAABgQMIcAADAgIQ5AACAAblmbhtZbW0zAAAwDjNzAAAAAxLmAAAABiTMAQAADEiYAwAAGJAwBwAAMCBhDgAAYEDCHAAAwICEOQAAgAEJcwAAAAMS5gAAAAa0Y6M7AAAAbC27D9x1yrpHD161jj3Z2szMAQAADEiYAwAAGJAwBwAAMCBhDgAAYEBzh7mqOqeqvlJV/3k6v7iq7q+qh6vqk1X1kqn8pdP5kal+97yvDQAAsF0tYmbu3UkeWnL+/iQf7O49SZ5Jct1Ufl2SZ7r7x5N8cGoHAADAGswV5qrqgiRXJfn16bySvDHJp6cmtye5ejreN51nqr98ag8AAMAZmvc+cx9K8i+S/IXp/BVJnu3uF6bzo0l2Tce7kjyeJN39QlU9N7X/7px9YImV7ukBAABsHWuemauqtyZ5qrsfWFq8TNM+jbqlz3t9VR2uqsPHjh1ba/cAAAC2tHmWWf5Ukp+pqkeT3JHZ8soPJTm3qo7P+F2Q5Inp+GiSC5Nkqv/RJE+f/KTdfUt37+3uvTt37pyjewAAAFvXmsNcd9/U3Rd09+4k1yb5XHf/4ySfT/K2qdn+JHdOx4em80z1n+vuF83MAQAAsLp5r5lbzo1J7qiqf53kK0luncpvTfKbVXUksxm5a8/Ca295rokDAACSBYW57v5Cki9Mx48kuXSZNn+c5JpFvB4AAMB2t4j7zAEAALDOhDkAAIABnY1r5gAAAJa12h4Qjx68ap16Mj5hDgAA2DSEvdNnmSUAAMCAhDkAAIABCXMAAAADEuYAAAAGJMwBAAAMSJgDAAAYkDAHAAAwIGEOAABgQMIcAADAgHZsdAc40Wp3vAcAAEjMzAEAAAxJmAMAABiQMAcAADAgYQ4AAGBAwhwAAMCAhDkAAIABCXMAAAADEuYAAAAGJMwBAAAMSJgDAAAYkDAHAAAwIGEOAABgQMIcAADAgIQ5AACAAQlzAAAAA9qx0R3YjnYfuGujuwAAAAxOmAMAAIax0sTIowevWseebLw1L7Osqgur6vNV9VBVPVhV757KX15V91TVw9PP86byqqoPV9WRqvpaVb1uUW8CAABgu5nnmrkXkvxSd/9EksuS3FBVlyQ5kOTe7t6T5N7pPEmuTLJn+nN9ko/M8doAAADb2prDXHc/2d1fno7/KMlDSXYl2Zfk9qnZ7Umuno73JflYz9yX5NyqevWaew4AALCNLWQ3y6raneS1Se5P8qrufjKZBb4kr5ya7Ury+JKHHZ3KAAAAOENzh7mq+pEkv5XkF7r7D1dqukxZL/N811fV4ao6fOzYsXm7BwAAsCXNFeaq6oczC3If7+7PTMXfOb58cvr51FR+NMmFSx5+QZInTn7O7r6lu/d2996dO3fO0z0AAIAta57dLCvJrUke6u4PLKk6lGT/dLw/yZ1Lyt8x7Wp5WZLnji/HBAAA4MzMc5+5n0ryT5N8vaq+OpW9J8nBJJ+qquuSPJbkmqnu7iRvSXIkyfNJ3jnHawMAAJxgpXvQJVvvPnRrDnPd/d+y/HVwSXL5Mu07yQ1rfT0AAAB+YCG7WQIAALC+hDkAAIABCXMAAAADEuYAAAAGJMwBAAAMSJgDAAAYkDAHAAAwIGEOAABgQGu+aTinttqd5wEAAOZlZg4AAGBAwhwAAMCAhDkAAIABCXMAAAADEuYAAAAGJMwBAAAMSJgDAAAYkDAHAAAwIGEOAABgQMIcAADAgIQ5AACAAQlzAAAAAxLmAAAABiTMAQAADGjHRncAAABgPew+cNcp6x49eNU69mQxzMwBAAAMSJgDAAAYkDAHAAAwINfMrcFKa20BAADWg5k5AACAAQlzAAAAAxLmAAAABrTuYa6qrqiqb1XVkao6sN6vDwAAsBWsa5irqnOS/FqSK5NckuTtVXXJevYBAABgK1jvmblLkxzp7ke6+/tJ7kiyb537AAAAMLz1DnO7kjy+5PzoVAYAAMAZWO/7zNUyZX1Cg6rrk1w/nf6fqvrWWe8VG+n8JN/d6E6wpRhTLJoxxaIZUyyaMbUA9f6N7sEp/ZVTVax3mDua5MIl5xckeWJpg+6+Jckt69kpNk5VHe7uvRvdD7YOY4pFM6ZYNGOKRTOmtq/1Xmb5pSR7quriqnpJkmuTHFrnPgAAAAxvXWfmuvuFqvq5JJ9Nck6S27r7wfXsAwAAwFaw3sss0913J7l7vV+XTcuSWhbNmGLRjCkWzZhi0Yypbaq6e/VWAAAAbCrrfc0cAAAACyDMcVZV1S9X1e9X1deq6j9W1blL6m6qqiNV9a2qevOS8iumsiNVdWBJ+cVVdX9VPVxVn5w20WGbqaprqurBqvrTqtp7Up0xxUKdauzAyarqtqp6qqq+saTs5VV1z/QZc09VnTeVV1V9eBpXX6uq1y15zP6p/cNVtX8j3gsbr6ourKrPV9VD0795757KjSlOIMxxtt2T5DXd/deT/PckNyVJVV2S2W6mP5nkiiT/rqrOqapzkvxakiuTXJLk7VPbJHl/kg92954kzyS5bl3fCZvFN5L8wyRfXFpoTLFoq4wdONlHM/vsWepAknunz5h7p/NkNqb2TH+uT/KRZPZFPcnNSV6f5NIkNx//ss6280KSX+run0hyWZIbps8fY4oTCHOcVd39O939wnR6X2b3FkySfUnu6O7vdfe3kxzJ7EPm0iRHuvuR7v5+kjuS7KuqSvLGJJ+eHn97kqvX632weXT3Q939rWWqjCkWbdmxs8F9YpPq7i8mefqk4n2ZfbYkJ37G7EvysZ65L8m5VfXqJG9Ock93P93dz2T2C9GTAyLbQHc/2d1fno7/KMlDSXbFmOIkwhzr6V1Jfns63pXk8SV1R6eyU5W/IsmzS4Lh8XI4zphi0U41duB0vaq7n0xmX86TvHIqP9PPK7axqtqd5LVJ7o8xxUnW/dYEbD1V9btJ/tIyVe/t7junNu/NbMnAx48/bJn2neV/wdArtGcLOp0xtdzDlikzppiHMcLZcqqxZcxxgqr6kSS/leQXuvsPZ4tKlm+6TJkxtQ0Ic8ytu396pfrpYtu3Jrm8f3AvjKNJLlzS7IIkT0zHy5V/N7MlAzummZSl7dliVhtTp2BMsWgrjSk4Hd+pqld395PTkrenpvJTja2jSd5wUvkX1qGfbEJV9cOZBbmPd/dnpmJjihNYZslZVVVXJLkxyc909/NLqg4lubaqXlpVF2d2we7vJflSkj3TLoMvyWxDi0NTCPx8krdNj9+f5FQzNGxPxhSLtuzY2eA+MZZDmX22JCd+xhxK8o5pB8LLkjw3LZn7bJI3VdV50yYVb5rK2Gam67pvTfJQd39gSZUxxQncNJyzqqqOJHlpkv89Fd3X3T871b03s+voXshs+cBvT+VvSfKhJOckua273zeV/1hmGxC8PMlXkvyT7v7eOr4dNoGq+gdJ/m2SnUmeTfLV7n7zVGdMsVCnGjtwsqr6RGYzIOcn+U5mOwj+pySfSnJRkseSXNPdT09f1H81s40onk/yzu4+PD3Pu5K8Z3ra93X3b6zn+2BzqKq/neS/Jvl6kj+dit+T2XVzxhR/RpgDAAAYkGWWAAAAAxLmAAAABiTMAQAADEiYAwAAGJAwBwAAMCBhDgAAYEDCHAAAwICEOQAAgAH9f/OH7EyNawBEAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 1080x360 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.hist(pKst[:,0],bins=100);\n",
"fig = plt.gcf()\n",
"fig.set_size_inches(15,5)"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"p1_q = lorentz_boost_np(p1, -q[:,0:3]/q[:,3:4])\n"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA3MAAAEwCAYAAAAdGiiRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAcqElEQVR4nO3df6yl9V0n8PendDrbKFlaoXUcYIfdTBOrcauZgAn+AUVbWo2z3ViXmlWsRDQLUbL+UQY3qatWx6zW4tplMworbGopq3U7UVyk/EgXs7RArbUU284WLAMToFIqpoQt9Lt/nGfay3DPj3vvOeee55zXK5ncc7/Pc+59TvJwuO/z/X4/n2qtBQAAgH55yXZfAAAAABsnzAEAAPSQMAcAANBDwhwAAEAPCXMAAAA9JMwBAAD00NgwV1X/pKo+VlV/XVX3V9V/7MbPqqqPVtXnquoDVfWybnxn9/2R7vieNT/rQDf+map646xeFAAAwLKrcX3mqqqSfFNr7R+rakeSu5L8fJJ/n+SDrbUbq+q/Jvnr1to1VfXvknxXa+1nq+qiJG9prf2bqnptkvcnOTvJtyX5cJLXtNaeH/a7Tz311LZnz54pvEwAAID+ue+++77YWjttvWMvHffkNkh7/9h9u6P715K8PsmPdePXJ/mlJNck2d89TpI/SvK7XSDcn+TG1tqzSR6sqiMZBLv/M+x379mzJ/fee++4SwQAAFhKVfV3w45NtGeuqk6qqk8keTzJrUn+b5KnWmvPdaccTbK7e7w7ycNJ0h3/cpJvWTu+znMAAADYgInCXGvt+dba65KcnsFs2revd1r3tYYcGzb+AlV1aVXdW1X3PvHEE5NcHgAAwMrZUDXL1tpTSe5M8r1JTqmq48s0T0/yaPf4aJIzkqQ7/k+TPLl2fJ3nrP0dh1pr+1pr+047bd2loQAAACtvkmqWp1XVKd3jlyf5/iQPJLkjyY90p12c5EPd48Pd9+mO397tuzuc5KKu2uVZSfYm+di0XggAAMAqGVsAJcmuJNdX1UkZhL+bWmt/WlWfTnJjVf1qkr9Kcm13/rVJ/ntX4OTJJBclSWvt/qq6KcmnkzyX5LJRlSwBAAAYbmxrgu20b9++ppolAACwqqrqvtbavvWObWjPHAAAAItBmAMAAOghYQ4AAKCHhDkAAIAeEuYAAAB6aJLWBACwWu749fXHzz8w3+sAgBGEOQBW07DABgA9YZklAABAD5mZA4CtsiwTgG1gZg4AAKCHhDkAAIAesswSACalaAoAC8TMHAAAQA8JcwAAAD0kzAEAAPSQPXMALDf73ABYUmbmAAAAekiYAwAA6CFhDgAAoIeEOQAAgB4S5gAAAHpImAMAAOghYQ4AAKCHhDkAAIAe0jQcAGZlVMPy8w/M7zoAWEpm5gAAAHrIzBwAy2HULBgALCEzcwAAAD1kZg6AfjEDBwBJzMwBAAD0kjAHAADQQ8IcAABADwlzAAAAPSTMAQAA9JAwBwAA0EPCHAAAQA+N7TNXVWckuSHJtyb5WpJDrbWrq+qXkvx0kie6U69qrd3cPedAkkuSPJ/k51prt3TjFya5OslJSX6/tXZwui8HAHpiWL+88w/M9zoA6K1JmoY/l+QXWmsfr6qTk9xXVbd2x367tfaba0+uqtcmuSjJdyT5tiQfrqrXdIffm+QHkhxNck9VHW6tfXoaLwQAAGCVjA1zrbVjSY51j5+uqgeS7B7xlP1JbmytPZvkwao6kuTs7tiR1trnk6SqbuzOFeYAAAA2aJKZua+rqj1JvjvJR5Ocm+TyqvqJJPdmMHv3pQyC3t1rnnY03wh/D58wfs6mrhqA5WCpIQBs2sRhrqq+OckfJ7mitfYPVXVNkl9J0rqvv5Xkp5LUOk9vWb/YSlvn91ya5NIkOfPMMye9PACWybCQtwoEXAAmNFE1y6rakUGQe19r7YNJ0lp7rLX2fGvta0l+L99YSnk0yRlrnn56kkdHjL9Aa+1Qa21fa23faaedttHXAwAAsBLGhrmqqiTXJnmgtfbuNeO71pz2liSf6h4fTnJRVe2sqrOS7E3ysST3JNlbVWdV1csyKJJyeDovAwAAYLVMsszy3CQ/nuRvquoT3dhVSd5WVa/LYKnkQ0l+Jklaa/dX1U0ZFDZ5LsllrbXnk6SqLk9ySwatCa5rrd0/xdcCAACwMiapZnlX1t8Hd/OI57wrybvWGb951PMAAACYzER75gAAAFgswhwAAEAPCXMAAAA9JMwBAAD00MRNwwGAbaSZOAAnMDMHAADQQ8IcAABADwlzAAAAPSTMAQAA9JACKADQZ8MKoySKowAsOTNzAAAAPSTMAQAA9JAwBwAA0EPCHAAAQA8JcwAAAD0kzAEAAPSQMAcAANBD+swBwAjX3vVgnn72q0OPn7xzRy75vrPmeEUAMCDMATB7oxpbz8GoQDYujD397FdzxQWvGXr8Pbd9dsvXBwCbIcwBsPRGBbJr73pwZCA7eeeOWV3W7A0L0ecfmO91ADATwhwA07HNs2+bZYkkAH0lzAHQe5PsawOAZSPMAdB74/a1AcAy0poAAACgh4Q5AACAHrLMEgBmRI86AGZJmAOAGdGjDoBZsswSAACgh8zMAbAxPe0nBwDLxswcAABAD5mZA4AtOHnnjqF73zQrB2CWhDkA2ALVKAHYLpZZAgAA9JAwBwAA0EPCHAAAQA8JcwAAAD00tgBKVZ2R5IYk35rka0kOtdaurqpXJvlAkj1JHkryo621L1VVJbk6yZuTfCXJT7bWPt79rIuT/IfuR/9qa+366b4cAPrq2rsezNPPfnVTz1U1EoBVNEk1y+eS/EJr7eNVdXKS+6rq1iQ/meS21trBqroyyZVJ3pHkTUn2dv/OSXJNknO68PfOJPuStO7nHG6tfWnaLwqA/nn62a/migtes92XsRqGNX4//8B8rwOALRkb5lprx5Ic6x4/XVUPJNmdZH+S87rTrk9yZwZhbn+SG1prLcndVXVKVe3qzr21tfZkknSB8MIk75/i6wGApTBupvLknTu0RQBYcRvqM1dVe5J8d5KPJnl1F/TSWjtWVa/qTtud5OE1TzvajQ0bBwBOMG6mclijcgBWx8QFUKrqm5P8cZIrWmv/MOrUdcbaiPETf8+lVXVvVd37xBNPTHp5AAAAK2Wimbmq2pFBkHtfa+2D3fBjVbWrm5XbleTxbvxokjPWPP30JI924+edMH7nib+rtXYoyaEk2bdv34vCHgAsi5N37hg6w6aoCwDjTFLNspJcm+SB1tq71xw6nOTiJAe7rx9aM355Vd2YQQGUL3eB75Ykv1ZVr+jOe0MSO60BVsQke8BWjT1vAGzFJDNz5yb58SR/U1Wf6MauyiDE3VRVlyT5QpK3dsduzqAtwZEMWhO8PUlaa09W1a8kuac775ePF0MBYPmpVgkA0zVJNcu7sv5+tyS5YJ3zW5LLhvys65Jct5ELBGCbDCtfDwAshIkLoAAAALA4NtSaAABYDKOKpxw/vuE9eZqJA/SKMAcAPTQuqE21D92oJbeCHsC2EeYAVp29cQDQS/bMAQAA9JAwBwAA0EOWWQIwNaMag69iU3AAmCVhDoCp0Ri8H0aF7mSTlTABmDthDoCJTRICWAyjWhecvHPHyNA91UqYAMyMMAfAxMy89YeZNYDlpwAKAABADwlzAAAAPSTMAQAA9JA9cwDL5I5fX3/8/APzvQ4AYObMzAEAAPSQMAcAANBDllkCAJtnaS/AthHmAFbBsD+4AYDeEuYAVsy1dz2Yp5/96rrHTt65Q7NpcvLOHXnPbZ/d9HPdQwDzIcwBrJinn/1qrrjgNese2+wf8CyXrYQx9xDA/CiAAgAA0ENm5gCWzKhllMlgGdyoY6NmVkY9F15AYRSAmRPmAJbMqGWU49jrBAD9IcwBAFMzbnb3f3z09vzlla+f4xUBLC9hDgCYmnGzu++55Zk5XQnA8lMABQAAoIeEOQAAgB4S5gAAAHpImAMAAOghYQ4AAKCHhDkAAIAe0poAoIfOPXh7HnnqxSXer3jpZ3Pyzh3bcEUwmd2nvDx7rvyzocf0oAOYnDAH0EOPPPVMHjr4gy8+cMcn538xsAGjwtqwkAfA+oQ5gAU0bObtuN2nvHyOVwOLYdR/F2b1gFUkzAEsoKEzb9B3d/z6+uPnHxj71FH/XZjVA1aRMAfQR8P+IAYAVsbYMFdV1yX5oSSPt9a+sxv7pSQ/neSJ7rSrWms3d8cOJLkkyfNJfq61dks3fmGSq5OclOT3W2sHp/tSAIA+G1Uc5fhxAL5hkpm5P0jyu0luOGH8t1trv7l2oKpem+SiJN+R5NuSfLiqXtMdfm+SH0hyNMk9VXW4tfbpLVw7ALBE7HkD2JixYa619pGq2jPhz9uf5MbW2rNJHqyqI0nO7o4daa19Pkmq6sbuXGEOAABgE7bSNPzyqvpkVV1XVa/oxnYneXjNOUe7sWHjL1JVl1bVvVV17xNPPLHeKQAAACtvs2HumiT/IsnrkhxL8lvdeK1zbhsx/uLB1g611va11vaddtppm7w8AACA5bapapattceOP66q30vyp923R5OcsebU05M82j0eNg6wcibuI6dqJati1L0+QdsCgFW0qTBXVbtaa8e6b9+S5FPd48NJ/rCq3p1BAZS9ST6Wwczc3qo6K8kjGRRJ+bGtXDhAn+kjBwBs1SStCd6f5Lwkp1bV0STvTHJeVb0ug6WSDyX5mSRprd1fVTdlUNjkuSSXtdae737O5UluyaA1wXWttfun/moAAABWxCTVLN+2zvC1I85/V5J3rTN+c5KbN3R1AAAArGsr1SwBAADYJsIcAABAD22qAAoA442qWPn1apXAVOw+5eXZc+WfjTz+l1e+fo5XBDB7whzAjKhYCVMyrG3BmpYF44LaqKAH0FeWWQIAAPSQmTmA7aYxOACwCcIcwCaN2hOXrLMvTmgDAKZImAPYJHviAIDtJMwBAEtPtUtgGQlzAMDSU+0SWEbCHLDSJtn35tN6WH6jZu68DwCLSpgDVtq4fW8+rYfVMCqseR8AFpU+cwAAAD1kZg5YahtuHwD0x7B2H+cfmO91AGwTYQ5YatoHAADLyjJLAACAHhLmAAAAesgyS4ARxpUrBxbQsL10if10wFIR5gBG0FsKGPWhzvHj3iuA7SDMAQCMMC6onXvwdmEP2BbCHADAFowLapqOA7MizAFM26j9OsDKGbf31qwdsFnCHACwOrah0fiosGbWDtgKrQkAAAB6yMwcsPDOPXh7HnnqmaHHLVMCAFaRMAcsvEeeeiYPHfzBocdHVZLTCw4AWFbCHNB7ZuUAgFUkzAELYdRSSrNrAAAvJswBC2HcUkqAZTSqbcEkz7UyAVabMAewWfrJAVu0lTCmrQGgNQEAAEAPCXMAAAA9JMwBAAD0kD1zAKPYFwerYdh/6+cfmO91AGyAmTkAAIAeGhvmquq6qnq8qj61ZuyVVXVrVX2u+/qKbryq6neq6khVfbKqvmfNcy7uzv9cVV08m5cDAACwGiaZmfuDJBeeMHZlkttaa3uT3NZ9nyRvSrK3+3dpkmuSQfhL8s4k5yQ5O8k7jwdAAAAANm7snrnW2keqas8Jw/uTnNc9vj7JnUne0Y3f0FprSe6uqlOqald37q2ttSeTpKpuzSAgvn/LrwAAYAWNaziuqTgsv80WQHl1a+1YkrTWjlXVq7rx3UkeXnPe0W5s2DgAAJswLqhpKg7Lb9rVLGudsTZi/MU/oOrSDJZo5swzz5zelQEzde7B2/PIU88MPe4TYgCA6dpsmHusqnZ1s3K7kjzejR9Ncsaa805P8mg3ft4J43eu94Nba4eSHEqSffv2rRv4gO0xKrDtPuXleejgD4587rjlQAAATG6zYe5wkouTHOy+fmjN+OVVdWMGxU6+3AW+W5L82pqiJ29IonEL9MwjTz0zMrCNsvCzcvrJAQA9MzbMVdX7M5hVO7WqjmZQlfJgkpuq6pIkX0jy1u70m5O8OcmRJF9J8vYkaa09WVW/kuSe7rxfPl4MBQCgd0Z9AKTRODAnk1SzfNuQQxesc25LctmQn3Ndkus2dHUAAMzEuKXzC7+iAph6ARQAABbAJK0Lhi2dVwkT+kGYAwAYpsf7ac2swfIT5gAAeAENyaEfhDkAAF5AQ3LoB2EOWC09XjIFALDWS7b7AgAAANg4M3OwYpSiBgBYDsIcrJhHnnpmaCnqcw/ePnbDOwAAi0GYA77OrBwAkxhV7dIqD5gfYQ4AgA0ZFdZUuoT5EeaA5aRqJQCw5FSzBAAA6CEzcwAA0zRsZcD5B+Z7HcDSE+YAAJiaUcVRjh9XIAWmQ5iDJTOqj1yivQAAszUuqCmQAtMjzMGSGdVHDoBtZPklMGXCHAAAc2MZJkyPMAcAwNxYhgnTI8zBAhq1780nlgAAJMIcLKRR+97OPXj72OUpAPSIvXQvMGoZpg804YWEOegZ/xNbY9gfQAD01qj/z1mCCS/0ku2+AAAAADZOmAMAAOghYQ4AAKCH7JkDFp+9ccAqGvXet6LFUYAXEuYAAFgKWvuwaoQ5YHGYgQNghFFtC44fH9baRyVMlpEwB9tg1CeHiV5xALAeM2vwQsIcbINRTcEBAGASwhzMgJk3AFgskyzRNPNH3whzMANm3gBgsYwLavbU0Uf6zAEAAPSQMAcAANBDllnCJo3rZQMAMzOslYtm4ptmTx19JMzBJtkXB8DCEfI2zZ46+kiYY6WNm13zCRwAAItqS2Guqh5K8nSS55M811rbV1WvTPKBJHuSPJTkR1trX6qqSnJ1kjcn+UqSn2ytfXwrvx+2atTsmk/gZmjYJ8cAAExsGjNz57fWvrjm+yuT3NZaO1hVV3bfvyPJm5Ls7f6dk+Sa7isspEnWzjOCwAYAMFOzWGa5P8l53ePrk9yZQZjbn+SG1lpLcndVnVJVu1prx2ZwDbBlllgCALDIttqaoCX5i6q6r6ou7cZefTygdV9f1Y3vTvLwmuce7cYAAADYoK3OzJ3bWnu0ql6V5Naq+tsR59Y6Y+1FJw1C4aVJcuaZZ27x8lh1owqcJJZKToXllACLT5XLLRu1/ULRNLbLlsJca+3R7uvjVfUnSc5O8tjx5ZNVtSvJ493pR5Ocsebppyd5dJ2feSjJoSTZt2/fi8IebIT2AQDANIwKa4qmsV02Heaq6puSvKS19nT3+A1JfjnJ4SQXJznYff1Q95TDSS6vqhszKHzyZfvlAADoOw3H2S5bmZl7dZI/GXQcyEuT/GFr7X9V1T1JbqqqS5J8Iclbu/NvzqAtwZEMWhO8fQu/GwAAFsK4oHbuwdst0WQmNh3mWmufT/Iv1xn/+yQXrDPekly22d8HbDN74wBgUyzRZFZm0ZoAAIA+GPVBneIosPCEOXpNtUoAAFaVMEevqVYJAPSZ4ilshTDHtptkds2bGADMmd50czHubxx76hhFmGPbjZtdG1cBCgBgWZm5YxRhjoXnDWqOVKwEgIVi5o5RhDkAAOipUTN3Zu2WnzAHAAA9pYfdahPmmItRRU7sewMAgI0T5pgLLQQAAGC6hDlYRQqdAAD0njAHAMDk9J+DhSHMMRWTNP4GAJaYkLdw9KhbfsIcU2FP3IKynBIAVpYedctPmAMAYHZGfbBo1g62RJhjIpZRAgDAYhHmmIhllAvOckoAYIPsqes/YQ4AAFbQuKB27sHbh4Y9QW8xCHN83aillJZRAgBTpwLmQhsV1kYFvUTYmxdhjq+zlBIAgEmolLkYhDlYNPa/AbDqzNjBRF6y3RcAAADAxpmZWyHaCwAAMA+jKmXaTzc9wtwKsScOAIB5GBXW7KebHmGuZyaZXfNJBwAALD9hrmfGza6N6wcCANBbmykSpmgKS0yYWzJm5XpE1UoAALZANUsAAIAeMjO3YFSc7CmzbAAAExlV6fL4cavNJiPMLRgVJwEApkgD8oUzLqiNqgGRCHtrCXMAAMDCGBfUtDb4BnvmAAAAesjMHKzHHjgAWG7aHLAEhDkAAKA3RhVQWbX9dMLcNhhVsVK1SgCABaWYykIYFdZWbT+dMLcNVKycM0smAYBZEvLYJnMPc1V1YZKrk5yU5PdbawfnfQ2zplfcNhHaAIBFIuTN3ar1sJtrmKuqk5K8N8kPJDma5J6qOtxa+/Q8r2PWzLzNmNAGAMA6Vq2twbxn5s5OcqS19vkkqaobk+xP0qswZ+ZtToQ2AGAZTfNvHLN8K23eYW53kofXfH80yTlzvoYtM/M2ggAGADA/WixsyLJVwpx3mKt1xtoLTqi6NMml3bf/WFWfmflVbUL9xnZfQW+dmuSL230RLCX3FrPgvmJW3FvMygT31lVzuZC++bsktZg5958NOzDvMHc0yRlrvj89yaNrT2itHUpyaJ4XxfxU1b2ttX3bfR0sH/cWs+C+YlbcW8yKe2u1vGTOv++eJHur6qyqelmSi5IcnvM1AAAA9N5cZ+Zaa89V1eVJbsmgNcF1rbX753kNAAAAy2DufeZaazcnuXnev5eFYQkts+LeYhbcV8yKe4tZcW+tkGqtjT8LAACAhTLvPXMAAABMgTDH1FTVf6qqv62qT1bVn1TVKWuOHaiqI1X1map645rxC7uxI1V15Zrxs6rqo1X1uar6QFcwhxVVVW+tqvur6mtVte+EY+4tZmLYPQTDVNV1VfV4VX1qzdgrq+rW7j3n1qp6RTdeVfU73f31yar6njXPubg7/3NVdfF2vBYWR1WdUVV3VNUD3f8Lf74bd28hzDFVtyb5ztbadyX5bJIDSVJVr82gcul3JLkwyX+pqpOq6qQk703ypiSvTfK27twk+Y0kv91a25vkS0kumesrYdF8Ksm/TvKRtYPuLWZlzD0Ew/xBBu9Fa12Z5LbuPee27vtkcG/t7f5dmuSaZPAHepJ3JjknydlJ3nn8j3RW1nNJfqG19u1JvjfJZd37kXsLYY7paa39RWvtue7buzPoI5gk+5Pc2Fp7trX2YJIjGbyJnJ3kSGvt8621/5fkxiT7q6qSvD7JH3XPvz7Jv5rX62DxtNYeaK19Zp1D7i1mZd17aJuviQXXWvtIkidPGN6fwXtN8sL3nP1JbmgDdyc5pap2JXljkltba0+21r6UwQelJwZEVkhr7Vhr7ePd46eTPJBkd9xbRJhjdn4qyZ93j3cneXjNsaPd2LDxb0ny1JpgeHwcTuTeYlaG3UOwUa9urR1LBn+UJ3lVN77R9y9IVe1J8t1JPhr3FtmG1gT0W1V9OMm3rnPoF1trH+rO+cUMlgS87/jT1jm/Zf0PE9qI81lik9xb6z1tnTH3FtPgXmHWht1j7j3WVVXfnOSPk1zRWvuHwWKT9U9dZ8y9taSEOTaktfb9o453m2l/KMkF7Rt9L44mOWPNaacnebR7vN74FzNYEvDSbgZl7fksqXH31hDuLWZl1L0FG/FYVe1qrR3rlro93o0Pu8eOJjnvhPE753CdLLCq2pFBkHtfa+2D3bB7C8ssmZ6qujDJO5L8cGvtK2sOHU5yUVXtrKqzMtiQ+7Ek9yTZ21UXfFkGhSwOdyHwjiQ/0j3/4iTDZmZYbe4tZmXde2ibr4l+OpzBe03ywvecw0l+oqs8+L1JvtwtlbslyRuq6hVdcYo3dGOsqG6/97VJHmitvXvNIfcWmoYzPVV1JMnOJH/fDd3dWvvZ7tgvZrCP7rkMlgf8eTf+5iTvSXJSkutaa+/qxv95BgUHXpnkr5L829bas3N8OSyQqnpLkv+c5LQkTyX5RGvtjd0x9xYzMewegmGq6v0ZzHycmuSxDCoH/s8kNyU5M8kXkry1tfZk9wf672ZQgOIrSd7eWru3+zk/leSq7se+q7X23+b5OlgsVfV9Sf53kr9J8rVu+KoM9s25t1acMAcAANBDllkCAAD0kDAHAADQQ8IcAABADwlzAAAAPSTMAQAA9JAwBwAA0EPCHAAAQA8JcwAAAD30/wHP8w7zGQbSEgAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 1080x360 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"#plt.hist(p1_B[:,0],alpha=0.3,bins=100);\n",
"plt.hist(p1[:,0],histtype=\"step\",bins=100);\n",
"plt.hist(p1_q[:,0],alpha=0.5,bins=100);\n",
"fig = plt.gcf()\n",
"fig.set_size_inches(15,5)"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"105.0"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.sqrt(scalar_product_4_np(p1_q,p1_q)).mean()"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [],
"source": [
"z=np.array([[0.,0.,1.,0] for i in range(n_events)])"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAV1UlEQVR4nO3df+xdd33f8ecbe0kZsCXU7po6MTad6Wa1CKpvQzUmIDRQp5VspqadU7ElXTYrXb0uZZ3qCIimVBWBSoNOtVRcyOhagSGZOr5djSxKgjaqJrPZaKgdOTFOt3zjrDGQ0FbtTAzv/XHP1zm+vvf7Pff7Pef++NznQ7J8zzmfe+77nnvv63zO5557vpGZSJJm30smXYAkqR0GuiQVwkCXpEIY6JJUCANdkgqxcVIPvGnTpty2bdukHl6SZtIXv/jFr2bm5kHLJhbo27Zt4/jx45N6eEmaSRHxv4ctc8hFkgphoEtSIQx0SSqEgS5JhTDQJakQjQI9InZFxKmIOB0RB4a0+amIOBkRJyLi4+2WKUlazaqnLUbEBuAg8DZgCTgWEYuZebLWZgdwF/DGzHwuIr6rq4IlSYM16aFfD5zOzDOZ+U3gMLCnr82/AA5m5nMAmflsu2VKklbTJNC3AE/VppeqeXWvAV4TEX8YEQ9HxK5BK4qIfRFxPCKOnzt3bm0VS5IGavJL0Rgwr/+vYmwEdgBvAa4F/ntEfH9mPn/JnTIPAYcAFhYW/Msa6s5D73vx9g13ze5jSCNo0kNfAq6rTV8LnB3Q5tOZ+UJmPgmcohfwkqQxadJDPwbsiIjtwNPAXuCn+9r8F+AW4GMRsYneEMyZNguV5ppHA2pg1UDPzAsRsR84CmwA7svMExFxD3A8MxerZW+PiJPAt4B/m5lf67JwaaoYuJoCja62mJlHgCN98+6u3U7gXdU/zTvDTZqIiV0+VyqWOzRNiIGu2WNgSgMZ6FKJ3OnNJQNd42PISJ0y0DUb6juDeTGPz1nrYqBrvqzlKMEjC80IA12C7kLbnYHGyEDXZHQRdG0NUTjUoRlloKt8BrTmhIGu+WXQqzAGujQuTYaZ3MloHQx0TZ4hJrXCQNfqPFOjW+vZofnaqMZAlyahi6MSj3TmnoGuwYaFw7z3CA1NTTEDXZ16470P8vTzfw3Ae1/2JLf/w+3tPoAB245531EXwkDX2jUI06ef/2v+9N4fB+BD73mg64qkVtU7JFuueil/eOCtE65oZQa6NE9G7Ym32XOfwaOAeodk24Hfn3A1qzPQSzXswzPBD9UrrvwbfOhzj1+83frwizTnZj7QZ+2QaJ7VA3w52LU+H/3Ck/zF+RcAd5IqINBn7ZBIKzOgRvMX51/gzh95DbCGnaRfKBdn5gNdZVlXQM2LUYO41r7RDrOwoJ+no3gDXTNtHnv0/c95FK326P1ic+rMZaDP0x57JB30zO7c+AA89Gjr611WWo++yQ6q/pz1Ij/XcxroY91jz+CpWqUZtRdfb9/0Pm2Z1R3UR39538Vtdv/L/8lEwnSeeuLDzEWg1/fc0Nt7j3KfWdzb10Pp/kcevFj/sPlrWe+sDHGMGpL9PeAm9+liu8zSaZ6XbOOjL37Who7ff2Hfi89nTjo648iURoEeEbuAXwM2AB/JzHv7lt8G/CrwdDXr1zPzIy3WuS71Pfda7vPGex+8uMevvxDTHPrDPmBDP3hrWe+U9SAnGYD17fLRLzw5cNuMWtOkTvMctnPqP3IZtUMwbe+dYZ/f/vltGccRxKqBHhEbgIPA24Al4FhELGbmyb6mn8zM/R3UOHH1N239hZiVQ7wtV730Yn3vfdmLX6K1Ob7ddc+9P0yW1b8UnJbz3Ic990nVNOqXqMOC97IjlzV0CAbWN8bhmvpnYctVLx34+V1LB3BaNOmhXw+czswzABFxGNgD9Ae6RtC0d7Dam7tJ+0vmdfQFZde9L78IHE1/iHex7d77sk9fvD7PqGfb1DU6ahzxu6j+4F427qPornr7wzQJ9C3AU7XpJeANA9r9RES8CXgc+IXMfKq/QUTsA/YBbN26dfRqC7LSkM4ovf5h66n3xNd79sp6TpMrQX0op3/+uIZ1+oeTVtPWDnClxx323PuDfhJj/+MM7v7v6Orqn+dxaBLoMWBe9k3/HvCJzDwfEXcAvwVctkUz8xBwCGBhYaF/HRPRvyfv8o0wbG99aQ/6ffDQI9XEay+dv2xID6WrnviwcBg1ZIaZ1A6j6TBR10MoTZ5/k1Bcz+sxrIa1hHEXQ1/r/TK/S9M0RNMk0JeA62rT1wJn6w0y82u1yd8E3r/+0oZr8zBm2Pj4qOo7hrr6NcDX88LX39AfOvri44zjMG6Ytnpe4xxO6Q+9QV9kjvsopK3nv57XY1qGtIZ1sNb7Zf68aBLox4AdEbGd3lkse4GfrjeIiGsy85lqcjfwWKtV9mkSjGsJ/WHjbk0M6zG0dQ3w+hv6zhuGPPcRh1b6hxKmeQihLaP2wts0bachjuM1GPUx6p+jD73nZ2pHqy8a51H1rFk10DPzQkTsB47SO23xvsw8ERH3AMczcxH4+YjYDVwAvg7c1mHNjaylNzxvb4z+QBnnWRhdDyFMo/pznuQRwaB6Zukxhh1VX3LW1pyc296v0XnomXkEONI37+7a7buA+dyCHer6Z/PTbtI92C6V/Nxm0XqOzqfJXPxSdJLqvcwtVw06OajZfYf24gq7Ml5bSuvdq1ulHJ0b6CMadU9e74ndeUPf2SzLhhwe2otbO7edujTu88ubKirQx3HY1MmefB3Xt55G89I7npfn2YV1bbtL3v+vHdqsS9N0qmJdUYFeymHTpLR1Fsa89I7n5Xl2Yei2m+LOyiyMsxcV6FqfabkWijSNZqHDaKBrIIcTpNljoI9qPX+wYooPJ/s5nCDNHgO9CzMU3NIsG3Zp6KEK/wtiBrqkqdRk2G8cl4aeJQa6pKnksN/oDHRJs8vhzUsY6G3xjSVpwgx0SWWbo87WSyZdgCSpHfbQJZVnjnrldfbQJakQ9tDrCv/RgaSyGejDNAn3OT2skzSdDHRJqpvhI3XH0CWpEPbQJWmYGeut20OXpELYQ5ekUU1pz91Ab8KzWaTyFPi5dshFkgrRqIceEbuAXwM2AB/JzHuHtLsZuB/4ocw83lqVkjRpM9CjXzXQI2IDcBB4G7AEHIuIxcw82dfuFcDPA490UWhnZuBFkqQmmvTQrwdOZ+YZgIg4DOwBTva1+2XgA8AvtlphFwxxSW3pz5MJfknaZAx9C/BUbXqpmndRRLweuC4z/+tKK4qIfRFxPCKOnzt3buRiJUnDNQn0GDAvLy6MeAnwQeDfrLaizDyUmQuZubB58+bmVUqSVtVkyGUJuK42fS1wtjb9CuD7gc9HBMB3A4sRsXuiX4xO6XmiktSVJj30Y8COiNgeEVcAe4HF5YWZ+Y3M3JSZ2zJzG/AwMNkwl6Q5tGoPPTMvRMR+4Ci90xbvy8wTEXEPcDwzF1dewxTwS1BJc6DReeiZeQQ40jfv7iFt37L+siRJo/KXopJUCANdkgox+xfn8mwWSQLsoUtSMQx0SSrETA653LnxAXjo0UmXIUnNXXL69Gs7eYiZDHRJmloT/F7PIRdJKoSBLkmFKGvIxZ/4S5omY84ke+iSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEI0CvSI2BURpyLidEQcGLD8joj4ckR8KSK+EBE72y9VkrSSVQM9IjYAB4GbgJ3ALQMC++OZ+QOZ+TrgA8C/b71SSdKKmvTQrwdOZ+aZzPwmcBjYU2+QmX9em3wZkO2VKElqosnfFN0CPFWbXgLe0N8oIn4OeBdwBfDWQSuKiH3APoCtW7eOWqskaQVNeugxYN5lPfDMPJiZ3wv8EvCeQSvKzEOZuZCZC5s3bx6tUknSipoE+hJwXW36WuDsCu0PA+9YT1GSpNE1CfRjwI6I2B4RVwB7gcV6g4jYUZv8ceCJ9kqUJDWx6hh6Zl6IiP3AUWADcF9mnoiIe4DjmbkI7I+IG4EXgOeAW7ssWpJ0uSZfipKZR4AjffPurt3+1y3XJUkakb8UlaRCGOiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhTDQJakQBrokFcJAl6RCNAr0iNgVEaci4nREHBiw/F0RcTIiHo2Iz0XEq9ovVZK0klUDPSI2AAeBm4CdwC0RsbOv2f8CFjLztcADwAfaLlSStLImPfTrgdOZeSYzvwkcBvbUG2TmQ5n5V9Xkw8C17ZYpSVpNk0DfAjxVm16q5g1zO/CZ9RQlSRrdxgZtYsC8HNgw4p3AAvDmIcv3AfsAtm7d2rBESVITTXroS8B1telrgbP9jSLiRuDdwO7MPD9oRZl5KDMXMnNh8+bNa6lXkjREk0A/BuyIiO0RcQWwF1isN4iI1wMfphfmz7ZfpiRpNasGemZeAPYDR4HHgE9l5omIuCcidlfNfhV4OXB/RHwpIhaHrE6S1JEmY+hk5hHgSN+8u2u3b2y5LknSiPylqCQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEI0CvSI2BURpyLidEQcGLD8TRHxPyPiQkTc3H6ZkqTVrBroEbEBOAjcBOwEbomInX3N/g9wG/DxtguUJDWzsUGb64HTmXkGICIOA3uAk8sNMvNPq2Xf7qBGSVIDTYZctgBP1aaXqnkji4h9EXE8Io6fO3duLauQJA3RJNBjwLxcy4Nl5qHMXMjMhc2bN69lFZKkIZoE+hJwXW36WuBsN+VIktaqSaAfA3ZExPaIuALYCyx2W5YkaVSrBnpmXgD2A0eBx4BPZeaJiLgnInYDRMQPRcQS8JPAhyPiRJdFS5Iu1+QsFzLzCHCkb97dtdvH6A3FSJImxF+KSlIhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKkSjQI+IXRFxKiJOR8SBAcuvjIhPVssfiYhtbRcqSVrZqoEeERuAg8BNwE7glojY2dfsduC5zPy7wAeB97ddqCRpZU166NcDpzPzTGZ+EzgM7Olrswf4rer2A8CPRES0V6YkaTUbG7TZAjxVm14C3jCsTWZeiIhvAN8JfLXeKCL2Afuqyb+MiFNrKRrY9Au/cum6p8QmsK4RWNfoprU26xrNpnj/mut61bAFTQJ9UE8719CGzDwEHGrwmCsXFHE8MxfWu562WddorGt001qbdY2mq7qaDLksAdfVpq8Fzg5rExEbgb8NfL2NAiVJzTQJ9GPAjojYHhFXAHuBxb42i8Ct1e2bgQcz87IeuiSpO6sOuVRj4vuBo8AG4L7MPBER9wDHM3MR+Cjw2xFxml7PfG+XRdPCsE1HrGs01jW6aa3NukbTSV1hR1qSyuAvRSWpEAa6JBViagM9In4yIk5ExLcjYujpPcMuS1B9iftIRDxRXZbgipbqemVEfLZa72cj4uoBbW6IiC/V/v2/iHhHtexjEfFkbdnrxlVX1e5btcderM2f5PZ6XUT8UfV6PxoR/7i2rNXttZ7LWETEXdX8UxHxo+upYw11vSsiTlbb53MR8arasoGv6Zjqui0iztUe/5/Xlt1ave5PRMSt/fftuK4P1mp6PCKery3rcnvdFxHPRsSfDFkeEfEfqrofjYgfrC1b//bKzKn8B/x94PuAzwMLQ9psAL4CvBq4AvhjYGe17FPA3ur2bwA/21JdHwAOVLcPAO9fpf0r6X1R/Der6Y8BN3ewvRrVBfzlkPkT217Aa4Ad1e3vAZ4Brmp7e630fqm1+ZfAb1S39wKfrG7vrNpfCWyv1rNhjHXdUHsP/exyXSu9pmOq6zbg1wfc95XAmer/q6vbV4+rrr72/4reyRydbq9q3W8CfhD4kyHLfwz4DL3f7vww8Eib22tqe+iZ+VhmrvZL0oGXJYiIAN5K7zIE0LsswTtaKq1+mYMm670Z+Exm/lVLjz/MqHVdNOntlZmPZ+YT1e2zwLPA5pYev249l7HYAxzOzPOZ+SRwulrfWOrKzIdq76GH6f0epGtNttcwPwp8NjO/npnPAZ8Fdk2orluAT7T02CvKzP/Gyr/B2QP8p+x5GLgqIq6hpe01tYHe0KDLEmyhd9mB5zPzQt/8NvydzHwGoPr/u1Zpv5fL30y/Uh1ufTAirhxzXd8REccj4uHlYSCmaHtFxPX0el1fqc1ua3sNe78MbFNtj+XLWDS5b5d11d1Or5e3bNBrOs66fqJ6fR6IiOUfIU7F9qqGprYDD9Zmd7W9mhhWeyvbq8lP/zsTEX8AfPeARe/OzE83WcWAebnC/HXX1XQd1XquAX6A3jn8y+4C/i+90DoE/BJwzxjr2pqZZyPi1cCDEfFl4M8HtJvU9vpt4NbM/HY1e83ba9BDDJjX9DIW63pPraLxuiPincAC8Oba7Mte08z8yqD7d1DX7wGfyMzzEXEHvaObtza8b5d1LdsLPJCZ36rN62p7NdHp+2uigZ6ZN65zFcMuS/BVeocyG6te1qDLFayproj4s4i4JjOfqQLo2RVW9VPA72bmC7V1P1PdPB8R/xH4xXHWVQ1pkJlnIuLzwOuB/8yEt1dE/C3g94H3VIeiy+te8/YaYJTLWCzFpZexaHLfLusiIm6kt5N8c2aeX54/5DVtI6BWrSszv1ab/E1evHT2EvCWvvt+voWaGtVVsxf4ufqMDrdXE8Nqb2V7zfqQy8DLEmTvW4aH6I1fQ++yBE16/E3UL3Ow2novG7urQm153PodwMBvw7uoKyKuXh6yiIhNwBuBk5PeXtVr97v0xhbv71vW5vZaz2UsFoG90TsLZjuwA/gf66hlpLoi4vXAh4Hdmflsbf7A13SMdV1Tm9wNPFbdPgq8varvauDtXHqk2mldVW3fR+8Lxj+qzetyezWxCPzT6myXHwa+UXVa2tleXX3bu95/wD+it9c6D/wZcLSa/z3AkVq7HwMep7eHfXdt/qvpfeBOA/cDV7ZU13cCnwOeqP5/ZTV/AfhIrd024GngJX33fxD4Mr1g+h3g5eOqC/gH1WP/cfX/7dOwvYB3Ai8AX6r9e10X22vQ+4XeEM7u6vZ3VM//dLU9Xl2777ur+50Cbmr5/b5aXX9QfQ6Wt8/iaq/pmOp6H3CievyHgL9Xu+8/q7bjaeBnxllXNf3vgHv77tf19voEvbO0XqCXX7cDdwB3VMuD3h8M+kr1+Au1+657e/nTf0kqxKwPuUiSKga6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKsT/B/eyn3hQzPzxAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"cos_theta_l_np=get_costheta_l_np(p1_q,z)\n",
"plt.hist(cos_theta_l_np,histtype=\"step\", density=True, bins=100);\n",
"cos_theta_l_np_B=get_costheta_l_np(p1,z)\n",
"plt.hist(cos_theta_l_np_B,alpha=0.5, density=True,bins=100);"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2.2160899504353184e-15"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"(p1+p2+pKst)[:,0].mean()"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2.1290421509112535e-15"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"(p1+p2+pKst)[:,1].mean()"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-5.32265759245476e-16"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"(p1+p2+pKst)[:,2].mean()"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5280.0"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"(p1+p2+pKst)[:,3].mean()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"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.6.8"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Davide Lancierini