{
 "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",
    "from utils.kin_utils_2 import *\n",
    "from utils.BtoKll_utils import *\n",
    "\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": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "mmu_mass = 105.0\n",
    "mB_mass = 5280.0\n",
    "mKst_mass = 892.0\n",
    "mgamma_mass = 0\n",
    "mBbar_mass = mB_mass-mgamma_mass\n",
    "\n",
    "minq=2*mmu_mass\n",
    "maxq=(mBbar_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",
    "def mBbar(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",
    "        mBbar_mass = tfp.distributions.Uniform(low=min_mass, high=max_mass).sample()\n",
    "        \n",
    "        return mBbar_mass\n",
    "    \n",
    "q=Particle('q', modq).set_children(el1, el2)\n",
    "Kst = Particle('Kst', mKst_mass)\n",
    "\n",
    "Bbar = Particle('Bbar', mBbar).set_children(Kst, q)\n",
    "k = Particle('gamma', mgamma_mass)\n",
    "B = Particle('B', mB_mass).set_children(Bbar, k)\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": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_events=100000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "105.00013069979116 5273.578115490983\n"
     ]
    }
   ],
   "source": [
    "Kst=particles_np['Kst']\n",
    "l1 =particles_np['l1']\n",
    "l2 =particles_np['l2']\n",
    "q_np=particles_np['q']\n",
    "Bbar_np = particles_np['Bbar']\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": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "44100.1097878929"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q2_array.min()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "Bbar_array= scalar_product_4_np(Bbar_np,Bbar_np)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "19.19822678214951"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(q2_array/1e6).max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO2de7BlaVnen3dfz62v9NyYaUACMRWt4uIUoCQWIYQATg1JxIgVFdSqKSwppZQyolVIUeUfxIpKnCjViiWjBDAIOloQncRMoUmgbMaZ4TIgA4wzPd10T9/O/Zx9e/PH3n2+53v3Xuvs03269+ndz6/qVK+917e+9e21d6/1fO/3XszdIYQQ4vqnMukBCCGE2B10QxdCiClBN3QhhJgSdEMXQogpQTd0IYSYEmqTOnHDmj6D+UmdXgghrkuWceGsu980at/EbugzmMfLK/9qUqcfiVVsa9t7cucUO8B7+WurXN6+IsY9xkom3eOeS+xp/qd//B+K9snkIoQQU8LEFLqZoVIfnL5MfVSrabvbpQ7Cs4iPK1MpDClyAACpcrscNVMyJlb8Fs97NbkcNThuf6HP7DPy9wbA+bvLurPidnSusmtWNpsae9ZVdG3G/S2hOv6+7Fy0r+w3Pfa5yrjc43ZGvM6X83vf8zPkPTrbkUIXQogpYWIKHWZArX96q4z5XKmVDJf76NHTsyS1gffibCApCavUCtsVjXeoXbWeXpDyjMfnqtQK2xWea9z0DXYVvm6n8RqNo5qPPVNprL6CesuUfZlKo+OMhWc8plIwviGok8uZ7e3Erk2/rWy8Q31Uadd4M7yyduOq3qL+y5R3Pjsr7q9oBjbopHgMY840iz7/5c4aCq/ZbvwurgJS6EIIMSVMVKFbfaBgq7v8XBmzv6FndMHT2KICJhWdKe+ZZt6uy4qjeEzG/RWtGQDZLCQbe7XENsp9lM1wuF0ca8EMYqgtzxpiH3E2NCDa1rPrVGqWJvXabqfj68F2T+flfUM2fRtvdpZB3+/Qd8DXqez3U+XrGRUrqU2MqWx5HHH2YwUzo/i799Gfa2iGU3TeElUa11by7i7jPmBxJjh6dpW9H48rWAfaEXvEpi6FLoQQU4Ju6EIIMSVM0OQC2CUTwLgmlzjdL4KnuHG6y2aHVjvfV6f+x1yozab7cSpco9cNWiDtdPJ2RdO8eq24HX0OL5vishmok5sZPJsyN9L70RzB0+Sy76BesNgX4Wl8/IxMm65TyVQ9M9NYySIr/xZiu3odo7BopuoUjOly6wpw/7EP/h7YxBTNQLwwXWbSsIKF9Pj/z0vMMcy4LsVl5h1iyCyS70zDG+0BO/rcRe8X/J8Z2z12j5hYIlLoQggxJYyl0M3sCQDLALoAOu5+Z9hvAN4P4A0A1gC81d0f2qbXpAxYpZQtppWp5iJVXguKhVXq7Mx4fZQtTs5QH1F5s+rjzxHVIKtFUvxDapgVP4+hHWYatFjnm63UrpGf19i/jK9LXEyrjVavQ5Sp95L+C/tgVRYVEasl/lzd2I5mIU2arZT1x8OJAVLxOx41hjiOeF3498TbZSqfZzUl+7IZYxh7kXofWpjm/qx4YT9bcG6UzPCysVIfRdcS2Gahv3jXbnM9qHJmJyaXf+HuZwv2vR7ACwd/Lwfw24N/hRBCXCN2y4b+RgD3eb9A6WfN7KCZ3ebupwqPoMCiIVsxw2quwkEnUTWzLbLkSUpKYqgdK6RMATbydjwmtsOXuC2ObW9l9V7WrsQezLOSzIY+5hqEVYLaLHJNBIoVZlRY/BWXjT2z15NLaBRzTXbPG0/1sXIcWu9o0GyAZjVDswkeX9naT9m1HtOun527zEWST5v9vott8pm6HjPQbVzGTfsw9Dl6JbNidkEscc30dsH3P2Yw0vWgwssY14buAP7SzD5vZveM2H87gKfo9YnBexlmdo+ZHTez463e+s5HK4QQopBxFfor3f2kmd0M4AEz+4q7f4b2j3qMD8lLdz8G4BgAHGje4ltKkm3oUREU2cN70QOEbZYl9rfMJl8v3lemTHhfiRdJNqNo0rmiiuD+WB1GWzsrQrabl80gSsLss9lFgZfH0HFxTQJ0HPdXFpzE1yxei6IAp+gBUbSeEq9Fi9YQ+DPG74qwZrHKLdwXvWH4WkT7Ou/j7zTY/7N1g5J22e8n++5R2C5XucXrVqUeRLSvTPGPa7vPPqMX//9jFW4exzQ6zcBQwNB1rsSLGEuhu/vJwb9nAHwSwMtCkxMAjtLrOwCc3I0BCiGEGI9tFbqZzQOouPvyYPu1AN4bmt0P4O1m9lH0F0MXS+3n/Y63VI2TgrFgo/UsFJwURrtEDRcdjxDGP66NuoxqsYLJ1CYr0bJZSKZeg/dKkX9wLSYyotccCj+UjKxA9UXbsBXMDIBiu/mQ7TmsL2yNr8TXnvuemy0+b6fEUyTzbClR1EyZ90WRzXvIE6pEYVL/1qDPNeShUzC7GtfWHij0RCnxcsm+nyL7dOi7dB1j3FlwPI5VOXvhdNqjmvf3lXnbTCnjmFxuAfDJQb6RGoD/5u7/w8zeBgDu/gEAn0LfZfFx9N0Wf+zqDFcIIUQR297Q3f0bAF404v0P0LYD+Kkdn33EE9qbub3R2e5HT/fYDqTYMxVey9Vmj/23o310bTP1z7OGqProdTYDqMfPky6vdWjs1WD32yRFUyvxhuFzzZa0YxXULugbCGsX7Cc/ptcRkKuxZoltnCmb1axvjN4XvVKKlFlJ0q2x4xNYocYZSdlaQ9F5467maDU75Gufpf4tSVqW2Y3Jlh1nIXTdeXRDc1H+TsdMVuW8VlHir57vKImMHUq6RZ+RvvsbxTY+LooUFUKIKUE3dCGEmBImXLGoP3f0meLQ/6HFz0vE6RqZMZyna518CmZkFommD8wlkwGbeuJiIptZKqvJRNBr5KkEjN2m4rmI3v60MGYtXjwtcyXk88RpZ8FzOkxPfYbMSivpc3gwRxi72cUpbtHCYJlJgxfuYjteFC5bqF2jMS3M03lKEp9lQVslrn/cx7hVsuLvsWxhtcBN18ZdjA3n8o1kKiw1lxUtmEYzVZaKg01Wxa6o2YIrmV/6b5QETxURvp/MBfEGXOwcFyl0IYSYEiaaPtcHLnXOlWSCEit64nhY7LQNcmvKUpyG0GBW8jHgoUKLncbt8svE52J1PeQiOWZiscoyqeNZDowpcUGjmUtvIcwM1miBKnMrDIuxG5y2IJ3XQtBNbyF9xspKiPAtSNw19P1kCcj4+wkLaLxovbSSdtSD2yKr8jLFlgXacNBSrDBUkFo2Kt6izzEUcDWmWySr4bjgGhdkLxF/ZzHJ3ICYSCxbWOXrEj8jK/luiboumtUEMjfDksVyqfArRwpdCCGmhMkpdNiW0siUbfTAI/WeuS0GteBzSd1UyPMtilxvUGGIOBsgpdzbl1RPZTXYBAsCnHw2qE0O6ol2bm63b7TCqqzlaoZVb/dAUqzVpY2sHStHn6e+o+21WRDQFZR8ZZOUYlCimRK3ktkPkc0aYjV2VtELc2l7MwY0jU5W5c089N9YEV6OLTeuGfQKUhOUBPQMwbOLNZrxxIIr/Psss+Vzugh2TYyuntUCm39QzUX28CHFX1Z0g2BVXqrCb3CXw91ACl0IIaaEiSl0N0PvklouC5Nmmzc9f2LunkorPd17syFBE/fHj7Co8klRVyjYp304t9/WSLFnNt/okcOCg5Rsrxls8q3R9sL2kfnsdXWd1wlS551Dc1k7Vs3VJfaACN4rrHp5BhHaeVZRoCSdLFepjzZVTuTEnz8qefac4DHFQDIb7a2UrQuE8+bFFaJnR0Fq5qi8s2Rs9BkrMYVBQdoHIA+e4mCseC2it8gl4uyCx8QeL3F9glMOsHoPNvQi+/VQoFKR5030pmJVXhKqL64cKXQhhJgSJmhDT3SpWEGvkT9jqhvkzVEvfv70qA8+Zsg/uMR2ysqZzxVD/7ldhdR1L6RuZeVYW0yqzOIYMh/6tM2KvH+CgllITGHQprWGgjUIAOjNJ2VXXUy2XC+7ZnFWU5CaYcirh2YUmV27WmJ7pZmGdWKYOPWHErXJ6n2V1OtMmMW1C8Ldh1IYk4ouS3XMir+siEe0mzNFaZXDOkHWjlT5ULk8Lk1YHf1+7AN03YeLU9O+Ajv54ECIa4MUuhBCTAmTU+gVoDcz8EOn9K/Vzfzp3qEkVFVWnsHuzjZ0Psai+KD+uR0A1FdGe0EMRXmSguuQd00ce4/t5rPkUdKO7UYn8erOhURlNI4qRUrG/tiTh2cX0XZfIXtzl3zZPcyEuF1cn2D7Pyv0oc/INvBasY4oTJfcze3JTjOywmhiAMbfCavyaENnlct+3XF9h+3f3C6q106Bd00ZQ33QGLlIxma4FqSAuYhHtHl7O3hDbTUsLnySqfLoJcUeMPIb3xNIoQshxJSgG7oQQkwJkzO5+LA5BAA68/mQeIEzC8ePabnJfFKJC2hEXHRlWgdGh7HXV0KAD0/DObakWRxMw0FMsR3vMzLnVNr5BcoWO8k00VsI14zcKtsHklmgtpYvwPUoORe7aaJnxe2iWyB9J5lZKbpIrtM0npKxDSUW4wXiknQJ2YIcm5ji+GhxMmsXFwKzFAl03rhoGWuWbnVYUoWqrC3/lqLpqKh6VXBHzBaF2UwTzUV8PctMKexmyIudcjnc80ihCyHElDC2QjezKoDjAJ5297vCvrcC+FUATw/eutfdf7e8Q6A3WAxtzxer69Y+chEkNTcUWNSmoBZq11kIaphzNYVanPVVWlzihc+4eLqUlEprjhRbUESsqHvkTje0oNthVU6LtnP511NbocVJmmnERUxrs/slp7cN7o20zS6CMSUCw2odyBN5VThp2Uw+9h4FP2WpFGI9S/r8FV6cC+fNApXqJQukBSH5HoO72qNngsMupnRtODArqmFW8jHJFgf/jJnArTTlcFzgHTC0UMkz3G5x8FRWEUiLndcVOzG5/AyAxwDsL9j/MXd/+5UPSQghxOUw1g3dzO4A8H0AfgXAz+72IIyjhoNqZlXO+6JCz1Qvex9GUyn3H9RXh5Jr1VdJsQUF3D6Q1Ber680DuXLiU/EMorZR4mbHRdFjQNMMzWR4bSHG3NA6BK9BdIcUPyvlYvdGtsl3FnIbcoVyKbAqt1BYJFOH7HLYjS6XdA3Znt4IP1Uab6aug0uks12fiocMpVvg3w/30Q4/NFazfEwM4uFAo2iHL1Li0W2R23FNzaD4OYw/q0sazssuiM6zH9nGp4Zxbei/AeDnUZagG/h+M3vUzD5uZkdHNTCze8zsuJkdb7VXdzpWIYQQJWyr0M3sLgBn3P3zZvaqgmZ/BuAj7r5pZm8D8CEAr46N3P0YgGMAsO/AHX4pXJ/Vdgziac+k1xXOJRWcDWrrSbG1yW5uQb2ybTyqfH7ZmaOApvUYJEMl6DbJ/h3s/5bNFNKL6E1TW6P0ASVBN+39aUwNPlUUkTSjaC4mdR0VOtueOzSm6kauNjv7KEXAelB97G3DKREs2rLZK6UkLL7gmCH7Pyti7iNOV3jGQ7OGWGYv65vt8GXj4xQBMXEVq2hOwBX3sSqPfRSkpI1pcTM7d0E4PiDb+I3AOAr9lQDuNrMnAHwUwKvN7A+5gbufc/dLKz2/A+C7dnWUQgghtmVbhe7u7wLwLgAYKPR3uvsPcxszu83dTw1e3o3+4um2bHm5LJDibUcllrY7pNZrG8F3lsPESTl1Z4IHSEkmWGuR2qZ2rNaBXHm3yQunthH8xkmxZ/7vQUVy/zWy3fNMAACqG+w1U5w8rL6aLlpnX6OwXVEfMQkav+7Vc7VZWxud0jf2YSzKyStjyDbOtvda8WfksnusqLMi28ht9D0us9eK6X3TptMPw9ol/urZuGN/9IGjl0uRZ0tRulwERV1iG+dx9Er6E9PJZQcWmdl7ARx39/sB/LSZ3Q2gA+A8gLfuzvCEEEKMy45u6O7+IIAHB9vvpve3VLwQQojJMLnQf8OWK1Z1s2AaC6Azy6aU9H4vFkif48XTYtNM1kec7ZPbIrfrNnLbTJVMM/W11LAa3BEzUw0vkIaP26H+K2RmiTnA2XxUX07nah0IboabNKYWf5D8vEWF7mOQDC+SxqCorEoR5ygP6QM43UFtkxemQzAWp1KYJfNOjO9ZolqcZTnuKfCLa69G98bsc6yTqSLWFOWMjWz6qMYfZHGofmaeqRYvEHtBrvS42JmZVpR7/IZGof9CCDElTLSm6CXly4udQ0qM1HCnWRwU1Fwkt8X5tK+1EML7yb2xF1wkOfSfx+HBw41VOQcTVZt5f505qu1JqrS2HhZPSQRWSF0PLYoWzGRYkfc74fwGIw8BALT2U9KtggAuAKi0eWaQL7RxPvg61S/10AcHMXFwUswhz4nBOKCrEqs3cVIwVs3hO62wKqcF2DgzqKxtYiRxRlJU2agVjudFzFhhiIOOKL96dCX0TbqeSpIlxkAKXQghpoQ9UVOUosdRCUKsTTb03H4b7MucI4tNlK283cbBdDIORgKA9SPVkfsy5Q6gtW90utssERZC2gGOkQk2ebb/V1vkIhhUbpWDomhfTEHL9nC2a/diel+e8FDftbUwg+BalEGxsirnCksxoVmDApz4XNF9sLM/uSNyENNQBSRS5ZwSoXZhLWuXuT62CoKRgNyVkPdtFCh3IE+SFWqZZl9C7KMoYdhm3q7Hx8k2LsZACl0IIaaEiSr0S6H3rNAzezqQ27JZyQfRw8rRM0+J4sCi+moM8Enb7AHDtvC4j2cD9ZVgAyV7LqcFiCq3VmDXj0FWLQr9nzudFO/modzDoko2+iwFb6yNSsqe0xTHz5HlPasUBydxw8y7BrmXS5aamBQ5ANQukk2ZjumF+qoc8FM/m/ICRe8VLnjBKXOHCmGwKuftudm83Rp513Da2himz4mwYtAR2cN7bEOXbVxcIVLoQggxJUxOoVeA7sArhG3PQ+lu6ZHD+9i2Hsn8q4PZmFVva1/w+63wOIp949kGzu02DodCGKujVXljKdiND/N0IG2252MIfupj+Wjy0Z47nSu7mO4gHZ+/z+OoXSR7dbDJczKx5vlQcZ5TLmwWp7EtIiYC44Rh7KFTCSqXE4H1mqNLBwKAVQoSVwWbt60X2Mqj7Zr9yzmZVlDovtmi7bxvJ79xJckSu4kUuhBCTAkTtaFf8r/O7NXBZNlYSdu5N0yuItn3fO6ZpJZaQeVy5Cm3A4CVZ6d9s+dGn3fUudP7+WteD6jSzIBL7gG5N0tjhRNcFfvJd8nnPSry5oWkItmGXolFHfhjlKT3bVxMijIW2WZvluYzZA+O/uBc4JoiJ4fK8XFUKpexC7mO+bNYQfpYIMwU2JGlIDXttnBpOVbaYQbBqry3vg4hrgVS6EIIMSXohi6EEFPCBEP/0+Iiu+1VunGhMm2zaYaDcQCgSmt1Lcqv3sm94lBLVoGhPrgiEp9rKOQ+uv8NiIuJKKhzunprbnJpLpKJaIFXgfPunMLiq5R0jOufAnnVIzbNDOWa5xzojeIqT5354kXHrNYl9zG0sJpMFfWldKE5KAiIFaVowbV4BFkCrkqoqJSZYPhUMdSf2/GCeMxlzgurlDyrt5KXVJQLopgEUuhCCDElTDawaPA46XDq27iwGBZJt9pFAVSQTIsVedzXDsm06iujXRp7h2K7tN2iGcD8t2KaAV4UTe83lkPagubozx9nBpv7U7vmEp3ncAgsism6Lr0fKip1qPYqL8xGJZ899sOuTG1zMq2QdKu+Qh+MA5DWQvUdCqevcFBQcIOsrJLCpqRbtpp/4U6JsWyzpIJPkZtqJ8x+aIGTg4IUmi/2AlLoQggxJYyt0M2sCuA4gKfd/a6wrwngPvSLQ58D8IPu/sTYoyBxtHkgnJeET5VEWSwqz0q+xrUPwiOLC2P0QtARi8BM2QeVz4nA2Ha/cTioSFK6PKPohiLwjSUKhSe3ym6YQTRoBsG23GqJ8GS3wPb+UAhjY7SqjO6SIPUeA4u4jiinu/VYUzSuL1xqF4s/ZDspiVcnjJVqcVbWWiPfB4Iqz+p3hike28pJlXcXl7Jmso2LvcxOFPrPoLj4808AuODuLwDw6wDed6UDE0IIsTPGUuhmdgeA7wPwKwB+dkSTNwJ4z2D74wDuNTNzLyqRDsBSIqosPD8cYZzxlFR5eyFvx+p981Dajjb5Kqntesi0ysFJfFz0+uDXHO8SU/Vyoq0aXYrZs/n0glMGsNpmmzkAzFFaW04LUAmzlZmzSUVy0NGQ9wolHeOZULShVyjRVjd4pWSQV4qHLzL7JL1i5c22ck4DUAkpArJj6uxdE5NksXvRCopwSlXbW0ntFJovrifGVei/AeDnUVz/5nYATwGAu3cALAJ4VmxkZveY2XEzO97ZWI27hRBCXAHbKnQzuwvAGXf/vJm9qqjZiPeG1Lm7HwNwDADmbj7qPjh7pg6DEGsdHH3CxmL+mm3jrMI9fEJWwBZFH9m2e2QnrwUlz77xPLuIhSt4RrFBnjKthXxQjeXRxRWil8uQbbugXYts5ZxyoBPHR+q1UpKMrLaePkh7Ph8728A5LUB9MdjaaYaSla27GMLi2Q1/szikPytWwao8+o23Ryt7X8vP22U/cnmsiOuUcRT6KwHcbWZPAPgogFeb2R+GNicAHAUAM6sBOADg/C6OUwghxDZsq9Dd/V0A3gUAA4X+Tnf/4dDsfgBvAfD/ALwJwF+V2s8D7fm0HdPn1lg4kRruzOftMr/0kjOzh0p7Lt+XeYuQmG3tC+di0cfl3krCGWcupIbRNs6eLTyGoXJ85K/PHi9DycNIlbONvxI8TWoUYZp5pQS1ninv1XxQXLiCVXRnX+4bz3b4rDB0uGhcao5VfSxBl5V14zS2wW+cFXtvKdnGe5vBdUmIKeCyA4vM7L0Ajrv7/QA+COAPzOxx9JX5m3dpfEIIIcZkRzd0d38QwIOD7XfT+xsAfmA3ByaEEGJnTDT0/1ICLF6A7AQzSOMCtadZfC2spfFxVdoXzRGt/Wk7pg9g9z9eMK3EmBbax0m8omlm7sxoM0t9LdiEOC8UpwAPY+dFVjbTDIXqO5twCqrZA9ikSkmc7CsusnLOcwv50OvLdBE5mVYw72SVjVrFlY3YPdE2OVta+Iy82LlWnG+8ez79gOSCKKYdhf4LIcSUMDmFbiQkeWExeozxI4cEZgyf5wVTdlWM7o3ZwmpI/MWpdlm9R6U8e44W/6gqUT241vO+vCpTrL4zOoyfjweAxupod7pYQYnHywFIzQu5QnVyQeTryUnK+gPhFAb5GHoc/EMLl1wpCQDqF0fX7LSQxCurzckzilZxulunJFnd1ehjKhdEceMghS6EEFPCRAtcbNnOOXw+eJO1yS5daY1+HwCaFzCStVtDO/KOjzb0LEiI1Hoc0/oRch+k3E2bIQiKbej8Idl1EsjdEZuL6ZhasGW3F8gOzUnLQnIutofXaLs7mz+/OUUAuzR2g528Rq6KXAgDABrnKZ3sTFrkqIZQfXZBzApXhFqcWSzyBn/hoWbncnJB7K4Uh/QLcSMhhS6EEFPCZL1cBmKPvVe8JDgnI2SaarEXyXJ6PwbnbFKGmXqeGTWzvXO4f1TybHtnVR9t6JwKgD1HLHzIbAZQkPgrvq6yF05JsQ9OEGbRU4SUMifxqoQCGT0KHqqt5CfrUhg/91c7H0qyURGK6somvR9K0LHHCtvGz1/M+1MaWyGGkEIXQogpYaIK/RKZQq8GFUmFkbsL7CsdPEXYN5yTbIX6xjUyt3aDlwvbpdmGHtMMsJ/7xuG0zTMDIKb75QLPebuZC+nEG4fSM7YWClBwsQr2ZGkt5NeCZwCc+IuPj31kKQKCJ0tnLqnoWDC6vpzs3FxoghX5EOQPbotxgYJU+cWkyuVDLsT2SKELIcSUMFk/9K1IUQ6VDO3q5B1RI0+M0Ky6QaXWqPhFtKGzKu+FT5+VrqMT1IMTBftsN0iVV4OrNCf/ivuyduy9UuKTz4Wm2Xfd53LlXV/lYteVke8Duf96hzxbopIvKlXXb5z66M0m9R79yysrdHHZlr+eK/QORXbKh1yInSGFLoQQU4Ju6EIIMSVMdlF0MPPuLtD0vJtP92tL1aH2o2gdovB0WjAdSuJFi51DppSZ0dux6hEHOGUujDGZFlsMON14CCxqHaDxkrskm0uAECTUHB1kBMQkZhxkFBZZ6VLXVygpVnBv5EXS4ZqdaZNzllc2QrQT5S93qg7UXQy5GYQQl40UuhBCTAmTVegDhVhZTyq8V49V4Ee7KnZnQkIqUpsc7MPpcvt9pO2h0P8CFW0lHnNZCH7wwOOUvhXax7OEOCZW1DGknxU6J/iKCp0XTGtrPMB89tNrkovkxeJAnTz1bVDvrYJ0t8t5YFHvYlLivbWQQEsIsStIoQshxJSwrUI3sxkAnwHQHLT/uLv/cmjzVgC/CuDpwVv3uvvvlvXrltwGvUL271b+jOkuJIXZ3UfqsJHL0srFZDjuzrLNN6aqTa/XbsvVZvMCF40oHju7O9ZIeUfbONvo2SYf1wIaywWqPJRl5aCjKinlxlI+hWjtz8PpLxHT7DaWkrrmdLeNc2GqwQo9jMnYHbFF9TvP5dnSVMNTiKvPOCaXTQCvdvcVM6sD+Bsz+7S7fza0+5i7v333hyiEEGIctr2hu7sDuKQ164O/En+TMeHAojkKBY9GIIq0seU03MpSrkKdgpN8hirMn88/InusDAUdkYpmu3m0UbO9ntV6bJcVzKBz1UKQUWsfh+pToE41vxic4KvaZo+XUBZuhWYo9E1V4loA7autpQEOlYUj27htxGxnaUrRe+Zs2m6FBQAhxFVnLBu6mVXN7GEAZwA84O6fG9Hs+83sUTP7uJkdLejnHjM7bmbHu6uro5oIIYS4TMbycnH3LoAXm9lBAJ80s+909y9Skz8D8BF33zSztwH4EIBXj+jnGIBjANA8etQ7c32JWF1Jw+jNhjJpPMJmUp69IOXrS+l1hyR0bzb6TXOa3ZCsikQ/K/RY7o6Tc7GfezsUuObi0tUxbfLsvcKFKoCQMbhkjsSh+j1K4cuFKiIVKtxcWc2nELZC9u/godI9lyqGKIGWEJNlR14u7kWdWskAABc2SURBVH4RwIMAXhfeP+ful+4CvwPgu3ZldEIIIcZm2xu6md00UOYws1kArwHwldDmNnp5N4DHdnOQQgghtmcck8ttAD5kZlX0HwB/5O5/bmbvBXDc3e8H8NNmdjeADoDzAN66XafmyYWwTaH/HGQEIAst5yCjaEppsfWEbBOVdv7MMrI6RFMK1xvljI21EAfD5pjNQ2l7/mRuB9mkkP4uhePHLI916p/7rgWTS69mI9tFOnPpMzeWixty3U/OZc4h/AAyMwubWACZWYTYS4zj5fIogJeMeP/dtP0uAO/a3aEJIYTYCRML/XcDeoNc516jRbxQYYjdEa1NC5rBHbFzE1WwX04qn4OWAMA413eJuMzSBxzI93Ht0C4FE60fKa6ixNRCjE3uIknulzEvebvAbbEe86GT2yL1V13PB2Qdnv7QdbqYF1vVwqcQ1wcK/RdCiClhshWLBmrcOIy/livq5nyy7W4uJqN3+1m5Dd0oZUCvSaq+FRJSUQWk+DRr76cAHxKpnVB7lNPncoqAaJNn2zunBYgujBySz6H/ng89t5uTomZFDgQlbqET7o/t5meTCu9QgJAQ4vpBCl0IIaaEiabPvWTfbswlydpazjNctdbSa2P1HgNrGhS6vpkM4JVe/sziwCXv5R41nP62vW/0+0Bel7RO9RnKPE8a1G795riPgok2k9qOofrsHcMpd6utEDxlHJyUrkv17HLebjllD+uczb1XhBDXH1LoQggxJUxOoXtKbdve4IxZwSulSl4f3fT8qTRy+dpbTu4xjXOpv6iuN26nY5q5st08TAUfVpPKjbZsDv1nr5x68FfnfWyHnz0T+mOPFQ7VD37ozSXyXiG7uwc7eXWT/MvXyWC/lCt02cqFmC6k0IUQYkqYrJfLQH1X60ltdz23a1uFFXp6v9HMfao3SKG3bqZ0r5shFewGFYlYy/dlZex4V1DonDDMSABH2/jCkwXHhCLMWRpbSpEbbfIVspVz4ebqWu42U12kKcT5i1ubUuRCTDdS6EIIMSXohi6EEFPCnjC51Cj0f3amOHH4ylIqAbSxOJPv5DzqnJAr5kNfpNzr9dz0UVtLtpXWQTJpbBY/97IaoGFRtEWujxzuH90ROc851wf1asjXTqkAKpvF+cs5dF9mFiFuHKTQhRBiSpis2+KgepBx3VDLVfPKSlLi9dm0ENoNKQI6q1RvdIZU7sU8UClLwTsf0gdQoFFthZ51YVG0Q5WJOB3vUAASnTqL2g9Jt2Yujg6Yigq9fi4p8epSWvi0M3lQkFS5EDcmUuhCCDElTE6hVxy9ub5u7XSSMo4KvUYBRK2VJHkb83lV+W6FbOObxf1lhTF6uQK+FOjU36ahBrM+7+MQ/Fi4okkJvrg+aH0tH1O3mc7LxSnmTuaSv7JJ7phn5Y4ohMiRQhdCiClhW4VuZjMAPgOgOWj/cXf/5dCmCeA+9ItDnwPwg+7+RGnHblspbyvVpJpr1WDXJhFdO5DsxusrMVctBd1wMFIrBBatknoPCr07k47j0P/WgTBroH1sJ4/eKxyc5CWPzplzlECLg4c286lB5fSFre3OGalyIUTOOAp9E8Cr3f1FAF4M4HVm9orQ5icAXHD3FwD4dQDv291hCiGE2I5xaoo6gEt5VuuDv5i89o0A3jPY/jiAe83MBseOptpD5UDfDt6sJ4U6U89D+quVpFjPryT3kko9FHUghd5aJdkcEnBVVsjWPhMSgXWS8m7vT8fNPJM/96I3y9apLoSPSxMAo2FYt/iy1BbT2kDlwmq2r3uasnp58K8XQtzwjGVDN7OqmT0M4AyAB9z9c6HJ7QCeAgB37wBYBPCsEf3cY2bHzex4d3k17hZCCHEFjOXl4u5dAC82s4MAPmlm3+nuX6Qmo+qcDclQdz8G4BgANL/tDu+1+/bsLhWhuGluJTvmyYuHtrZvOZDSv6638mrS55fmt7YrlOyrt5m36x5I+xqn84+flYnboPS5eb4wdClItUbPpV4tvwyz5yndLalyLmIBAPXlZCuvPpO8V7onTmbtVKBZCFHGjrxc3P0igAcBvC7sOgHgKACYWQ3AAQAqgSOEENeQbW/oZnbTQJnDzGYBvAbAV0Kz+wG8ZbD9JgB/VWo/F0IIseuMY3K5DcCHzKyK/gPgj9z9z83svQCOu/v9AD4I4A/M7HH0lfmbt+21B2CQm3zfLSk65/GzR7JmR/Ylm8bZ5WRW6XRzO0iPEnJ5h55TwW2xvpxed2eL85JXN8k1MXpIUrpxTs7VXC5eqOQc6I2LeVBU7Xz6jL1Tp9NwZGIRQuyAcbxcHgXwkhHvv5u2NwD8wO4OTQghxE6YXOh/1VFZ6C8Gsmvi4YU8B+3yRpLHM43k0njhQpDNDKlyO5ir4dYMBRZt5CqfE3JtHk6KevZbIUUACfG8ylHejtPzNi8mtV1dzv0ee08+nbZb+XiFEGJcFPovhBBTwuQUOgw+kLdtsodHt8WDM8lgfXFjdmt7biEv6tDtkiqnxF0ba3n63Mp68EHkfWyypiCjXt5FZkPPAoZCrdAG2dTrF2m8J05n7Xrr6xBCiCtFCl0IIaaECZagc1QG4fqNWrKNt3rBrk329Q6p8PXVXDY3ZlIfG7Svt5F/RJun0m2L+T5Ok8upcKuhwhu/nn2GEott5F4uDSpIUXkqqfLOxYsQQojdRgpdCCGmhIkpdDOgPgjRr5PxeqFe7OVxeC55wDRDEq9TZw9ubTsp+cpM3i5T7MENvUdZAiqkwmMyrsYSpdklVV5fytPd1r5FRSjOPAMhhLiaSKELIcSUMEEvl8Tz9qW0L2fW92X7Zsi+fmF9DkVYhd1NKGq0lz+zrJ28VzgBF5AXoZil+hEWAjYrlGjLOhQBenIxa9d96kTheIUQYreRQhdCiClBN3QhhJgSJrcoCkdlYCZ5ciXlPGc3RQBYaqUQfw46emrxYNZuYT6tYi5Rci6r5iufPXJV7DXzffXFUWnd8+AhIDfBzJyixFpPKX+5EGJySKELIcSUsCcWRW+eTcq75/kz5vTawtb2IoX+b4RKRHOzyd2xS66J1g2JtWi7upbvy9LiblBiraVcyc+eTn6M9tSpdN7NgmKjQghxDZBCF0KIKWGiCv1Scq4jjWSHPt/OXRNvIbv506v7t7YX5vJ4/IsXR7s02mweWGSrySY/FNJPMU1sN6+t5bbw+pPJp7Fz/sLI8wohxLVGCl0IIaaEbRW6mR0FcB+AW9EvHHfM3d8f2rwKwJ8C+ObgrU+4+3vL+3XMNvuh8vM1SmIVXErOtpINfbOdhnv7vjyIZ2U9Ke/eUkrOVTud29q786l/r4aCFPR4ay6ldrNP5OfqBG8WIYTYC4xjcukA+Dl3f8jM9gH4vJk94O5fDu3+2t3v2v0hCiGEGIdxaoqeAnBqsL1sZo8BuB1AvKHviIoBjWrfNn1bIyWxevDct2ftFki9b3bScGMx6Xot2blbJLy5DBwAVDbIRz34l89c4ELOyfbuTzyVN/TiYtBCCDEpdmRDN7PnoV8w+nMjdn+3mT1iZp82s+8oOP4eMztuZsfbi2ujmgghhLhMxr6hm9kCgD8G8A53Xwq7HwLwXHd/EYDfBPAno/pw92Pufqe731k/UJxoSwghxM4Zy23RzOro38w/7O6fiPv5Bu/unzKz3zKzI+5+NrZNfTpmR+Q+f978uez1tzaSq+J8M7Vf8WbWbmVpJr2okJklPLLIuoPZM/m+udPJzNL8alr47KjmpxDiOmBbhW5mBuCDAB5z918raHProB3M7GWDfs+NaiuEEOLqMI5CfyWAHwHwBTN7ePDeLwJ4DgC4+wcAvAnAT5pZB8A6gDe7u4/q7BIV861c51WkRcae566EM9VUBegFB5Lg//zaHVk7I1XuTaobejb/iFx9qLaeD3HmdApw6pw6DSGEuJ4Yx8vlb5CnQBnV5l4A9+7WoIQQQuyciYX+uwOdQSTP+U4KHrq1uVh0CB45f/vWdrudD73XrqYXq7QdHkWs0GfO52kB8JVvpm25JgohrjMU+i+EEFPCxBR6vdLDrbPLAIDDtZXCdp9dff7WdruXlHejkavrdkHKXPZqAYDmYrKbz335W9k+ebMIIa5npNCFEGJKmJhC73gF5zf7wUVtT8OYq+Q5bdnrpUrl6TY386F7Kz2bZs6m7VoQ3QtPJSN658mnL2PkQgixN5FCF0KIKWFiCr1qPSwMIkXr1sneZ/Y3kqK+sJlK0FUqIenWehWjmDuT99d4LCXa6siTRQgxRUihCyHElKAbuhBCTAkTM7k0Kl0cnevX4zxQTal0F7t5FsalVkq6VSNzTGs5T85VW0vPJo5N2vdEnqa380xhvjAhhLiukUIXQogpYWIKvesVXGz3FznPdlKK3FOtA1m7WwbBRwDw1NLBrW2rd7N21k21Q2fOpQXType+kZ/3CsYshBB7GSl0IYSYEiYX+m8d3DHTt6F/z9zXtt7/a/zjrN3/OfN8jMLO5TZ0DvE/+JWk6rsrxWkFhBBimpBCF0KIKWFy6XNh2Oj17d4bnuzfi53cy+U5+y5sbT/yrWenHZYHFs1zPYovf333BiqEENcJUuhCCDElbKvQzewogPsA3AqgB+CYu78/tDEA7wfwBgBrAN7q7g+V9gvHTKVfXu4w+aHfFvLd/t+z37a1vbHe2NqeOZ0/iw59MTmf95QGVwhxAzKOyaUD4Ofc/SEz2wfg82b2gLt/mdq8HsALB38vB/Dbg3+FEEJcI8apKXoKwKnB9rKZPQbgdgB8Q38jgPsGhaE/a2YHzey2wbEjqZpjYVAP7iT5oa/1GkWHoLeY9tWXw87HvgEhhLiR2ZEN3cyeB+AlAD4Xdt0O4Cl6fWLwXjz+HjM7bmbHV863djZSIYQQpYx9QzezBQB/DOAd7r4Ud484xIfecD/m7ne6+50Lh4uVuBBCiJ0zltuimdXRv5l/2N0/MaLJCQBH6fUdAE6W9dmwDo7WzwMA5iyp9Vtqi1m7jU5yaayupufPzQ/lSbe0ECqEuNHZVqEPPFg+COAxd/+1gmb3A/hR6/MKAItl9nMhhBC7zzgK/ZUAfgTAF8zs4cF7vwjgOQDg7h8A8Cn0XRYfR99t8ce267TrFSz3+qlxG5ZSZp3u5Mm5LqymKkUzZ6m+6Je+mfc3xgcRQohpZhwvl7/BaBs5t3EAP7VbgxJCCLFzJpqc69n1flj/wcrm1vvfWL8pa9dqpSHe+tWkw7uLua1dCCFudBT6L4QQU8LEFHrbazjZPgQAeF4tJeBqVjpZu875VIJu36MpA1feSgghhBS6EEJMCRNT6M1KG89vnAEAbHh16/2T6wezdvseT/s633zy2gxOCCGuQ6TQhRBiStANXQghpoSJmVwqcOwbuCteCjACgDPr81m7I49SEi/vXZOxCSHE9YgUuhBCTAkTVeiNQcD+BpLyfuLUkazdP3nkia1tuSoKIUQxUuhCCDElTEyhNwx47uDsf7F2aOv95t/PZO06z5y9lsMSQojrFil0IYSYEiam0LvuuNjrW8W/uH7H1vs3/50s5UIIcTlIoQshxJQwMYUOA6qDLOvn2sn3fP7hp7Nm0utCCDEeUuhCCDElbKvQzez3ANwF4Iy7f+eI/a8C8KcALtWE+4S7v3e7fqsw7LP+6R988gVb799+4kvjjFsIIURgHJPL7wO4F8B9JW3+2t3v2pURCSGEuCy2Nbm4+2cAnL8GYxFCCHEF7Nai6Heb2SMATgJ4p7uPtJuY2T0A7gGA59xew/7KbH/H3x7YpWEIIcSNy24sij4E4Lnu/iIAvwngT4oauvsxd7/T3e+86VnVomZCCCEugytW6O6+RNufMrPfMrMj7j52zP6zvty90mEIIcQNzxUrdDO71cxssP2yQZ/nrrRfIYQQO2Mct8WPAHgVgCNmdgLALwOoA4C7fwDAmwD8pJl1AKwDeLO7+3b9dtHDUm8dALDwdymYSIFEQghxeWx7Q3f3H9pm/73ouzUKIYSYIBML/V/s1fGptWcDADpPnZzUMIQQYmpQ6L8QQkwJE1Po59oL+PDJl/df+KlJDUMIIaYGKXQhhJgSJqbQNzbr+NI3bwcAvBBS6EIIcaVIoQshxJSgG7oQQkwJEzO5VDYNs483JnV6IYSYOqTQhRBiSpigQgf2P7FthgAhhBBjIoUuhBBTwsQUenWjhwNfXQEASKcLIcSVI4UuhBBTwsQUum22UPnGCQCAylsIIcSVI4UuhBBTgo1Ri+LqnNjsGQD/MJGT75wjAMYuqbcH0HivHtfTWAGN92oyqbE+191vGrVjYjf06wkzO+7ud056HOOi8V49rqexAhrv1WQvjlUmFyGEmBJ0QxdCiClBN/TxODbpAewQjffqcT2NFdB4ryZ7bqyyoQshxJQghS6EEFOCbuhCCDEl6IY+wMyOmtn/NrPHzOxLZvYzI9q8yswWzezhwd+7JzFWGs8TZvaFwViOj9hvZvZfzOxxM3vUzF46oXF+O12zh81syczeEdpM9Nqa2e+Z2Rkz+yK9d9jMHjCzrw3+PVRw7FsGbb5mZm+Z4Hh/1cy+MviuP2lmBwuOLf3dXMPxvsfMnqbv/A0Fx77OzL46+B3/woTG+jEa5xNm9nDBsdf82ma4u/766wi3AXjpYHsfgL8H8E9Dm1cB+PNJj5XG8wSAIyX73wDg0wAMwCsAfG4PjLkK4FvoB0fsmWsL4HsBvBTAF+m9/wTgFwbbvwDgfSOOOwzgG4N/Dw22D01ovK8FUBtsv2/UeMf53VzD8b4HwDvH+L18HcDzATQAPBL/X16LsYb9/xnAu/fKteU/KfQB7n7K3R8abC8DeAzA7ZMd1RXzRgD3eZ/PAjhoZrdNeEz/EsDX3X1PRQm7+2cAnA9vvxHAhwbbHwLwb0Yc+q8BPODu5939AoAHALzuqg10wKjxuvtfuntn8PKzAO642uMYl4LrOw4vA/C4u3/D3VsAPor+93LVKBurmRmAfw/gI1dzDJeLbugjMLPnAXgJgM+N2P3dZvaImX3azL7jmg5sGAfwl2b2eTO7Z8T+2wE8Ra9PYPIPqTej+D/DXrq2AHCLu58C+g98ADePaLMXrzEA/Dj6s7NRbPe7uZa8fWAi+r0Ck9Zeu77/HMBpd/9awf6JXlvd0ANmtgDgjwG8w92Xwu6H0DcVvAjAbwL4k2s9vsAr3f2lAF4P4KfM7HvDfhtxzMT8VM2sAeBuAP99xO69dm3HZU9dYwAws18C0AHw4YIm2/1urhW/DeAfAXgxgFPomzIie+36/hDK1flEr61u6ISZ1dG/mX/Y3T8R97v7kruvDLY/BaBuZkeu8TB5PCcH/54B8En0p6fMCQBH6fUdAE5em9GN5PUAHnL303HHXru2A05fMlEN/j0zos2eusaDRdm7APwHHxh1I2P8bq4J7n7a3bvu3gPwOwXj2DPX18xqAP4dgI8VtZn0tdUNfcDANvZBAI+5+68VtLl10A5m9jL0r9+5azfKbCzzZrbv0jb6C2JfDM3uB/CjA2+XVwBYvGRCmBCF6mYvXVvifgCXvFbeAuBPR7T5CwCvNbNDA5PBawfvXXPM7HUA/iOAu919raDNOL+ba0JYz/m3BeP4WwAvNLNvG8zw3oz+9zIJXgPgK+5+YtTOPXFtJ7Uau9f+APwz9KdyjwJ4ePD3BgBvA/C2QZu3A/gS+ivtnwXwPRMc7/MH43hkMKZfGrzP4zUA/xV9L4EvALhzguOdQ/8GfYDe2zPXFv0HzSkAbfRV4U8AeBaA/wXga4N/Dw/a3gngd+nYHwfw+ODvxyY43sfRtzdf+v1+YND22QA+Vfa7mdB4/2Dwu3wU/Zv0bXG8g9dvQN/r7OvXYryjxjp4//cv/V6p7cSvLf8p9F8IIaYEmVyEEGJK0A1dCCGmBN3QhRBiStANXQghpgTd0IUQYkrQDV0IIaYE3dCFEGJK+P8gbEKctzgPdwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist2d(q2_array/1e6,np.sqrt(Bbar_array)/1e3, weights=weights_np, density=True,bins=100);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3AAAAEvCAYAAAAErSPcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAYtklEQVR4nO3df6ydd30f8PcH08AErKONtVE7Xkxrpoa2gnFx/kCjE+WHmau4f4AwiCmoVBZTojIhtIGKQDNDMlRi8EfWYYVItB3z+LEhq02XUQHb0BrwDT/KnCzCGEMuRiIlGQzBkjp89sc9tKc3ju9zc6997nPO6yVZPs/zfJ/jz5FOHL/v5/v9PtXdAQAAYPt7wqwLAAAAYBgBDgAAYCQEOAAAgJEQ4AAAAEZCgAMAABgJAQ4AAGAknjjrAta6+uqr+9prr511GQAAADNx1113/UV377zYtW0X4K699tosLy/PugwAAICZqKpvPNY1UygBAABGQoADAAAYCQEOAABgJAQ4AACAkRDgAAAARkKAAwAAGAkBDgAAYCQEOAAAgJEQ4AAAAEZCgAMAABgJAQ4AAGAknjhkUFUdSPL+JDuS3Nrdxx5j3CuSfDTJ87t7eXLurUlen+SRJL/d3XdsReEAAJtx7Vv++DGvnTt28ApWAjDcugGuqnYkuSXJS5KsJDlVVSe7++41456W5LeTfG7q3HVJDid5dpKfS/KnVfWs7n5k6z4CADBLlwpCiTAEsJWGdOD2JznT3WeTpKpOJDmU5O41496Z5D1J3jx17lCSE939UJKvV9WZyfv92WYLBwDGbzPh73IGx/Xeez1CK3C5DAlwu5LcN3W8kuT66QFV9dwk13T3H1XVm9fce+eae3c9zloBAAbbbAgD2I6GBLi6yLn+q4tVT0jyb5K8bqP3Tr3HkSRHkmTPnj0DSgIAxmIRg5T1dcDlMiTArSS5Zup4d5LzU8dPS/JLST5TVUny95KcrKobBtybJOnu40mOJ8nS0tKjAh4AwKIQ/oBLGRLgTiXZV1V7k3wrq5uSvOYnF7v7e0mu/slxVX0myZu7e7mqfpTkw1X13qxuYrIvyee3rnwAYJ7NY/duHj8TcOWsG+C6+0JV3Zzkjqw+RuC27j5dVUeTLHf3yUvce7qqPpLVDU8uJLnJDpQAAACPz6DnwHX37UluX3Pu7Y8x9h+vOX5Xknc9zvoAgC2wma6PaXsA28cTZl0AAAAAwwhwAAAAIzFoCiUAsLhsurF9XM6HlwPjoAMHAAAwEjpwAAALwjPmYPx04AAAAEZCBw4A5oS1avgOwPzTgQMAABgJAQ4AAGAkBDgAAICRsAYOAADPmIOREOAAYJvwD2gA1iPAAcBI2GGQWfIMOdgeBDgAuIKEMAA2wyYmAAAAI6EDBwBbSIeNRWT9Jlw5OnAAAAAjIcABAACMhAAHAAAwEtbAAQBwWVkjB1tHBw4AAGAkBDgAAICREOAAAABGYtAauKo6kOT9SXYkubW7j625/oYkNyV5JMkPkhzp7rur6tok9yS5dzL0zu5+w9aUDgCz4VlvsLUu9d+U9XHwN60b4KpqR5JbkrwkyUqSU1V1srvvnhr24e7+d5PxNyR5b5IDk2tf6+7nbG3ZAAAAi2fIFMr9Sc5099nufjjJiSSHpgd09/enDp+SpLeuRAAAAJJhUyh3Jblv6nglyfVrB1XVTUnelOSqJC+aurS3qr6Y5PtJ3tbd/+PxlwsAl58pkgBsV0MCXF3k3KM6bN19S5Jbquo1Sd6W5MYk306yp7u/W1XPS/KJqnr2mo5dqupIkiNJsmfPng1+BADYGAENgLEaEuBWklwzdbw7yflLjD+R5PeSpLsfSvLQ5PVdVfW1JM9Ksjx9Q3cfT3I8SZaWlky/BGDThDQA5tGQAHcqyb6q2pvkW0kOJ3nN9ICq2tfdX50cHkzy1cn5nUke6O5HquqZSfYlObtVxQMAMN/W+2GMXSpZNOsGuO6+UFU3J7kjq48RuK27T1fV0STL3X0yyc1V9eIkf5nkwaxOn0ySFyY5WlUXsvqIgTd09wOX44MAAADMu0HPgevu25Pcvubc26dev/Ex7vt4ko9vpkAAAABWDXmMAAAAANuAAAcAADASAhwAAMBICHAAAAAjIcABAACMxKBdKAFgu/GgbgAWkQAHAMBoXeqHOR7yzTwS4ADYtnTZAOBvsgYOAABgJHTgAACYS+t18U2xZIwEOABmxhRJANgYUygBAABGQoADAAAYCQEOAABgJAQ4AACAkRDgAAAARsIulABcNnaZBICtpQMHAAAwEjpwAAAsJA/6Zox04AAAAEZCgAMAABgJAQ4AAGAkBDgAAICRGBTgqupAVd1bVWeq6i0Xuf6GqvpKVX2pqj5bVddNXXvr5L57q+plW1k8AADAIlk3wFXVjiS3JHl5kuuSvHo6oE18uLt/ubufk+Q9Sd47ufe6JIeTPDvJgST/dvJ+AAAAbNCQDtz+JGe6+2x3P5zkRJJD0wO6+/tTh09J0pPXh5Kc6O6HuvvrSc5M3g8AAIANGvIcuF1J7ps6Xkly/dpBVXVTkjcluSrJi6buvXPNvbseV6UAXBbrPQdpPZ6TBMyrS/396O8+ZmVIB64ucq4fdaL7lu7++ST/MsnbNnJvVR2pquWqWr7//vsHlAQAALB4hnTgVpJcM3W8O8n5S4w/keT3NnJvdx9PcjxJlpaWHhXwANi+NtvBAwCGG9KBO5VkX1XtraqrsropycnpAVW1b+rwYJKvTl6fTHK4qp5UVXuT7Evy+c2XDQAAsHjW7cB194WqujnJHUl2JLmtu09X1dEky919MsnNVfXiJH+Z5MEkN07uPV1VH0lyd5ILSW7q7kcu02cBAACYa0OmUKa7b09y+5pzb596/cZL3PuuJO96vAUCAACwatCDvAEAAJg9AQ4AAGAkBDgAAICRGLQGDgAA+GvrPULFg765XAQ4gAXgWW0AMB9MoQQAABgJAQ4AAGAkTKEEmAOmSALAYhDgAABgi13qB2s2OGEzTKEEAAAYCQEOAABgJAQ4AACAkRDgAAAARsImJgAjYadJAEAHDgAAYCQEOAAAgJEQ4AAAAEbCGjiAbcIaNwBgPTpwAAAAIyHAAQAAjIQABwAAMBICHAAAwEgIcAAAACMxaBfKqjqQ5P1JdiS5tbuPrbn+piS/leRCkvuT/GZ3f2Ny7ZEkX5kM/WZ337BFtQMAwOist+vwuWMHr1AljNG6Aa6qdiS5JclLkqwkOVVVJ7v77qlhX0yy1N0/rKp/luQ9SV41ufaj7n7OFtcNAACwcIZModyf5Ex3n+3uh5OcSHJoekB3f7q7fzg5vDPJ7q0tEwAAgCFTKHcluW/qeCXJ9ZcY//okfzJ1/OSqWs7q9Mpj3f2JDVcJMAc8qBsA2KwhAa4ucq4vOrDqtUmWkvzq1Ok93X2+qp6Z5FNV9ZXu/tqa+44kOZIke/bsGVQ4AADAohkyhXIlyTVTx7uTnF87qKpenOR3ktzQ3Q/95Hx3n5/8fjbJZ5I8d+293X28u5e6e2nnzp0b+gAAAACLYkgH7lSSfVW1N8m3khxO8prpAVX13CQfSHKgu78zdf7pSX7Y3Q9V1dVJXpDVDU4A5pJpkgDA5bRugOvuC1V1c5I7svoYgdu6+3RVHU2y3N0nk/xukqcm+WhVJX/9uIBfTPKBqvpxVrt9x9bsXgkAAMBA1X3R5Wwzs7S01MvLy7MuA+Bx0YED4HLznLj5V1V3dffSxa4NWQMHAADANjBkDRwAEzpsAMAs6cABAACMhAAHAAAwEgIcAADASAhwAAAAIyHAAQAAjIRdKAHWsNMkALBd6cABAACMhA4cAACMyKVmipw7dvAKVsIs6MABAACMhAAHAAAwEgIcAADASAhwAAAAI2ETE2DheEwAADBWOnAAAAAjoQMHzB0dNgBgXunAAQAAjIQABwAAMBKmUAIAwJxYbxnBuWMHr1AlXC46cAAAACMhwAEAAIyEAAcAADASg9bAVdWBJO9PsiPJrd19bM31NyX5rSQXktyf5De7+xuTazcmedtk6L/u7g9tUe3AAvOoAABgEa3bgauqHUluSfLyJNcleXVVXbdm2BeTLHX3ryT5WJL3TO79mSTvSHJ9kv1J3lFVT9+68gEAABbHkCmU+5Oc6e6z3f1wkhNJDk0P6O5Pd/cPJ4d3Jtk9ef2yJJ/s7ge6+8Ekn0xyYGtKBwAAWCxDAtyuJPdNHa9Mzj2W1yf5k43cW1VHqmq5qpbvv//+ASUBAAAsniFr4Ooi5/qiA6tem2Qpya9u5N7uPp7keJIsLS1d9L2BxWKNGwDAow3pwK0kuWbqeHeS82sHVdWLk/xOkhu6+6GN3AsAAMD6hgS4U0n2VdXeqroqyeEkJ6cHVNVzk3wgq+HtO1OX7kjy0qp6+mTzkpdOzgEAALBB606h7O4LVXVzVoPXjiS3dffpqjqaZLm7Tyb53SRPTfLRqkqSb3b3Dd39QFW9M6shMEmOdvcDl+WTAAAAzLnq3l5LzpaWlnp5eXnWZQAzZg0cAFx5544dnHUJJKmqu7p76WLXhkyhBAAAYBsQ4AAAAEZCgAMAABiJIc+BA7gsrHMDANgYHTgAAICREOAAAABGQoADAAAYCWvggMvGGjcAGJdL/b/bM+K2Bx04AACAkRDgAAAARkKAAwAAGAkBDgAAYCQEOAAAgJGwCyXwuNllEgDgytKBAwAAGAkBDgAAYCQEOAAAgJEQ4AAAAEZCgAMAABgJu1ACl2SnSQCA7UMHDgAAYCQEOAAAgJEQ4AAAAEZi0Bq4qjqQ5P1JdiS5tbuPrbn+wiTvS/IrSQ5398emrj2S5CuTw2929w1bUTiwNaxxAwCGWO/fDOeOHbxClSy2dQNcVe1IckuSlyRZSXKqqk52991Tw76Z5HVJ3nyRt/hRdz9nC2oFAABYaEM6cPuTnOnus0lSVSeSHEryVwGuu89Nrv34MtQIAABAhq2B25Xkvqnjlcm5oZ5cVctVdWdV/cbFBlTVkcmY5fvvv38Dbw0AALA4hgS4usi53sCfsae7l5K8Jsn7qurnH/Vm3ce7e6m7l3bu3LmBtwYAAFgcQ6ZQriS5Zup4d5LzQ/+A7j4/+f1sVX0myXOTfG0DNQKbZKMSAID5MKQDdyrJvqraW1VXJTmc5OSQN6+qp1fVkyavr07ygkytnQMAAGC4dQNcd19IcnOSO5Lck+Qj3X26qo5W1Q1JUlXPr6qVJK9M8oGqOj25/ReTLFfVl5N8OsmxNbtXAgAAMNCg58B19+1Jbl9z7u1Tr09ldWrl2vv+Z5Jf3mSNAAAAZNgUSgAAALaBQR04YHuzSQkAwGLQgQMAABgJAQ4AAGAkTKEEAAA27VJLOs4dO3gFK5lvOnAAAAAjIcABAACMhCmUMAJ2mQQAINGBAwAAGA0BDgAAYCRMoYRtwjRJAADWowMHAAAwEgIcAADASAhwAAAAIyHAAQAAjIQABwAAMBJ2oYQrxC6TAABslg4cAADASOjAAQAAl9V6M5HOHTt4hSoZPx04AACAkRDgAAAARsIUStgiNikBAOByG9SBq6oDVXVvVZ2pqrdc5PoLq+oLVXWhql6x5tqNVfXVya8bt6pwAACARbNugKuqHUluSfLyJNcleXVVXbdm2DeTvC7Jh9fc+zNJ3pHk+iT7k7yjqp6++bIBAAAWz5AO3P4kZ7r7bHc/nOREkkPTA7r7XHf/eZIfr7n3ZUk+2d0PdPeDST6Z5MAW1A0AALBwhgS4XUnumzpemZwbYjP3AgAAMGXIJiZ1kXM98P0H3VtVR5IcSZI9e/YMfGu48mxUAgDALA3pwK0kuWbqeHeS8wPff9C93X28u5e6e2nnzp0D3xoAAGCxDAlwp5Lsq6q9VXVVksNJTg58/zuSvLSqnj7ZvOSlk3MAAABs0LoBrrsvJLk5q8HrniQf6e7TVXW0qm5Ikqp6flWtJHllkg9U1enJvQ8keWdWQ+CpJEcn5wAAANigQQ/y7u7bk9y+5tzbp16fyur0yIvde1uS2zZRIwAAABkY4GBR2KQEAIDtbMgaOAAAALYBAQ4AAGAkTKEEAABmar1lLOeOHbxClWx/AhwLxzo3AADGyhRKAACAkRDgAAAARkKAAwAAGAkBDgAAYCRsYsLcsUkJAADzSgcOAABgJAQ4AACAkRDgAAAARsIaOEbHGjcAABaVDhwAAMBICHAAAAAjIcABAACMhAAHAAAwEjYxAQAAtrVLbWJ37tjBK1jJ7AlwbEt2mgQAgEczhRIAAGAkBDgAAICREOAAAABGYtAauKo6kOT9SXYkubW7j625/qQkv5/keUm+m+RV3X2uqq5Nck+SeydD7+zuN2xN6YyZNW4AALBx6wa4qtqR5JYkL0mykuRUVZ3s7runhr0+yYPd/QtVdTjJu5O8anLta939nC2uGwAAYOEMmUK5P8mZ7j7b3Q8nOZHk0Joxh5J8aPL6Y0l+rapq68oEAABgSIDbleS+qeOVybmLjunuC0m+l+RnJ9f2VtUXq+q/VdU/2mS9AAAAC2vIGriLddJ64JhvJ9nT3d+tqucl+URVPbu7v/83bq46kuRIkuzZs2dASYyBdW4AALC1hnTgVpJcM3W8O8n5xxpTVU9M8tNJHujuh7r7u0nS3Xcl+VqSZ639A7r7eHcvdffSzp07N/4pAAAAFsCQAHcqyb6q2ltVVyU5nOTkmjEnk9w4ef2KJJ/q7q6qnZNNUFJVz0yyL8nZrSkdAABgsaw7hbK7L1TVzUnuyOpjBG7r7tNVdTTJcnefTPLBJH9QVWeSPJDVkJckL0xytKouJHkkyRu6+4HL8UEAAADm3aDnwHX37UluX3Pu7VOv/1+SV17kvo8n+fgma2SbssYNAACurCFTKAEAANgGBDgAAICRGDSFksVkiiQAANvdev9mPXfs4BWq5MrQgQMAABgJAQ4AAGAkBDgAAICRsAZuwVnnBgAA46EDBwAAMBICHAAAwEiYQjnnTJEEAID5oQMHAAAwEgIcAADASJhCOQdMkwQAgMWgAwcAADASAhwAAMBImEI5AqZIAgAAiQ4cAADAaAhwAAAAI2EKJQAAMLcutRzp3LGDV7CSrSHAbQPWuAEAAEOYQgkAADASOnBXiC4bAACwWTpwAAAAIzEowFXVgaq6t6rOVNVbLnL9SVX1HyfXP1dV105de+vk/L1V9bKtKx0AAGCxrDuFsqp2JLklyUuSrCQ5VVUnu/vuqWGvT/Jgd/9CVR1O8u4kr6qq65IcTvLsJD+X5E+r6lnd/chWf5BZM0USAAC43IZ04PYnOdPdZ7v74SQnkhxaM+ZQkg9NXn8sya9VVU3On+juh7r760nOTN4PAACADRqyicmuJPdNHa8kuf6xxnT3har6XpKfnZy/c829ux53tTOkwwYAAMzakABXFznXA8cMuTdVdSTJkcnhD6rq3gF1DXV1kr/YwveD7cD3mnnke8088r1mHs3N97rePesKHtPff6wLQwLcSpJrpo53Jzn/GGNWquqJSX46yQMD7013H09yfEAtG1ZVy929dDneG2bF95p55HvNPPK9Zh75Xs/WkDVwp5Lsq6q9VXVVVjclOblmzMkkN05evyLJp7q7J+cPT3ap3JtkX5LPb03pAAAAi2XdDtxkTdvNSe5IsiPJbd19uqqOJlnu7pNJPpjkD6rqTFY7b4cn956uqo8kuTvJhSQ3zeMOlAAAAFdCrTbK5ldVHZlM0YS54XvNPPK9Zh75XjOPfK9na+4DHAAAwLwYsgYOAACAbWBuA1xV3VZV36mq/zXrWmCrVNU1VfXpqrqnqk5X1RtnXRNsVlU9uao+X1Vfnnyv/9Wsa4KtUlU7quqLVfVHs64FtkJVnauqr1TVl6pqedb1LKK5nUJZVS9M8oMkv9/dvzTremArVNUzkjyju79QVU9LcleS3+juu2dcGjxuVVVJntLdP6iqn0ry2SRv7O47Z1wabFpVvSnJUpK/3d2/Put6YLOq6lySpe6ei+fAjdHcduC6+79ndUdMmBvd/e3u/sLk9f9Nck+SXbOtCjanV/1gcvhTk1/z+dNFFkpV7U5yMMmts64FmB9zG+Bg3lXVtUmem+Rzs60ENm8yzexLSb6T5JPd7XvNPHhfkn+R5MezLgS2UCf5r1V1V1UdmXUxi0iAgxGqqqcm+XiSf97d3591PbBZ3f1Idz8nye4k+6vK1HdGrap+Pcl3uvuuWdcCW+wF3f0Pk7w8yU2TZUtcQQIcjMxkjdDHk/z77v5Ps64HtlJ3/58kn0lyYMalwGa9IMkNk/VCJ5K8qKr+cLYlweZ19/nJ799J8p+T7J9tRYtHgIMRmWz28MEk93T3e2ddD2yFqtpZVX9n8vpvJXlxkv8926pgc7r7rd29u7uvTXI4yae6+7UzLgs2paqeMtlELVX1lCQvTWLH9ytsbgNcVf2HJH+W5B9U1UpVvX7WNcEWeEGSf5rVn+R+afLrn8y6KNikZyT5dFX9eZJTWV0DZ8t1gO3n7yb5bFV9Ocnnk/xxd/+XGde0cOb2MQIAAADzZm47cAAAAPNGgAMAABgJAQ4AAGAkBDgAAICREOAAAABGQoADAAAYCQEOAABgJAQ4AACAkfj/xkVBvO7zEZ4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(np.sqrt(Bbar_array)/1e3, weights=weights_np, bins=100,density=True);\n",
    "fig=plt.gcf()\n",
    "fig.set_size_inches(15,5);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3UAAAEvCAYAAADihOiYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAY+ElEQVR4nO3df6xfZ30f8PendpNqMNFAvK4kMQ6tOzUMFoqbbEKjaEAwpIrpBMN03cKWyaMjWytUCTNQQGZMBtZqbMtWMrBEUWn4NTpLMUuzQrVNXagdSKEOzXBSl1yMIMURGYKSOXz2x/2Gfrnc7/Vx7s/j+3pJls85z3O+fq4fn+/3+/bznOdUdwcAAIBx+oH1bgAAAACPn1AHAAAwYkIdAADAiAl1AAAAIybUAQAAjJhQBwAAMGJb17sBC1188cW9Y8eO9W4GAADAurjrrrv+rLu3Da2/4ULdjh07cuzYsfVuBgAAwLqoqj89l/qmXwIAAIyYUAcAADBiQh0AAMCICXUAAAAjJtQBAACMmFAHAAAwYkIdAADAiAl1AAAAIybUAQAAjJhQBwAAMGJCHQAAwIhtXe8GjMWO/bctWX7y4LVr1BIAAIC/YKQOAABgxIQ6AACAETP9coWYngkAAKwHI3UAAAAjJtQBAACMmFAHAAAwYkIdAADAiAl1AAAAIybUAQAAjJhQBwAAMGJCHQAAwIgJdQAAACMm1AEAAIyYUAcAADBiQh0AAMCICXUAAAAjJtQBAACMmFAHAAAwYkIdAADAiAl1AAAAIzYo1FXV7qq6t6pOVNX+RcpfU1Wfq6q7q+p/VdUVU2VvmJx3b1W9eCUbDwAAsNmdNdRV1ZYkNyd5SZIrkrxqOrRNfKC7n9ndVyZ5R5Jfm5x7RZK9SZ6RZHeS/zh5PQAAAFbAkJG6q5Kc6O77u/uRJLcm2TNdobsfntp9QpKebO9Jcmt3f7u7/yTJicnrAQAAsAK2DqhzSZIHpvbnkly9sFJVvTbJ65JckOTvTJ1754JzL3lcLQUAAOD7DBmpq0WO9fcd6L65u38syeuTvOlczq2qfVV1rKqOPfjggwOaBAAAQDIs1M0luWxq/9Ikp5aof2uSl53Lud19S3fv6u5d27ZtG9AkAAAAkmHTL48m2VlVlyf5UuYXPvn56QpVtbO7vzDZvTbJY9uHk3ygqn4tyVOT7EzyByvR8LHZsf+2mWUnD167hi0BAADOJ2cNdd19pqpuTHJ7ki1JDnX38ao6kORYdx9OcmNVvTDJ/0vyUJLrJ+cer6oPJbknyZkkr+3uR1fpZwEAANh0hozUpbuPJDmy4NhNU9u/tMS5b0vytsfbQAAAAGYb9PBxAAAANiahDgAAYMSEOgAAgBET6gAAAEZMqAMAABgxoQ4AAGDEhDoAAIARE+oAAABGTKgDAAAYMaEOAABgxIQ6AACAERPqAAAARkyoAwAAGDGhDgAAYMSEOgAAgBET6gAAAEZMqAMAABgxoQ4AAGDEhDoAAIAR27reDSDZsf+2JctPHrx2jVoCAACMjZE6AACAERPqAAAARkyoAwAAGDGhDgAAYMSEOgAAgBET6gAAAEZMqAMAABgxoQ4AAGDEhDoAAIARE+oAAABGbFCoq6rdVXVvVZ2oqv2LlL+uqu6pqs9W1e9W1dOmyh6tqrsnvw6vZOMBAAA2u61nq1BVW5LcnORFSeaSHK2qw919z1S1zyTZ1d3frKpfTPKOJK+clH2ru69c4XYDAACQYSN1VyU50d33d/cjSW5Nsme6Qnd/sru/Odm9M8mlK9tMAAAAFjMk1F2S5IGp/bnJsVluSPLxqf0fqqpjVXVnVb3scbQRAACAGc46/TJJLXKsF61Y9QtJdiX5manD27v7VFU9Pcknqupz3X3fgvP2JdmXJNu3bx/UcAAAAIaN1M0luWxq/9IkpxZWqqoXJnljkuu6+9uPHe/uU5Pf70/ye0mevfDc7r6lu3d1965t27ad0w8AAACwmQ0JdUeT7Kyqy6vqgiR7k3zPKpZV9ewk7858oPvq1PGLqurCyfbFSZ6bZHqBFQAAAJbhrNMvu/tMVd2Y5PYkW5Ic6u7jVXUgybHuPpzknUmemOTDVZUkX+zu65L8ZJJ3V9V3Mh8gDy5YNZMBduy/bcnykwevXaOWAAAAG82Qe+rS3UeSHFlw7Kap7RfOOO/3kzxzOQ0EAABgtkEPHwcAAGBjEuoAAABGTKgDAAAYMaEOAABgxIQ6AACAERPqAAAARkyoAwAAGDGhDgAAYMSEOgAAgBET6gAAAEZMqAMAABgxoQ4AAGDEhDoAAIAR27reDWD5duy/bWbZyYPXrmFLAACAtWakDgAAYMSEOgAAgBET6gAAAEZMqAMAABgxoQ4AAGDEhDoAAIARE+oAAABGTKgDAAAYMaEOAABgxIQ6AACAERPqAAAARkyoAwAAGLGt690AVteO/bctWX7y4LVr1BIAAGA1GKkDAAAYMaEOAABgxIQ6AACAERsU6qpqd1XdW1Unqmr/IuWvq6p7quqzVfW7VfW0qbLrq+oLk1/Xr2TjAQAANruzhrqq2pLk5iQvSXJFkldV1RULqn0mya7uflaSjyR5x+TcJyd5c5Krk1yV5M1VddHKNR8AAGBzGzJSd1WSE919f3c/kuTWJHumK3T3J7v7m5PdO5NcOtl+cZI7uvt0dz+U5I4ku1em6QAAAAwJdZckeWBqf25ybJYbknz8cZ4LAADAORjynLpa5FgvWrHqF5LsSvIz53JuVe1Lsi9Jtm/fPqBJAAAAJMNG6uaSXDa1f2mSUwsrVdULk7wxyXXd/e1zObe7b+nuXd29a9u2bUPbDgAAsOkNCXVHk+ysqsur6oIke5Mcnq5QVc9O8u7MB7qvThXdnuSaqrposkDKNZNjAAAArICzTr/s7jNVdWPmw9iWJIe6+3hVHUhyrLsPJ3lnkicm+XBVJckXu/u67j5dVW/NfDBMkgPdfXpVfhIAAIBNaMg9denuI0mOLDh209T2C5c491CSQ4+3gQAAAMw2KNRx/tqx/7aZZScPXruGLQEAAB6PIffUAQAAsEEJdQAAACMm1AEAAIyYUAcAADBiQh0AAMCICXUAAAAjJtQBAACMmFAHAAAwYh4+zkxLPZg88XByAADYCIzUAQAAjJhQBwAAMGJCHQAAwIgJdQAAACMm1AEAAIyYUAcAADBiQh0AAMCIeU4dj5vn2AEAwPozUgcAADBiQh0AAMCICXUAAAAjJtQBAACMmFAHAAAwYkIdAADAiAl1AAAAI+Y5dayapZ5j5xl2AACwMozUAQAAjJhQBwAAMGJCHQAAwIgJdQAAACM2KNRV1e6qureqTlTV/kXKn1dVn66qM1X18gVlj1bV3ZNfh1eq4QAAAAxY/bKqtiS5OcmLkswlOVpVh7v7nqlqX0zy6iS/sshLfKu7r1yBtgIAALDAkEcaXJXkRHffnyRVdWuSPUm+G+q6++Sk7Dur0EYAAABmGDL98pIkD0ztz02ODfVDVXWsqu6sqpedU+sAAABY0pCRulrkWJ/Dn7G9u09V1dOTfKKqPtfd933PH1C1L8m+JNm+ffs5vDQAAMDmNiTUzSW5bGr/0iSnhv4B3X1q8vv9VfV7SZ6d5L4FdW5JckuS7Nq161wCIyO1Y/9tS5afPHjtGrUEAADGbcj0y6NJdlbV5VV1QZK9SQatYllVF1XVhZPti5M8N1P34gEAALA8Zw113X0myY1Jbk/y+SQf6u7jVXWgqq5Lkqr66aqaS/KKJO+uquOT038yybGq+sMkn0xycMGqmQAAACzDkOmX6e4jSY4sOHbT1PbRzE/LXHje7yd55jLbCAAAwAyDHj4OAADAxiTUAQAAjNig6Zew1pZaHdPKmAAA8BeM1AEAAIyYUAcAADBiQh0AAMCICXUAAAAjJtQBAACMmFAHAAAwYh5pwOgs9biDxCMPAADYXIzUAQAAjJhQBwAAMGJCHQAAwIgJdQAAACMm1AEAAIyY1S8571gdEwCAzcRIHQAAwIgJdQAAACMm1AEAAIyYUAcAADBiFkph01lqIRWLqAAAMDZG6gAAAEZMqAMAABgxoQ4AAGDEhDoAAIARE+oAAABGzOqXMGWplTETq2MCALDxGKkDAAAYMaEOAABgxIQ6AACAERt0T11V7U7yriRbkrynuw8uKH9ekn+b5FlJ9nb3R6bKrk/ypsnuv+ru961Ew2E9LHXPnfvtAABYD2cdqauqLUluTvKSJFckeVVVXbGg2heTvDrJBxac++Qkb05ydZKrkry5qi5afrMBAABIhk2/vCrJie6+v7sfSXJrkj3TFbr7ZHd/Nsl3Fpz74iR3dPfp7n4oyR1Jdq9AuwEAAMiwUHdJkgem9ucmx4ZYzrkAAACcxZBQV4sc64GvP+jcqtpXVceq6tiDDz448KUBAAAYEurmklw2tX9pklMDX3/Qud19S3fv6u5d27ZtG/jSAAAADFn98miSnVV1eZIvJdmb5OcHvv7tSf711OIo1yR5wzm3EkZgqZUxE6tjAgCwOs46UtfdZ5LcmPmA9vkkH+ru41V1oKquS5Kq+umqmkvyiiTvrqrjk3NPJ3lr5oPh0SQHJscAAABYAYOeU9fdR5IcWXDspqnto5mfWrnYuYeSHFpGGwEAAJhhUKgDls/0TAAAVsOQhVIAAADYoIQ6AACAETP9EjaIpaZnmpoJAMAsRuoAAABGTKgDAAAYMaEOAABgxIQ6AACAERPqAAAARszqlzACHlwOAMAsRuoAAABGzEgdnAeM5AEAbF5G6gAAAEZMqAMAABgx0y9hE1hqeqapmQAA42akDgAAYMSEOgAAgBET6gAAAEbMPXWwyXkcAgDAuBmpAwAAGDEjdcCSrJwJALCxGakDAAAYMaEOAABgxIQ6AACAEXNPHfC4WTkTAGD9GakDAAAYMSN1wKoxkgcAsPqM1AEAAIyYkTpg3XgGHgDA8hmpAwAAGLFBoa6qdlfVvVV1oqr2L1J+YVV9cFL+qaraMTm+o6q+VVV3T379+so2HwAAYHM76/TLqtqS5OYkL0oyl+RoVR3u7numqt2Q5KHu/vGq2pvk7UleOSm7r7uvXOF2A+c5i6wAAAwzZKTuqiQnuvv+7n4kya1J9iyosyfJ+ybbH0nygqqqlWsmAAAAixmyUMolSR6Y2p9LcvWsOt19pqq+nuQpk7LLq+ozSR5O8qbu/p/LazKARVYAAB4zJNQtNuLWA+t8Ocn27v5aVT0nyW9X1TO6++HvOblqX5J9SbJ9+/YBTQIAACAZFurmklw2tX9pklMz6sxV1dYkT0pyurs7ybeTpLvvqqr7kvxEkmPTJ3f3LUluSZJdu3YtDIwA58T9eADAZjLknrqjSXZW1eVVdUGSvUkOL6hzOMn1k+2XJ/lEd3dVbZsstJKqenqSnUnuX5mmAwAAcNaRusk9cjcmuT3JliSHuvt4VR1Icqy7Dyd5b5L3V9WJJKczH/yS5HlJDlTVmSSPJnlNd59ejR8EYCgjeQDA+WTI9Mt095EkRxYcu2lq+8+TvGKR8z6a5KPLbCMAAAAzDAp1AJuJlTUBgDEZck8dAAAAG5RQBwAAMGKmXwKcA4usAAAbjVAHsILcjwcArDWhDmCNGOUDAFaDe+oAAABGzEgdwAZhJA8AeDyEOoCRcL8eALAY0y8BAABGzEgdwHnA1E0A2LyEOoBNwNRNADh/CXUAm5xRPgAYN6EOgCUJfQCwsVkoBQAAYMSM1AGwLGcbyVuKUT4AWD4jdQAAACNmpA6AdeN+PQBYPqEOgA3L1E4AODvTLwEAAEbMSB0A5yUPXAdgsxDqANh03MsHwPlEqAOABYQ+AMZEqAOAc2QBFwA2EqEOANaQQAjAShPqAGAkBEIAFiPUAcAmIBACnL+EOgBgSRaOAdjYhDoAYFmMAgKsL6EOAFg3ywmEZyMwApvFoFBXVbuTvCvJliTv6e6DC8ovTPIbSZ6T5GtJXtndJydlb0hyQ5JHk/yL7r59xVoPADCDwAhsFmcNdVW1JcnNSV6UZC7J0ao63N33TFW7IclD3f3jVbU3yduTvLKqrkiyN8kzkjw1yX+vqp/o7kdX+gcBAFgrqxkYlyJMAosZMlJ3VZIT3X1/klTVrUn2JJkOdXuSvGWy/ZEk/6GqanL81u7+dpI/qaoTk9f73yvTfACAzWO9wuTZCJuwvoaEukuSPDC1P5fk6ll1uvtMVX09yVMmx+9ccO4lj7u1AABsOBs1bG5UQjArbUioq0WO9cA6Q85NVe1Lsm+y+42qundAu9bSxUn+bL0bQRJ9sZHoi41Ff2wc+mLj0Bcbi/6YqLevdwv0xQYyqy+edi4vMiTUzSW5bGr/0iSnZtSZq6qtSZ6U5PTAc9PdtyS5ZXiz11ZVHevuXevdDvTFRqIvNhb9sXHoi41DX2ws+mPj0Bcbx0r1xQ8MqHM0yc6quryqLsj8wieHF9Q5nOT6yfbLk3yiu3tyfG9VXVhVlyfZmeQPlttoAAAA5p11pG5yj9yNSW7P/CMNDnX38ao6kORYdx9O8t4k758shHI688Evk3ofyvyiKmeSvNbKlwAAACtn0HPquvtIkiMLjt00tf3nSV4x49y3JXnbMtq4EWzYqaGbkL7YOPTFxqI/Ng59sXHoi41Ff2wc+mLjWJG+qPlZkgAAAIzRkHvqAAAA2KCEuomq2l1V91bViarav0j5hVX1wUn5p6pqx9q3cnOoqsuq6pNV9fmqOl5Vv7RInedX1der6u7Jr5sWey2Wr6pOVtXnJn/PxxYpr6r6d5Nr47NV9VPr0c7zXVX9tal/73dX1cNV9csL6rguVlFVHaqqr1bVH00de3JV3VFVX5j8ftGMc6+f1PlCVV2/WB2Gm9EX76yqP568D32sqn54xrlLvqdx7mb0x1uq6ktT70cvnXHukt+/ODcz+uKDU/1wsqrunnGua2MFzfo+u1qfG6ZfJqmqLUn+T5IXZf4xDEeTvKq775mq88+SPKu7X1NVe5P8XHe/cl0afJ6rqh9N8qPd/emq+stJ7krysgX98fwkv9LdP7tOzdw0qupkkl3dvejzbCYf1P88yUuTXJ3kXd199dq1cPOZvGd9KcnV3f2nU8efH9fFqqmq5yX5RpLf6O6/Pjn2jiSnu/vg5AvpRd39+gXnPTnJsSS7Mv+s1ruSPKe7H1rTH+A8MqMvrsn86ttnquafArawLyb1TmaJ9zTO3Yz+eEuSb3T3v1nivLN+/+LcLNYXC8p/NcnXu/vAImUn49pYMbO+zyZ5dVbhc8NI3byrkpzo7vu7+5EktybZs6DOniTvm2x/JMkLqmqxh6uzTN395e7+9GT7/yb5fJJL1rdVLGFP5j88urvvTPLDkzcyVs8Lktw3HehYfd39PzK/wvO06c+G92X+A3uhFye5o7tPTz6Q70iye9Uaugks1hfd/TvdfWaye2fmn43LGphxbQwx5PsX52Cpvph8b/17SX5rTRu1SS3xfXZVPjeEunmXJHlgan8u3x8ivltn8qHx9SRPWZPWbWI1P8312Uk+tUjx36qqP6yqj1fVM9a0YZtLJ/mdqrqrqvYtUj7k+mFl7c3sD2XXxdr6ke7+cjL/AZ7kryxSxzWy9v5xko/PKDvbexor58bJdNhDM6aYuTbW1t9O8pXu/sKMctfGKlnwfXZVPjeEunmLjbgtnJc6pA4rqKqemOSjSX65ux9eUPzpJE/r7r+R5N8n+e21bt8m8tzu/qkkL0ny2snUjmmujTVUVRckuS7Jhxcpdl1sTK6RNVRVb8z8s3F/c0aVs72nsTL+U5IfS3Jlki8n+dVF6rg21tarsvQonWtjFZzl++zM0xY5tuS1IdTNm0ty2dT+pUlOzapTVVuTPCmPb6oBA1TVD2b+AvjN7v4vC8u7++Hu/sZk+0iSH6yqi9e4mZtCd5+a/P7VJB/L/HSZaUOuH1bOS5J8uru/srDAdbEuvvLYdOPJ719dpI5rZI1MFhP42SR/v2csGjDgPY0V0N1f6e5Hu/s7Sf5zFv97dm2skcl317+b5IOz6rg2Vt6M77Or8rkh1M07mmRnVV0++V/wvUkOL6hzOMljK8+8PPM3Y/vfpFUwmfP93iSf7+5fm1Hnrz52T2NVXZX5f8tfW7tWbg5V9YTJzb2pqickuSbJHy2odjjJP6x5fzPzN2B/eY2bupnM/J9W18W6mP5suD7Jf12kzu1JrqmqiyZT0K6ZHGMFVdXuJK9Pcl13f3NGnSHvaayABfdW/1wW/3se8v2LlfHCJH/c3XOLFbo2Vt4S32dX5XNj6/KbPH6TlbJuzPxf1pYkh7r7eFUdSHKsuw9nvlPeX1UnMj9Ct3f9Wnzee26Sf5Dkc1PL7v7LJNuTpLt/PfPB+her6kySbyXZK2Svih9J8rFJTtia5APd/d+q6jXJd/viSOZXvjyR5JtJ/tE6tfW8V1V/KfOrxP3TqWPTfeG6WEVV9VtJnp/k4qqaS/LmJAeTfKiqbkjyxSSvmNTdleQ13f1Puvt0Vb01819gk+RAd5vpsQwz+uINSS5McsfkPevOyYrVT03ynu5+aWa8p63Dj3BemdEfz6+qKzM/ZexkJu9b0/0x6/vXOvwI543F+qK735tF7sV2bay6Wd9nV+VzwyMNAAAARsz0SwAAgBET6gAAAEZMqAMAABgxoQ4AAGDEhDoAAIARE+oAAABGTKgDAAAYMaEOAABgxP4/badQVrXL6YsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(q2_array/1e6, weights=weights_np, bins=100,density=True);\n",
    "fig=plt.gcf()\n",
    "fig.set_size_inches(15,5);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3IAAAE6CAYAAABTQ2DJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd3hUZf7+8fuTRug19AChN0EkAgKCKAgoiquoWLAsLqJg15+4xa/rurvYKyqo6OqqqNhQUESlCgKhS69CpEU60lKe3x8Z3DEGCGQyZ87k/bquuZhzznNm7sTI8MnTzDknAAAAAIB/xHgdAAAAAABwcijkAAAAAMBnKOQAAAAAwGco5AAAAADAZyjkAAAAAMBnKOQAAAAAwGc8LeTMbLSZbTezHwrQ9mkzWxh4rDKz3eHICAAAAACRxrzcR87MukjaL+lN51zLk7jvNkltnHN/LLJwAAAAABChPO2Rc85Nk7Qz+JyZNTCzL81snplNN7Om+dx6laR3wxISAAAAACJMnNcB8jFK0mDn3Gozay/pRUnnHr1oZnUlpUj61qN8AAAAAOCpiCrkzKyMpI6SPjCzo6dL5GnWX9JY51x2OLMBAAAAQKSIqEJOuUM9dzvnTj9Om/6ShoQpDwAAAABEnIjafsA5t1fSejO7XJIsV+uj182siaSKkmZ5FBEAAAAAPOf19gPvKrcoa2Jm6WY2UNI1kgaa2SJJSyX1DbrlKkljnJdLbQIAAACAxzzdfgAAAD8ws16SnpUUK+lV59zwPNcHK3fYf7Zyt9UZ5JxbZmb1JC2XtDLQ9Hvn3OBw5QYARC8KOQAAjsPMYiWtktRDUrqkuZKucs4tC2pTLjA9QGZ2saRbnXO9AoXc5yezVyoAAAURUXPkAACIQO0krXHOrXPOHZE0Rr8d9n90jvdRpSXxW1IAQJHybNXKKlWquHr16nn19gCAMJo3b97Pzrkkr3OcolqSNgUdp0tqn7eRmQ2RdLekBAXtfyopxcwWSNor6a/Ouen5vYmZDZI0SJJKly7dtmnTpqFJDwCIWIX5fPSskKtXr57S0tK8ensAQBiZ2Y9eZygEy+fc73rcnHMjJI0ws6sl/VXS9ZK2SKrjnNthZm0lfWJmLfL04B29f5SkUZKUmprq+IwEgOhXmM9HhlYCAHB86ZKSg45rS9p8nPZjJF0iSc65w865HYHn8yStldS4iHICAIoRCjkAAI5vrqRGZpZiZgmS+ksaF9zAzBoFHV4oaXXgfFJgsRSZWX1JjSStC0tqAEBU82xoJQAAfuCcyzKzoZImKnf7gdHOuaVm9rCkNOfcOElDzay7pExJu5Q7rFKSukh62MyylLs1wWDn3M7wfxUAgGhDIQcAwAk45yZImpDn3INBz+84xn0fSvqwaNMBAIojhlYCAAAAgM9QyAEAAACAz1DIAQAAAIDPUMgBAAAAgM9QyAEAAACAzxSLVSuPZOVoxda92rDjgI5k5SgxPkY1yieqfpUyqlg6wet4AAAAAHBSorqQ23MgUy9MXq3309K152Bmvm3qVS6l1HqVdG7TqurSOEllSkT1twQAAABAFDhh1WJmyZLelFRdUo6kUc65Z/O0OUfSp5LWB0595Jx7OLRRT87yLXs18I252rr3kC5sVVO9WlRXo2plVDI+Vr8cydKW3Ye0fOteLdy4W5OWbdPYeemKjzWd3ShJl55RS92bVVNifKyXXwIAAKes3rDxvz7fMPxCD5MAAIpCQbqfsiTd45ybb2ZlJc0zs0nOuWV52k13zvUJfcSTl77rgAa8NltxMTH6ZEgntapd4XdtmlYvp25Nq0qSsrJzNO/HXfp6+TZ9vniLhr6zXeUS43RR65q6tkNdNatRLtxfAgAAAAAc0wkLOefcFklbAs/3mdlySbUk5S3kIoJzTvd+sEiHM3M0ZshZali1zAnviYuNUfv6ldW+fmUN691MM9f+rA/npevD+el6e/ZGtUuppBs71lOP5tUUF8v6MAAAfwnunZPooQOAaHBSE8LMrJ6kNpJm53P5LDNbJGmzpHudc0sLne4UfPnDVn2/bqceuaRlgYq4vGJjcodXnt0oSbsPHNH7aZv05qwfdcvb81WzfKKu61hP17Svo7KJ8UWQHgAAAABOrMDdS2ZWRtKHku50zu3Nc3m+pLrOudaSnpf0yTFeY5CZpZlZWkZGxqlmPq5R09epXuVSuqpdnUK/VoVSCRrUpYGm3tdNr1yXqpSk0hr+xQp1HP6tHp+4Qj/vPxyCxAAAAABwcgpUyJlZvHKLuLedcx/lve6c2+uc2x94PkFSvJlVyafdKOdcqnMuNSkpqZDRf2/p5j1asHG3ru9YT7ExFrLXjY0x9WheTW/f1EGfDe2ssxtV0YtT1qrT8G/14Kc/aNPOAyF7LwAAAAA4kYKsWmmSXpO03Dn31DHaVJe0zTnnzKydcgvEHSFNWgATl25TjEkXt65ZZO9xWu3yevGatlqbsV+jpq7Tu3M26u3ZG/WHNrV027kNVbdy6SJ7bwAAAACQCjZHrpOkAZKWmNnCwLk/S6ojSc65lyX1k3SLmWVJOiipv3POFUHe45q0bJva1q2oymVKFPl7NUgqo0f7tdKdPRrplWnr9fbsH/Xxgp902Rm1dNu5jZRcqVSRZwAAAABQPBVk1coZko47TtE594KkF0IV6lTs+uWIlm/Zq/t6Ngnr+9YoX1IPXtRcg7vW14tT1uqdORv10fyfdHlqbQ3p1lC1K1LQAQAAAAitqFlLf+Gm3ZKkM+pU9OT9q5ZL1EMXt9C0+7rpmvZ19OG8n9TtiSn688dL9NPug55kAgAAABCdoqaQW7Bxl2JMap1c3tMc1csn6u99W2rKfefoyjOT9UHaJnV7fIoeGreUVS4BAAAAhETUFHI/bN6rxtXKqlTCSW2NV2RqViipRy45TVPu66ZLz6ilt77/UV0em6ynJq3SvkOZXscDAAAA4GNRU8itzdivBqewAXhRq1WhpIZf1kpf3dVF3ZpU1XPfrFaXxybr1enrdCgz2+t4AAAAAHwoKgq5w1nZ2rTzgBpUidyl/xskldGIa87QZ0M7q2Wt8npk/HJ1e2KK3pu7UVnZOV7HAwAAAOAjUVHIbdxxQDlOqp8UeT1yeZ1Wu7zeGthe79zUXlXLJer+D5fo/Gem6YslW+TBjg0AAAAAfCgqCrl1P/8iSUqJ4B65vDo2rKJPbu2okQPaKsZMt7w9X31HfKeZa372OhoAAACACBcVhdzWPYckSbUqlvQ4yckxM/VsUV0T7+yix/u10o79R3T1q7N1w+tztHLrPq/jAQAAAIhQUVHIbdt7SHExpkqlEryOckpiY0yXpybrm3u66oHeTTXvx13q/ew03T928a9FKgAAAAAcFSWF3GFVLVtCMTHmdZRCSYyP1c1dG2jafd10Q8cUfbQgXec8MVlPfrVS+w9neR0PAAAAQISIikJu+75Dqlou0esYIVOxdIIevKi5vrn7HHVvVk3Pf7tGXR+brLdmbVAmK1wCAAAAxV5UFHJb9xxStXIlvI4RcnUql9ILV5+hT4Z0UoOkMvrbp0vV8+lpmrh0KytcAgAAAMVYVBRy2/YeUtWy0dMjl9fpyRX03s0d9Mp1qTKTbn5rnq4YOUvzN+7yOhoAAAAAD/i+kMvOcdp7KEuVSvtzoZOCMjP1aF5NE+/son/+oaXW/3xAl744U0Penq8fd/zidTwAiGpm1svMVprZGjMbls/1wWa2xMwWmtkMM2sedO2BwH0rzaxneJMDAKKV7wu5fYcyJUnlS8Z7nCQ84mJjdE37uppy3zm647xG+nbFdnV/aqoe/myZ9hzI9DoeAEQdM4uVNEJSb0nNJV0VXKgFvOOcO805d7qkxyQ9Fbi3uaT+klpI6iXpxcDrAQBQKL4v5PYcLF6F3FFlSsTprh6NNfW+c3TZGbX1+sz16vrEZL3+3XoWRAGA0GonaY1zbp1z7oikMZL6Bjdwzu0NOiwt6ehE5r6SxjjnDjvn1ktaE3g9AAAKhULO56qWS9Twy1pp/G1nq0XNcvr7Z8vU8+lp+nrZNhZEAYDQqCVpU9BxeuDcb5jZEDNbq9weudtP5l4AAE5W9BRypYpnIXdU85rl9N+B7fXa9amSSTe9maZrX5utZZv3nvhmAMDx5LdJ6e9+U+acG+GcayDpfkl/PZl7JcnMBplZmpmlZWRknHJYAEDxED2FXDHtkQtmZjqvWe6CKA9d1FxLN+/Vhc9P1/1jF2v7vkNexwMAv0qXlBx0XFvS5uO0HyPpkpO91zk3yjmX6pxLTUpKKkRcAEBxQCEXheJjY3RDpxRNvbebBnZK0UcL0nXO41P0wrerdSgz2+t4AOA3cyU1MrMUM0tQ7uIl44IbmFmjoMMLJa0OPB8nqb+ZlTCzFEmNJM0JQ2YAQJSL8zpAYVHIHVv5UvH6a5/muqZDXQ3/Yrme+GqV3pm9Uff3bqqLWtVUTEx+I34AAMGcc1lmNlTSREmxkkY755aa2cOS0pxz4yQNNbPukjIl7ZJ0feDepWb2vqRlkrIkDXHOef4btXrDxv/6fMPwCz1MAgA4Vb4v5PYfylJcjCkxntWcjyWlSmmNHJCq79ft0CPjl+mOMQs1+rsN+tuFzZRar5LX8QAg4jnnJkiakOfcg0HP7zjOvf+U9M+iSwcAKI58P7TywJFslUygiCuIDvUra9yQznri8tbauueg+r08S0Pema9NOw94HQ0AAADASfB9IXcoM1ulKOQKLCbG1K9tbU2+N3dD8W+Wb9N5T07Vv79Yrr2H2FAcAAAA8APfF3IHjmSrJMMqT1qphNwNxSffe476tK6hkVPXqdvjU/TO7I3KzmH/OQAAACCS+b6QO5iZrZIJvp/q55ka5UvqqStO17ihndQgqYz+/PES9Xl+hmat3eF1NAAAAADH4P9C7ki2Ssb7/svwXKvaFfTezR30wtVttPdgpq565Xvd8t95zJ8DAAAAIpDvK6CDmdkqRY9cSJiZ+rSqqW/u6aq7ezTWlJUZOu+pqXp84gr9cjjL63gAAAAAAnxfyB04ks3WAyGWGB+r289rpG/v7aoLWlbXiMlr1e2JKfpwXrpymD8HAAAAeM73hRyrVhadGuVL6pn+bfThLR1Vo3yi7vlgkf7w0kzN37jL62gAAABAseb7Qu7AkSxWrSxibetW1Me3dtKTl7fWlt0HdemLM3XXewu1dc8hr6MBAAAAxZLvC7mDbAgeFjExpssC+88N6dZA45dsUbcnpui5b1brUGa21/EAAACAYsX/hVwmhVw4lS4Rp/t6NtU3d3fVOU2S9NSkVTrvyakav3iLnGP+HAAAABAOvi7ksrJzlJntlBhHIRduyZVK6aVr2+rdP3VQ2cQ4DXlnvq4c+b1++GmP19EAAACAqOfrQi4zO7cHqAT7yHnmrAaVNf72s/WvP5ymNRn7ddELMzTsw8XK2HfY62gAAABA1PJ1BXQkO0eSFB/r6y/D92JjTFe3r6PJ956jgZ1SNHZeuro9MUWjpq3Vkawcr+MBAAAAUcfXFVBmoJBLiDWPk0CSypeM11/7NNfEu7qoXUol/WvCCp3/9FR9vWwb8+cAAACAEIqKQo4eucjSIKmMRt9wpt648UzFxphuejNN178+V2u27/c6GgAAABAVfF0BZWbl9vJQyEWmc5pU1Zd3dtHf+jTXgo271OuZafrXhOXadyjT62gAAACAr/m6Avp1jlycr7+MqBYfG6OBnVM0+d5zdNkZtfXK9HXq9sRUjZ2XrpwchlsCAAAAp+KEFZCZJZvZZDNbbmZLzeyOfNqYmT1nZmvMbLGZnVE0cX/r6EIazJGLfFXKlNCj/Vrpk1s7qXbFkrr3g0W69KWZWrRpt9fRAAAAAN8pSFdWlqR7nHPNJHWQNMTMmudp01tSo8BjkKSXQpryGJgj5z+tkyvoo1s66snLWyt910Fd8uJ3un/sYv28n+0KAAAAgIKKO1ED59wWSVsCz/eZ2XJJtSQtC2rWV9KbLndpwu/NrIKZ1QjcW2Qo5PwpJsZ0WdvaOr9FNT3/7RqNnrFeE37Yoju7N9Z1Z9XlvycAhFG9YeN/c7xh+IUeJQEAnIyT+hezmdWT1EbS7DyXaknaFHScHjhXpNhHzt/KJsbrzxc005d3dlGbOhX1j8+X6YJnp2vG6p+9jgYAAABEtAJXQGZWRtKHku50zu3NezmfW363koWZDTKzNDNLy8jIOLmk+cjMzn2LBBY78bWGVcvoPzeeqVeuS9XhrBxd+9psDX5rnjbtPOB1NAAAACAiFagCMrN45RZxbzvnPsqnSbqk5KDj2pI2523knBvlnEt1zqUmJSWdSt7fyPx1sRMKOb8zM/VoXk1f3dVF9/VsoqmrMtT9qal6atIqHTyS7XU8AAAAIKIUZNVKk/SapOXOuaeO0WycpOsCq1d2kLSnqOfHSUFz5OJYtTJaJMbHaki3hvr23q7q2aK6nvtmtbo/NVUTlmxR7hRMAAAAAAXpyuokaYCkc81sYeBxgZkNNrPBgTYTJK2TtEbSK5JuLZq4v8UcuehVo3xJPXdVG703qIPKJsbp1rfn6+pXZmvl1n1eRwMAAAA8V5BVK2co/zlwwW2cpCGhClVQv86Ro5CLWu3rV9bnt3XWu3M36cmvVuqC56ZrQIe6uqt7Y5UvFe91PAAAAMATvq6A2H6geIiLjdGADnU1+Z5zdHW7Onpz1gZ1e3KK3p2zUdk5DLcEAABA8ePrCuh/hRxz5IqDiqUT9I9LWuqz2zqrYVIZPfDREvUdMUPzftzpdTQAAAAgrHxdyB0JrFoZR49csdKiZnm9d3MHPXdVG/2874gue2mW7nl/kTL2HfY6GgAAABAWvq6AsgLD6uiRK37MTBe3rqlv7umqW85poHGLftK5T07R69+tV1agpxYAQsXMepnZSjNbY2bD8rl+t5ktM7PFZvaNmdUNupYdtFjYuPAmBwBEK18XckfnR8UYhVxxVbpEnO7v1VRf3tlFpydX0N8/W6Y+z8/QnPUMtwQQGmYWK2mEpN6Smku6ysya52m2QFKqc66VpLGSHgu6dtA5d3rgcXFYQgMAop6vC7mcQCEXG0MhV9w1SCqjN//YTi9fe4b2HcrSFSNn6a73Fmr73kNeRwPgf+0krXHOrXPOHZE0RlLf4AbOucnOuQOBw+8l1Q5zRgBAMePrQi47sEF0LD1yUO5wy14ta+jru7tqaLeGGr94i859cqpenb7u14VxAOAU1JK0Keg4PXDuWAZK+iLoONHM0szsezO7pCgCAgCKH18Xckd75GLokUOQkgmxurdnE028q4tS61XUI+OXq89zM/T9uh1eRwPgT/l9yOS794mZXSspVdLjQafrOOdSJV0t6Rkza3CMewcFCr60jIyMwmYGAEQ5Xxdy2c4xrBLHlFKltF6/4UyNGtBW+w9nqf+o73XHmAXaxnBLACcnXVJy0HFtSZvzNjKz7pL+Iuli59yvy+g65zYH/lwnaYqkNvm9iXNulHMu1TmXmpSUFLr0J6nesPG/PgAAkcvXhVxWDoUcjs/MdH6L6vr67q66/bxG+uKHrTr3iSl6ZRrDLQEU2FxJjcwsxcwSJPWX9JvVJ82sjaSRyi3itgedr2hmJQLPq0jqJGlZ2JIDAKKWrwu5nBzH/DgUSMmEWN3do7Em3dVFHepX1j8nLFfvZ6dr5pqfvY4GIMI557IkDZU0UdJySe8755aa2cNmdnQVyscllZH0QZ5tBppJSjOzRZImSxrunKOQAwAUWpzXAQojO4cVK3Fy6lYurdduOFNfL9umv3++VFe/Olt9WtXQXy5sphrlS3odD0CEcs5NkDQhz7kHg553P8Z9MyWdVrTpAADFkb975JwTdRxORffm1TTprq66q3tjTVq2Tec9OVUvT12rI1kMtwQAAEDk83Uhl80cORRCYnys7ujeSF/f3VWdGlbR8C9WqPez0zRjNcMtAQAAENn8Xcg5p9gYX38JiADJlUrpletS9foNZyorx+na12br1rfnafPug15HAwAAAPLl7zly2U6x1HEIkW5Nq+qsBpX1yrR1GjFljSavyNDQcxvqT2fXV0IcP2gAAACIHL7+12m2Y9VKhFZifKxuOy93uGWXxlX0+MSV6v3sNFa3BAAAQETxdSGXk+MUwxw5FIHaFUtp5IBUvX7jmcrMdrr61dm6Y8wCbWczcQAAAEQAXxdyuXPkKORQdLo1qaqv7uqiO85rpC+WbNV5T07VG9+tVxabiQMAAMBD/i7k2BAcYZAYH6u7ejTWxLu66PQ6FfTQZ8vUd8R3WrBxl9fRAAAAUEz5v5CjRw5hklKltN78Yzu9eM0Z2rH/iC59aaYe+GiJdh844nU0AAAAFDMUcsBJMDNdcFoNfX1PV93UOUXvp23SuU9O1ftpm5ST47yOBwAAgGLC14VcjnOKYWglPFCmRJz+cmFzfX5bZ6VUKa3/N3axrhg5S8u37PU6GgAAAIoBXxdy9MjBa81qlNMHN5+lx/q10tqM/erz/Aw98vky7T+c5XU0AAAARDF/bwjuxPYD8FxMjOmK1GT1aFZNj01cqVdnrNdnizfrwT4tdMFp1WX0GgPwqXrDxv/meMPwCz1KAgDIy9c9cjk5TnEUcogQFUsn6N+XnqaPbu2oyqVLaMg783Xd6Dla//MvXkcDAABAlPF1IZeVk8P2A4g4Z9SpqHFDO+mhi5pr4cbd6vn0ND01aZUOZWZ7HQ0AAABRwteFXE6OFOPrrwDRKi42Rjd0StE393RV79Oq67lvVqvnM9M0eeV2r6MBAAAgCvi6DMp2LHaCyFa1XKKe7d9G79zUXrExphtfn6tb/jtPm3cf9DoaAAAAfMzfhVwO2w/AHzo2rKIv7jhb9/Vsoskrt6vHU1P12oz1ysrO8ToaAAAAfMjXhVwOPXLwkRJxsRrSraEm3dVV7VIq6R+fL1PfEd9pcfpur6MBAADAZ3xdyGVls2ol/Ce5UimNvuFMvXjNGcrYd1iXjPhOD41bqn2HMr2OBgAAAJ/wdSGX4xhaCX8yM11wWg19fU9XDehQV/+ZtUHdn5qqL5ZskXPO63gAAACIcBRygIfKJcbr731b6uNbO6ly6RK65e35GvifNG3aecDraAAAAIhgPi/k2H4A0eH05AoaN7ST/nphM32/bofOf3qaRk5dq0wWQwEAAEA+fF0GOedkokcO0SEuNkY3nV1fk+7uqs6NqujfX6zQRc/P0Lwfd3kdDQAAABHG34WcJEZWItrUqlBSr1yXqpED2mrPwUz1e3mm/vLxEu05yGIoAAAAyOXvQs7lLhoBRKOeLapr0t1d9cdOKXp3zkad9+RUjVu0mcVQAAAA4PdCzondBxDNypSI09/6NNe4oZ1Vs0Kibn93ga4bPUc/7vjF62gAiqF6w8b/+gAAeOuEhZyZjTaz7Wb2wzGun2Nme8xsYeDxYOhj5i/HiRlyKBZa1iqvj2/tpIcuaq4FG3fr/KenacTkNTqSxWIoAAAAxVFBeuTekNTrBG2mO+dODzweLnysgnFyDK1EsREbY7qhU4q+vrurzmtWVY9PXKkLn5uuOet3eh0NAAAAYXbCQs45N01SRP5LMXeOnNcpgPCqXj5RL17TVqNvSNWBI9m6YuQsPfDRYhZDAQAAKEZCNUfuLDNbZGZfmFmLEL3mCTknth9AsXVu02qadHcXDepSX+/N3aTuT03VhCVbWAwFKAJm1svMVprZGjMbls/1u81smZktNrNvzKxu0LXrzWx14HF9eJMDAKJVKAq5+ZLqOudaS3pe0ifHamhmg8wszczSMjIyCv3GLHaC4q5UQpz+fEEzjRvaWVXLltCtb8/XoLfmaeueQ15HA6KGmcVKGiGpt6Tmkq4ys+Z5mi2QlOqcayVprKTHAvdWkvR/ktpLaifp/8ysYriyAwCiV6ELOefcXufc/sDzCZLizazKMdqOcs6lOudSk5KSCvvWuYudUMgBalmrvD4d0kl/vqCppq/OUPenpuqtWRuUk0PvHBAC7SStcc6tc84dkTRGUt/gBs65yc65A4HD7yXVDjzvKWmSc26nc26XpEk68bxzAABOqNCFnJlVt8CKI2bWLvCaOwr7ugXh5BhaCQTExcZoUJcG+urOrjo9uYL+9ulSXT5yllZv2+d1NMDvaknaFHScHjh3LAMlfXGK9wIAUCAF2X7gXUmzJDUxs3QzG2hmg81scKBJP0k/mNkiSc9J6u/CNEnHOSnG1zvhAaFXp3IpvTWwnZ68vLXWZezXBc9N11OTVulwVrbX0QC/yu83hvl+zpnZtZJSJT1+CveGdPoBACC6xZ2ogXPuqhNcf0HSCyFLdBJyR43RIwfkZWa6rG1tndMkSf/4fJme+2a1xi/erOGXtdKZ9Sp5HQ/wm3RJyUHHtSVtztvIzLpL+oukrs65w0H3npPn3in5vYlzbpSkUZKUmprKuGgAwHH5vD+LxU6A46lcpoSe6d9G//ljOx3KzNHlL8/SXz5eor2H2KoAOAlzJTUysxQzS5DUX9K44AZm1kbSSEkXO+e2B12aKOl8M6sYWOTk/MA5AAAKxdeFHIudAAXTtXGSJt3dRTd1TtG7czaqx1NT9eUPW72OBfiCcy5L0lDlFmDLJb3vnFtqZg+b2cWBZo9LKiPpAzNbaGbjAvfulPQP5RaDcyU9HDgHAEChnHBoZSRzjsVOgIIqlRCnv/ZprotPr6lhHy7R4P/OU88W1fRw35aqVi7R63hARAusyjwhz7kHg553P869oyWNLrp0AIDiyN+FnMTQSuAktapdQZ8O7aTXZqzX05NWqfuTU3V/76a6ul0dxfA/FIACqjds/G+ONwy/0KMkAFA8+XtoZY6TMbYSOGnxsTEa3LWBvrqri1oll9dfP/lBV46apXUZ+72OBgAAgALwdSHnxBw5oDDqVi6t/w5sr8f7tdKqbfvV69npennqWmVl53gdDQAAAMfh70LOiTlyQCGZmS5PTdaku7vo3CZVNfyLFfrDizO1fMter6MBAADgGHxeyDl65IAQqVo2US8PaKsXrzlDW/Yc1EXPz9DTk1bpSBa9cwAAAJHG34WcWOwECLULTquhSXd11UWta+rZb1broudnaNGm3V7HAgAAQPnWWyMAACAASURBVBBfF3I5jsVOgKJQsXSCnr7ydI2+IVV7DmbqDy9+p39NWK5DmdleRwMAAIB8Xsg5NgQHitS5Tavpq7u76Moz62jUtHXq9cw0zV63w+tYAAAAxZ7/CzkWOwGKVLnEeP370tP0zp/aK8dJV476Xn/75AftP5zldTQAAIBiy+cbgrPYCRAuHRtU0Zd3nq0nJq7S6zPX69sV2/WvS09T18ZJXkcDEAGCNwhnc3AAKHq+75FjsRMgfEolxOnBi5pr7OCOKpkQq+tHz9G9HyzSngOZXkcDAAAoVnxdyOU4x9BKwANt61bU57d11tBuDfXxgp/U/empmrh0q9exAAAAig1fF3JsPwB4JzE+Vvf2bKJPh3RSUpkSuvmtebpzzALtPnDE62gAAABRz9+FnBPLVgIea1mrvD4d2kl3dW+szxdvUY+np+mb5du8jgUAABDVfFvIOeckiYGVQASIj43RHd0b6ZMhnVS5dIIG/idN97y/SHsOMncOAACgKPi4kMv9M4YeOSBitKxVXuOGdtZt5zbUJwt/Us+np2nyyu1exwIAAIg6vi3kco72yFHHARElIS5G95zfRB/f2lFlE+N04+tzdf/Yxdp7iN45AACAUPHtPnKBDjkWOwEiVKvaFfT57Z31zNerNXLqWk1fnaFH+7XS2Y3Ydw6IdsF7yknsKwcARSEKeuSo5IBIVSIuVvf3aqoPb+moxIRYDXhtjv788RLtP5zldTQAAABf820hd3SOHIDI16ZORU24/WwN6lJf787ZqJ5PT9PMNT97HQsAAMC3fFvIHcViJ4A/JMbH6s8XNNMHN5+lhLgYXf3qbD346Q/6hd45AACAk+bbQo7FTgB/Sq1XSRNuP1t/7JSit77/URc8N13zftzpdSwAAABf8W0h97/tB7zNAeDklUyI1YMXNdeYP3VQdo7T5S/P0mNfrtCRrByvowEAAPiCbwu5X3vk2BIc8K329SvrizvOVr+2tfXilLW6ZMR3WrVtn9exAAAAIp5vC7mja50wtBLwt7KJ8XqsX2uNGtBW2/YeUp/nZ+jV6euUk8OKRgAAAMfi30Iu8G88th8AosP5Lapr4l1d1KVRkh4Zv1zXvDpbP+0+6HUsAACAiOTjQu7o0EoA0aJKmRJ65bq2euyyVlqcvlu9np6mD+el//r/OwB/qjds/K8PAEBo+LiQy/2TxU6A6GJmuuLMZH1xRxc1rVFW93ywSLe+PV87fznidTQAAICI4dtC7n/bD1DJAdGoTuVSGjPoLA3r3VRfL9+mns9M0+QV272OBQAAEBF8W8ix2AkQ/WJjTIO7NtCnQzqrcukE3fjGXP354yVsIg4AAIo9/xZyLHYCFBvNa5bTp0M76eYu9fXunI3q8/wMLU7f7XUsFCNm1svMVprZGjMbls/1LmY238yyzKxfnmvZZrYw8BgXvtQAgGjm40KOxU6A4qREXKweuKCZ3v1TBx3KzNalL87US1PWsk0BipyZxUoaIam3pOaSrjKz5nmabZR0g6R38nmJg8650wOPi4s0LACg2PBvIRf4M4YeOaBY6RDYRPz8FtX06JcrdO1rs7V1zyGvYyG6tZO0xjm3zjl3RNIYSX2DGzjnNjjnFkvK8SIgAKD4ifM6wKn632InHgcBEHYVSiVoxNVn6IO0dP3fuKXq9ew0Db+0lXq1rO51NESnWpI2BR2nS2p/EvcnmlmapCxJw51zn+TXyMwGSRokSXXq1DnFqJEv7xYEG4Zf6FESAPA3//bIHZ0j520MAB45uk3B+Ns7K7liKQ3+7zw98NESHTjCQigIufw+ak5mTG8d51yqpKslPWNmDfJr5Jwb5ZxLdc6lJiUlnUpOAEAx4t9CLvAnPXJA8VY/qYw+vKWjbu5aX2Pm5i6E8sNPe7yOheiSLik56Li2pM0Fvdk5tznw5zpJUyS1CWU4AEDx5NtC7iijTw4o9hLiYvRA72Z6e2B7/XI4S3948Tu9Mm0dC6EgVOZKamRmKWaWIKm/pAKtPmlmFc2sROB5FUmdJC0rsqQAgGLjhIWcmY02s+1m9sMxrpuZPRdYknmxmZ0R+pi/d3TVSgA4qmPDKvryji7q1qSq/jlhua5/fY627WUhFBSOcy5L0lBJEyUtl/S+c26pmT1sZhdLkpmdaWbpki6XNNLMlgZubyYpzcwWSZqs3DlyFHIAgEIryGInb0h6QdKbx7jeW1KjwKO9pJd0cpPAC4cOOQBBKpZO0MgBbfXunE16+POl6v3sdD15RWt1a1LV62jwMefcBEkT8px7MOj5XOUOucx730xJpxV5QABAsXPCHjnn3DRJO4/TpK+kN12u7yVVMLMaoQoIACfLzHR1+zr6/LbOqlq2hG58fa7+/cVyZWazMjwAAIgOoZgjl9+yzLVC8LrHxchKACfSsGpZfTKkk65uX0cjp67TFSNnKX3XAa9jAQAAFFooCrkCL8tsZoPMLM3M0jIyMkLw1oysBHB8ifGx+tcfTtPzV7XR6m37deFzM/TV0q1exwIQUG/Y+F8fAICCC0UhV+BlmdkjB4BXLmpdU+Nv76w6lUpp0Fvz9PfPlupwVrbXsQAAAE5JKAq5cZKuC6xe2UHSHufclhC8LgCEVN3KpTX2lrN0Y6d6ev27Der30iz9uOMXr2MBAACctIJsP/CupFmSmphZupkNNLPBZjY40GSCpHWS1kh6RdKtRZY2/3zhfDsAPlciLlb/d1ELjRzQVj/u+EUXPjdDny8u8N7OAAAAEeGE2w845646wXUnaUjIEgFAGPRsUV0tapbTbe8u0NB3Fmjm2h16sE9zJcbHeh0NAADghEIxtNITrFoJoLBqVyyl928+Szd3ra93Zm9Uv5dnatNOVrUEAACRz7eF3FEMrARQGPGxMXqgdzO9cl2qftxxQH2en6HJK7Z7HQsAAOC4fFvIufx3OACAU9KjeTV9fltn1apQUje+MVdPfbVS2Tn8PQOEE1sRAEDB+baQO4q1TgCESt3KpfXRrR11edvaeu7bNbrh9Tna+csRr2MBAAD8ju8LOQAIpcT4WD1+eWs9etlpmr1+p/o8N10LNu7yOhYAAMBv+LaQY7ETAEXpyjPr6KNbOiomxnTFyFl6a9YGOf7iAQAAEcK3hdxRDK0EUFRa1iqv8bedrbMbJelvny7Vne8t1MEj2V7HAgAAOPE+cgBQnJUvFa9Xr0vVi1PW6MlJq7R6236NHNBWyZVKeR0NiGp5FzzZMPxCj5IAQGTybY8cA5wAhEtMjGnouY00+oYztWnXAfUd8Z1mrv3Z61gAAKAY820hd5SxkxyAMOnWpKrGDe2sSqUTNOC1OXr9u/XMmwMAAJ7wfSEHAOGUUqW0Pr61o7o1qaq/f7ZM941drEOZzJsDAADh5dtCjt+CA/BK2cR4jRrQVnec10hj56XrylHfa+ueQ17HAgAAxYhvC7mjWLUSgBdiYkx39WiskQPaas22ferz/AzN+3Gn17EAAEAx4ftCDgC81LNFdX08pJPKlIhV/1Hf64O0TV5HAqJSvWHjf30AAHxcyDGwEkCkaFytrD4d0lntUyrrvrGL9eiXK5STw99SAACg6Pi2kAOASFK+VLxev/FMXd2+jl6asla3vj1fB45keR0LAABEKd8Wcqx1AiDSxMfG6J+XtNTf+jTXxGVbdeXI77VtL4ugAACA0PNtIXeUsdoJgAhiZhrYOUWvXpeqtRn71feF7/TDT3u8jgUAAKKM7ws5AIhE5zWrprGDOyrGpMtfnqWvlm71OhIQNYIXPmHxEwDFlY8LOcZWAohszWuW0ydDOqlxtTK6+b/z9Or0dV5HAgAAUcLHhVwuBlYCiGRVyyVqzKCz1LN5dT0yfrke+XwZK1oCAIBC830hBwCRrmRCrEZcc4auP6uuXp2xXne8t1CHs7K9jgUAAHzMt4Ucq1YC8JPYGNNDF7fQsN5N9dmizbp+9BztOZjpdSwUkJn1MrOVZrbGzIblc72Lmc03sywz65fn2vVmtjrwuD58qQEA0SzO6wCFxaKVAPzCzDS4awNVL5eo+8Yu0pUjZ+mNG9upevlEr6PhOMwsVtIIST0kpUuaa2bjnHPLgpptlHSDpHvz3FtJ0v9JSlXu5O55gXt3hSN7cRG84MmG4Rd6mAQAwse3PXIA4FeXtKml129op/RdB3Xpi99p1bZ9XkfC8bWTtMY5t845d0TSGEl9gxs45zY45xZLyslzb09Jk5xzOwPF2yRJvcIRGgAQ3XxbyDGyEoCfdW5URe/d3EGZOU79Xpqp+RvpoIlgtSRtCjpOD5wr6nsBADgm3xZyRxnrVgLwqRY1y+ujWzqqYukEXfvqbM1c87PXkZC//D5oCvr7xALfa2aDzCzNzNIyMjIKHA4AUDz5vpADAD9LrlRKH9x8lpIrltINb8zVpGXbvI6E30uXlBx0XFvS5lDf65wb5ZxLdc6lJiUlnVJQAEDx4dtCjlUrAUSLquUS9d7NHdSsRjkN/u88fbrwJ68j4bfmSmpkZilmliCpv6RxBbx3oqTzzayimVWUdH7gHIpIvWHjf/MAgGjl20LuKFatBBANKpRK0Ns3tVe7epV053sL9d/vf/Q6EgKcc1mShiq3AFsu6X3n3FIze9jMLpYkMzvTzNIlXS5ppJktDdy7U9I/lFsMzpX0cOAcAACF4tvtBxzLnQCIMmVKxOn1G8/U0Hfm66+f/KBfDmfp5q4NvI4FSc65CZIm5Dn3YNDzucodNpnfvaMljS7SgACAYse3hdxRdMgBiCaJ8bF66dq2uvv9Rfr3FyuU46RbzqGYA04Ve8wBiFa+L+QAINrEx8bo6StaK8akR79cIYliDgAA/JZvCzkWOwEQzeJiY/Tk5a0lUcwBAIDf820hdxSLnQCIVhRzAADgWHxfyAFANKOYA0In73YEzJkD4Ge+LeQYWgmguDhazDmXW8yVKRGrAWfV8zoWAADwkG8Luf9hbCWA6BcXG6Mnr2itA0ey9bdPl6pMYpz+0Cbf1e4BAEAx4PsNwQGguIiPjdELV7fRWfUr694PFmvSsm1eRwIAAB4pUCFnZr3MbKWZrTGzYflcv8HMMsxsYeBxU+ij/hYbggMojhLjY/XK9alqWau8hrwzXzPX/Ox1JAAA4IETFnJmFitphKTekppLusrMmufT9D3n3OmBx6shznmcfOF6JwCIDGVKxOk/N56plMqlddObaVq0abfXkQBfqjds/K8PAPCbgvTItZO0xjm3zjl3RNIYSX2LNtaJsdgJgOKsQqkEvTWwnSqXSdDA/8zVpp0HvI4EAADCqCCFXC1Jm4KO0wPn8rrMzBab2VgzSw5JOgDAMVUtl6g3bmynzGyn61+fo90HjngdCQAAhElBCrn8Bi/m7Q/7TFI951wrSV9L+k++L2Q2yMzSzCwtIyPj5JKeRDgAKC4aJJXRK9elKn3nQQ16a54OZ2V7HQnwpeBhlgy1BOAHBSnk0iUF97DVlrQ5uIFzbodz7nDg8BVJbfN7IefcKOdcqnMuNSkp6VTyAgDyaJdSSU9c0Vpz1u/UfR8sVk4OY88BAIh2BSnk5kpqZGYpZpYgqb+kccENzKxG0OHFkpaHLuLxGaudAIAubl1T9/dqqnGLNuuZb1Z7HQcAABSxE24I7pzLMrOhkiZKipU02jm31MwelpTmnBsn6XYzu1hSlqSdkm4owswAgHwM7lpf63/er+e+Wa3mNcqqV8saJ74JQL6Ch1duGH6hh0kAIH8nLOQkyTk3QdKEPOceDHr+gKQHQhvtRJnC+W4AEPnMTA/3balV2/br7vcXKaVKGTWpXtbrWAAAoAgUaEPwSMbASgD4n8T4WI0c0FalS8Rp0FtprGQJAECUKlCPHADAP6qVS9TL17ZV/1GzdNu7C/TGje0UG8OvvYBTlXcVS4ZaAogEvu2Rc7/bAQEAcFTbuhX1cN+Wmr76Z704eY3XcQAAQIj5tpA7ikUrASB//c9MVt/Ta+rpr1dpzvqdXscBAAAh5PtCDgCQPzPTP/9wmupUKqXb312gnb8wXw4AgGjh20KOVSsB4MTKlIjTC1efoZ2/HNG9HyyS4y9PoNDqDRv/6wMAvOLbQu4ohlYCwPG1rFVef76gqb5dsV3vzNnodRwAABACvl21kt8pA0DBXd+xnr5evl3/HL9cXRolKblSKa8jAVGBjcMBeMX/PXLsJAcAJ2RmerRfK8WY6f+NXaycHH4dBgCAn/m+kAMAFEytCiX1tz7NNGvdDv139o9exwEAAIXg36GVTNgHgJN2RWqyJizZquFfrFD3ZtVUs0JJryMBUYONwwGEk/975BhZCQAFZmZ65JKWynFOj4xf5nUcAABwinzbIwcAODXJlUppaLeGeuKrVZq2KkNdGid5HQmISiyEAqAo+bZHjoGVAHDq/tSlvupXKa0HP/1BhzKzvY4T8cysl5mtNLM1ZjYsn+slzOy9wPXZZlYvcL6emR00s4WBx8vhzg4AiE6+75FjZCUAnLwScbF66OIWum70HL05a4MGdWngdaSIZWaxkkZI6iEpXdJcMxvnnAsemzpQ0i7nXEMz6y/pUUlXBq6tdc6dHtbQiDjMnwMQar7tkQMAFE6Xxknq0jhJIyav1Z6DmV7HiWTtJK1xzq1zzh2RNEZS3zxt+kr6T+D5WEnnmRm/awQAFBnfFnIsWgkAhXd/rybaczBTI6eu9TpKJKslaVPQcXrgXL5tnHNZkvZIqhy4lmJmC8xsqpmdXdRh4Q/1ho3/9QEAp8K3hdxR/MITAE5di5rl1ff0mhr93Xpt23vI6ziRKr8Pmry/TjxWmy2S6jjn2ki6W9I7ZlYu3zcxG2RmaWaWlpGRUajAAIDo5/tCDgBQOPf0aKLMbKdXpq3zOkqkSpeUHHRcW9LmY7UxszhJ5SXtdM4dds7tkCTn3DxJayU1zu9NnHOjnHOpzrnUpCRWEgUAHJ+PFzthbCUAhEKdyqV0ceuaemfORg3p1lAVSyd4HSnSzJXUyMxSJP0kqb+kq/O0GSfpekmzJPWT9K1zzplZknILumwzqy+pkSQqZvwGC6EAOBW+75FjYCUAFN7NXevrwJFsvTnrR6+jRJzAnLehkiZKWi7pfefcUjN72MwuDjR7TVJlM1uj3CGUR7co6CJpsZktUu4iKIOdczvD+xUAAKKRb3vkWOwEAEKnafVyOq9pVb0xc71u7lpfifGxXkeKKM65CZIm5Dn3YNDzQ5Iuz+e+DyV9WOQBEVWOtwAKvXUAjvJ/jxxdcgAQEgM7p2jXgUxNWLLF6ygAAOAEfF/IAQBC46wGlZVSpbTemb3R6ygAAOAE/Du00usAABBlzExXtUvWvyas0Kpt+9S4WlmvIwHII3jYJcMsgeLN9z1yxnInABAy/domKz7WNHZeutdRAADAcfi2Rw4AEHqVSifo7EZJGr94ix7o3VTGRGQgYrFtAVC8+bZHjlUrAaBoXHBaDf20+6AWpe/xOgoAADgG3/fI8ctiAAitHs2rKT7WNGHJFp2eXMHrOAAKiPlzQPHi2x45AEDRKF8yXu1TKmvqygyvowA4RfWGjf/1ASA6+bZHzjG2EgCKTKeGVfTolyu0fd8hVS2b6HUcAIXAXDogOvm+R46RlQAQep0bVpEkzVyzw+MkAAAgP77tkQMAFJ0WNcupXGKcZq/fqUva1PI6DoAQYi4dEB18W8gxsBIAik5MjKlV7Qpa8tNur6MAKEIMuwT8y7eF3K8YWwkAReK02uX16vR1OpSZrcT4WK/jAAgDeusA//BtIcdaJwBQtE5PrqCUKqWVse+wkiuV8joOgDCjtw6IbL4t5I4yuuQAoEj0bFFdPVtU9zoGAADIh+8LOQAAABQ9hl0CkcW3hZxjuRMAAABPMOwS8J7/95FjZCUAAACAYqZAPXJm1kvSs5JiJb3qnBue53oJSW9Kaitph6QrnXMbQhsVAAAAkShvD10weuuAonHCQs7MYiWNkNRDUrqkuWY2zjm3LKjZQEm7nHMNzay/pEclXVkUgX/FyEoAAICIx9w6oGgUpEeunaQ1zrl1kmRmYyT1lRRcyPWV9FDg+VhJL5iZOVf0mwQwshIAAMAfmFsHhE5BCrlakjYFHadLan+sNs65LDPbI6mypJ9DERIAAADR51hDMinwgBMrSCGXX6dX3p62grSRmQ2SNEiS6tSpU4C3PrZWyRX0+W2dVa9K6UK9DgAAACILPXfAiRWkkEuXlBx0XFvS5mO0STezOEnlJe3M+0LOuVGSRklSampqoYZdlikRp5a1yhfmJQAAAOADzLMDfq8ghdxcSY3MLEXST5L6S7o6T5txkq6XNEtSP0nfhmN+HAAAAIoXeuuAXCcs5AJz3oZKmqjc7QdGO+eWmtnDktKcc+MkvSbpLTNbo9yeuP5FGRoAAACQjr/1QTAKPkSbAu0j55ybIGlCnnMPBj0/JOny0EYDAAAAQoOePESbAhVyAAAAQDRhE3P4HYUcAAAAEIThmvADCjkAAADgFNCrBy9RyAEAcAJm1kvSs8pd9OtV59zwPNdLSHpTUltJOyRd6ZzbELj2gKSBkrIl3e6cmxjG6AA8Qq8eihqFHAAAx2FmsZJGSOqh3H1T55rZOOfcsqBmAyXtcs41NLP+kh6VdKWZNVfuSs4tJNWU9LWZNXbOZYf3qwAQqQpa8OVFAQgKOQAAjq+dpDXOuXWSZGZjJPWVFFzI9ZX0UOD5WEkvmJkFzo9xzh2WtD6wTU875e67CgCn7FQLwGAUg/5GIQcAwPHVkrQp6DhdUvtjtQnsv7pHUuXA+e/z3Fur6KICQMGFohgMFYrKk+dZITdv3ryfzezHQr5MFUk/hyJPGJE5PMgcHmQODz9mln6bu66XQQrJ8jnnCtimIPfmvoDZIEmDAof7zWxlgRP+nh9/ZsgcPn7MTebw8CyzPXrKt/r9+3zKn4+eFXLOuaTCvoaZpTnnUkORJ1zIHB5kDg8yh4cfM0v+zZ2PdEnJQce1JW0+Rpt0M4uTVF7SzgLeK0lyzo2SNCoUgf34vSdz+PgxN5nDg8zhEarMMaEIAwBAFJsrqZGZpZhZgnIXLxmXp804SdcHnveT9K1zzgXO9zezEmaWIqmRpDlhyg0AiGLMkQMA4DgCc96GSpqo3O0HRjvnlprZw5LSnHPjJL0m6a3AYiY7lVvsKdDufeUujJIlaQgrVgIAQsHvhVxIhqCEGZnDg8zhQebw8GNmyb+5f8c5N0HShDznHgx6fkjS5ce495+S/lmkAX/Pj997MoePH3OTOTzIHB6hGUafO/IDAAAAAOAXzJEDAAAAAJ+J+ELOzHqZ2UozW2Nmw/K5XsLM3gtcn21m9cKf8neZks1sspktN7OlZnZHPm3OMbM9ZrYw8Hgwv9cKJzPbYGZLAnnS8rluZvZc4Hu92MzO8CJnUJ4mQd+/hWa218zuzNPG8++zmY02s+1m9kPQuUpmNsnMVgf+rHiMe68PtFltZtfn1yaMmR83sxWB//Yfm1mFY9x73J+jMGd+yMx+Cvrvf8Ex7j3u3zNhzvxeUN4NZrbwGPd69X3O9++3SP+ZjlZ++4zk8zE8/PL5GMjBZ6R3mfmMDH3m8H5GOuci9qHcSeVrJdWXlCBpkaTmedrcKunlwPP+kt6LgNw1JJ0ReF5W0qp8cp8j6XOvs+bJtEFSleNcv0DSF8rdF6mDpNleZ87zs7JVUt1I+z5L6iLpDEk/BJ17TNKwwPNhkh7N575KktYF/qwYeF7Rw8znS4oLPH80v8wF+TkKc+aHJN1bgJ+d4/49E87Mea4/KenBCPs+5/v3W6T/TEfjw4+fkXw+evZzEpGfj4EcfEZ6l5nPyNBnDutnZKT3yLWTtMY5t845d0TSGEl987Tpq//f3p282FFFcRz/HkhAHIgT0RhdGHHlwoEgzptITETigEhEUIwgAbNw5yI7/wA3Ii4cMEoQce5FxARduGoVGxOViGlXNmk6oJIoboweF/eWFtVVr0vtvkPz+8DjDXUfHE7ddw+33q0q2BdfvwVsMbO+G7Am4+7z7j4TX/8CHAU25oxpmdwNvOrBNHCumW3IHVS0Bfje3f/vTeaXnbt/QriKXVu73+4D7un56h3AIXf/yd1/Bg4B21Ys0Ja+mN39oLufjm+nCffDKsZAnscYM86siEkxx3HsAeD1FLGMNWF8K7pPr1LV1UjVxyyKrY+gGpmKamQaqWtk6RO5jcAPrfdzLB7w/24Tf0AngQuSRDdCXMZyLfBpz+YbzeywmX1gZlclDayfAwfN7Asze7xn+5j9kctOhn/MpeUZ4CJ3n4fwowfW97QpOd+7CEef+yzVj1LbE5e6vDywlKHUPN8KLLj7sYHt2fPcGd9q79M1qrpGqj4mU1t9hPrHE9XIlacaSfkTub6jht3LbI5pk4WZnQ28DTzp7qc6m2cIyxyuBp4F3ksdX4+b3f06YDvwhJnd1tleZK4t3KB3B/Bmz+YS8zxWqfneS7gf1v6BJkv1o5SeB64ArgHmCcswuorMM/Agk480Zs3zEuPb4Nd6Pish17WqtkaqPqaxiusjlJtz1cg0VCMpfyI3B1zWen8pcHyojZmtAdbx3/46XlZmtpawA/e7+zvd7e5+yt1/ja8PAGvN7MLEYXZjOh6fTwDvEv5ObxuzP3LYDsy4+0J3Q4l5jhaaZTfx+URPm+LyHU+8vQt4yOOC7q4R/SgZd19w9z/c/U/ghYFYSszzGuA+4I2hNjnzPDC+VdmnK1dljVR9TKrG+giVjieqkWmoRv6j9Inc58CVZnZ5PKq0E5jqtJkCmqu63A98PPTjSSWu230JOOruzwy0ubg5T8HMrifsix/TRbkonrPM7JzmNeGk3a87zaaAhy24ATjZ/E2c2eBRmdLy3NLut48A7/e0+RDYambnxeUOW+NnWZjZNuApYIe7/zbQZkw/SqZzjsq9A7GMGWdSux341t3n+jbmzPOE8a26Pr0KVFcjYnE4RwAAAUBJREFUVR+Tq7E+QoXjiWpkUqqRDU98NZd/+yBcCeo7whVz9sbPnib8UADOICwZmAU+AzYVEPMthL9CjwBfxsedwG5gd2yzB/iGcPWfaeCmzDFvirEcjnE1uW7HbMBzcV98BWwuINdnEgrPutZnReWZUETngd8JR1seI5yj8hFwLD6fH9tuBl5sfXdX7NuzwKOZY54lrN1u+nRzJbxLgAOT+lHGmF+LffUIYRDd0I05vl80zuSKOX7+StOHW21LyfPQ+FZ0n16tj76+S8E1ckL/KWrc7sSs+riycapG5otZNXL5Y05aIy1+SURERERERCpR+tJKERERERER6dBETkREREREpDKayImIiIiIiFRGEzkREREREZHKaCInIiIiIiJSGU3kREREREREKqOJnIiIiIiISGU0kRMREREREanMX8HnS2UTFtTzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "val=np.sqrt(lambda_function_np_2(mB_mass, mKst_mass, q2_range))*beta_l_np(q2_range, mmu_mass)\n",
    "plt.subplot(1,2,1)\n",
    "plt.plot(q2_range/1e6, val);\n",
    "fig=plt.gcf()\n",
    "plt.subplot(1,2,2)\n",
    "plt.hist(q2_array/1e6, weights=weights_np, bins=100,density=True);\n",
    "fig=plt.gcf()\n",
    "fig.set_size_inches(15,5);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#def matrix_elt(q2, m, m1, m2, cos_theta_l):\n",
    "#    \n",
    "#    out = a(q2, m, m1, m2)+b(q2, m, m1, m2)*cos_theta_l+c(q2, m, m1, m2)*ztf.square(cos_theta_l)\n",
    "#    \n",
    "#    return out \n",
    "#\n",
    "#def a(q2, m, m1, m2):\n",
    "#    \n",
    "#    step1=q2*(1+ztf.square(beta_l(q2, m)))+ lambda_function(m1, m2, q2)/2 + 2*m*(ztf.square(m1)-ztf.square(m2)+q2)+4*ztf.square(m)*ztf.square(m1)\n",
    "#    #step2 = beta_l(q2, m)*ztf.sqrt(lambda_phase(m1, m2, q2))*step1\n",
    "#    \n",
    "#    return step1\n",
    "#\n",
    "#def b(q2, m, m1, m2):\n",
    "#    \n",
    "#    step1=2*(q2*(1+ztf.square(beta_l(q2, m)))+ m*(ztf.sqrt(lambda_function(m1, m2, q2))*beta_l(q2,m)+(ztf.square(m1)-ztf.square(m2)+q2)))\n",
    "#    #step2 = beta_l(q2, m)*ztf.sqrt(lambda_phase(m1, m2, q2))*step1\n",
    "#    \n",
    "#    return step1\n",
    "#\n",
    "#def c(q2, m, m1, m2):\n",
    "#    \n",
    "#    step1= q2*(1+ztf.square(beta_l(q2, m)))- lambda_function(m1, m2, q2)/2 + 2*m*ztf.sqrt(lambda_function(m1, m2, q2))*beta_l(q2, m)\n",
    "#    #step2 = beta_l(q2, m)*ztf.sqrt(lambda_phase(m1, m2, q2))*step1\n",
    "#    \n",
    "#    return step1\n",
    "#\n",
    "#def phase_space(q2, m, m1, m2):\n",
    "#    \n",
    "#    phase_space_term = beta_l(q2, m)*ztf.sqrt(lambda_function(m1, m2, q2))\n",
    "#    \n",
    "#    return phase_space_term\n",
    "#    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "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": 9,
   "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, p1+p2)\n",
    "        \n",
    "        costheta_l= tf.expand_dims(get_costheta_l(p1q, z), axis=1)\n",
    "        \n",
    "        pdf = matrix_elt(q2, ml, mB, mKst, costheta_l)\n",
    "        \n",
    "        return pdf[:, 0]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "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": 11,
   "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": 12,
   "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": 13,
   "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": 14,
   "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"
     ]
    }
   ],
   "source": [
    "sampler = pdf.create_sampler(n=n_events)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "zfit.settings.set_verbosity(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "for i in range(1):\n",
    "    sampler.resample()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3kAAAEvCAYAAAD4uAgWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3df7DddX3n8edLFGtrW4IGNs2Phta0FZ0psneBXXc6LLQQ0GnsjHSwnZpSZtLtwK6uzpbQdhaL0oHdVlpXSyctqaFjjRR1yLgoTRHGcUaQoBSFlOUWWHJNFuIGUJctFn3vH+cTPUnOvSH3nvvjfM/zMXPmnO/7+/me+/l+c5NP3t/Pj2+qCkmSJElSN7xksSsgSZIkSRoekzxJkiRJ6hCTPEmSJEnqEJM8SZIkSeoQkzxJkiRJ6hCTPEmSJEnqkJcudgVm69WvfnWtXbt2sashSZpn991339eravli12NU2D5K0viYro0c2SRv7dq17Nq1a7GrIUmaZ0n+12LXYZTYPkrS+JiujXS4piRJkiR1iEmeJEnzJMlxSb6c5FNt+5Qk9yR5JMnHkhzf4i9v25Nt/9q+77iyxR9Ocv7inIkkaZSY5EmSNH/eAezu274OuL6q1gFPA5e2+KXA01X1GuD6Vo4kpwIXA68D1gN/muS4Baq7JGlEmeRJkjQPkqwC3gT8RdsOcA5wSyuyDXhL+7yhbdP2n9vKbwC2V9XzVfUYMAmcsTBnIEkaVSZ5kiTNjz8Gfhv4btt+FfBMVb3QtqeAle3zSmAPQNv/bCv/vfiAYyRJGsgkT5KkIUvyZuCpqrqvPzygaB1l30zH9P+8TUl2Jdm1f//+Y66vJKlbTPIkSRq+NwK/mORxYDu9YZp/DJyQ5ODji1YBe9vnKWA1QNv/o8CB/viAY76nqrZU1URVTSxf7iMFJWncmeRJkjRkVXVlVa2qqrX0Fk75bFX9KnAn8NZWbCNwa/u8o23T9n+2qqrFL26rb54CrAO+uECnIUkaUSP7MHRJkkbQFcD2JO8Dvgzc2OI3An+VZJJeD97FAFX1YJKbgYeAF4DLquo7C19tSdIoMcmTJGkeVdVdwF3t86MMWB2zqv4JuGia468Brpm/GkqSusbhmpIkSZLUIfbkSUvQ2s3/45jKP37tm+apJpIkLS2D2kjbQelQ9uRJkiRJUoeY5EmSJElSh5jkSZIkSVKHHHVOXpKtwJuBp6rq9S12IvAxYC3wOPDLVfV0kgB/AlwIPAf8elV9qR2zEfi99rXvq6ptLf4vgQ8DrwBuA97Rng0kSZIkHdV0c9mdq6dx9WIWXvkw8EHgpr7YZuCOqro2yea2fQVwAb0Hta4DzgRuAM5sSeFVwARQwH1JdlTV063MJuBuekneeuDTcz81aek71gVWjuV7bNgkSZLG01GHa1bV5+g9mLXfBmBb+7wNeEtf/KbquRs4IckK4HxgZ1UdaIndTmB92/cjVfWF1nt3U993SZIkSZKO0Wzn5J1cVfsA2vtJLb4S2NNXbqrFZopPDYhLkiRJkmZh2AuvZECsZhEf/OXJpiS7kuzav3//LKsoSZIkSd012yTvyTbUkvb+VItPAav7yq0C9h4lvmpAfKCq2lJVE1U1sXz58llWXZIkSZK668UsvDLIDmAjcG17v7UvfnmS7fQWXnm2qvYluR34gyTLWrnzgCur6kCSbyY5C7gHeDvw32dZJ2lJG9YiK5IkSdJMXswjFD4KnA28OskUvVUyrwVuTnIp8ARwUSt+G73HJ0zSe4TCJQAtmXsvcG8rd3VVHVzM5bf4/iMUPo0ra0qSJI09b45Ks3fUJK+q3jbNrnMHlC3gsmm+ZyuwdUB8F/D6o9VDkiRJknR0sx2uKWmJ88GwkqRxZ1uocTXs1TUlSZIkSYvIJE+SJEmSOsQkT5IkSZI6xCRPkiRJkjrEJE+SJEmSOsQkT5IkSZI6xCRPkiRJkjrEJE+SJEmSOsSHoUtDNt2DV5cKHwwrSZLUbfbkSZI0ZEl+IMkXk/x9kgeT/H6LfzjJY0nub6/TWjxJPpBkMskDSU7v+66NSR5pr42LdU6SpNFhT54kScP3PHBOVX0rycuAzyf5dNv3n6vqlsPKXwCsa68zgRuAM5OcCFwFTAAF3JdkR1U9vSBnIUkaSfbkSZI0ZNXzrbb5svaqGQ7ZANzUjrsbOCHJCuB8YGdVHWiJ3U5g/XzWXZI0+uzJkyRpHiQ5DrgPeA3woaq6J8lvAdck+S/AHcDmqnoeWAns6Tt8qsWmi0udsdTnskujyJ48SZLmQVV9p6pOA1YBZyR5PXAl8DPAvwJOBK5oxTPoK2aIHyLJpiS7kuzav3//UOovSRpdJnmSJM2jqnoGuAtYX1X72pDM54G/BM5oxaaA1X2HrQL2zhA//GdsqaqJqppYvnz5PJyFJGmUOFxTkqQhS7Ic+OeqeibJK4CfB65LsqKq9iUJ8Bbgq+2QHcDlSbbTW3jl2VbuduAPkixr5c6j1xsoaQ58nJC6ziRPkqThWwFsa/PyXgLcXFWfSvLZlgAGuB/49638bcCFwCTwHHAJQFUdSPJe4N5W7uqqOrCA5yFJGkEmedIsOVFc0nSq6gHgDQPi50xTvoDLptm3Fdg61ApKkjrNOXmSJEmS1CEmeZIkSZLUIQ7XlAQ4CV2SJKkr7MmTJEmSpA4xyZMkSZKkDnG4piRJkoZqVFegHlRvpy1oFNmTJ0mSJEkdYpInSZIkSR3icE1JM3LoiiRJ0mixJ0+SJEmSOsSePOkoRnXyuCRJksaTPXmSJEmS1CEmeZIkSZLUISZ5kiRJktQhJnmSJEmS1CEmeZIkSZLUISZ5kiRJktQhJnmSJEmS1CEmeZIkSZLUIT4MXWp86LkkSTrcdP8/ePzaNy1wTaQXzyRP0jGzwZMkSVq6HK4pSZIkSR0ypyQvyX9K8mCSryb5aJIfSHJKknuSPJLkY0mOb2Vf3rYn2/61fd9zZYs/nOT8uZ2SJEmSJI2vWSd5SVYC/xGYqKrXA8cBFwPXAddX1TrgaeDSdsilwNNV9Rrg+laOJKe2414HrAf+NMlxs62XJEmSJI2zuQ7XfCnwiiQvBX4Q2AecA9zS9m8D3tI+b2jbtP3nJkmLb6+q56vqMWASOGOO9ZIkadG0kS1fTPL3bcTL77e4o10kSfNu1kleVX0N+EPgCXrJ3bPAfcAzVfVCKzYFrGyfVwJ72rEvtPKv6o8POEaSpFH0PHBOVf0scBqwPslZONpFkrQAZr26ZpJl9HrhTgGeAf4GuGBA0Tp4yDT7posP+pmbgE0Aa9asOcYaS5K0MKqqgG+1zZe1V9Eb7fIrLb4NeA9wA7329D0tfgvwwcNHuwCPJTk42uUL838W0tH5+CFpaZrLIxR+HnisqvYDJPkE8G+AE5K8tPXWrQL2tvJTwGpgqg3v/FHgQF/8oP5jDlFVW4AtABMTEwMTQUmLx0crSN/XetzuA14DfAj4R17kaJck/aNd7u77Wke7SJKOai5J3hPAWUl+EPh/wLnALuBO4K3AdmAjcGsrv6Ntf6Ht/2xVVZIdwF8neT/wY8A64ItzqJd0VN55lDTfquo7wGlJTgA+Cbx2ULH2PqfRLo50kST1m8ucvHvoDSn5EvCV9l1bgCuAd7UhJa8CbmyH3Ai8qsXfBWxu3/MgcDPwEPAZ4LLWMEqSNPKq6hngLuAs2miXtmvQaBdmM9qlqrZU1URVTSxfvnw+TkOSNELm0pNHVV0FXHVY+FEGrI5ZVf8EXDTN91wDXDOXukiStFQkWQ78c1U9k+QV9KY4XIejXSRJC2BOSZ4kSRpoBbCtzct7CXBzVX0qyUPA9iTvA77MoaNd/qqNdjlAb0VNqurBJAdHu7yAo10kSS+CSZ4kSUNWVQ8AbxgQd7SLJGnemeRJkiRJx8gVpbWUzXrhFUmSJEnS0mNPnjrNRyVIkiRp3NiTJ0mSJEkdYpInSZIkSR3icE1J887J6ZIkSQvHnjxJkiRJ6hB78iRJkjQjFzKTRos9eZIkSZLUISZ5kiRJktQhDteUJEmShsTFxrQU2JMnSZIkSR1ikidJkiRJHeJwTUmSJAGuoil1hUmepEUz6D8TzlmQJEmaG5M8dYJ3HiVJkqQe5+RJkiRJUoeY5EmSJElSh5jkSZIkSVKHmORJkiRJUoeY5EmSJElSh5jkSZIkSVKHmORJkiRJUof4nDyNFJ+H133T/Rn7kHRJ0igb1L7Ztmm+2JMnSZIkSR1ikidJ0pAlWZ3kziS7kzyY5B0t/p4kX0tyf3td2HfMlUkmkzyc5Py++PoWm0yyeTHOR5I0WhyuKUnS8L0AvLuqvpTkh4H7kuxs+66vqj/sL5zkVOBi4HXAjwF/l+Sn2u4PAb8ATAH3JtlRVQ8tyFmos5z+IHWbSZ4kSUNWVfuAfe3zN5PsBlbOcMgGYHtVPQ88lmQSOKPtm6yqRwGSbG9lTfIkSdNyuKYkSfMoyVrgDcA9LXR5kgeSbE2yrMVWAnv6Dptqsenih/+MTUl2Jdm1f//+IZ+BJGnUmORJkjRPkrwS+Djwzqr6BnAD8JPAafR6+v7oYNEBh9cM8UMDVVuqaqKqJpYvXz6UukuSRpfDNSVJmgdJXkYvwftIVX0CoKqe7Nv/58Cn2uYUsLrv8FXA3vZ5urgkSQOZ5EmSNGRJAtwI7K6q9/fFV7T5egC/BHy1fd4B/HWS99NbeGUd8EV6PXnrkpwCfI3e4iy/sjBnIWm++WxYzReTPEkjwYZQI+aNwK8BX0lyf4v9DvC2JKfRG3L5OPCbAFX1YJKb6S2o8gJwWVV9ByDJ5cDtwHHA1qp6cCFPRJI0ekzytGS5vLOkUVVVn2fwfLrbZjjmGuCaAfHbZjpOkqTDufCKJEmSJHWISZ4kSZIkdYhJniRJkiR1iHPyJEmSpCXExcY0V/bkSZIkSVKHzCnJS3JCkluS/EOS3Un+dZITk+xM8kh7X9bKJskHkkwmeSDJ6X3fs7GVfyTJxrmelCRJkiSNq7n25P0J8Jmq+hngZ4HdwGbgjqpaB9zRtgEuoPdw13XAJuAGgCQnAlcBZwJnAFcdTAwlSZIkScdm1nPykvwI8HPArwNU1beBbyfZAJzdim0D7gKuADYAN1VVAXe3XsAVrezOqjrQvncnsB746GzrJml8OG9BkiTpUHPpyfsJYD/wl0m+nOQvkvwQcHJV7QNo7ye18iuBPX3HT7XYdHFJkiRJ0jGaS5L3UuB04IaqegPwf/n+0MxBMiBWM8SP/IJkU5JdSXbt37//WOsrSZIkSZ03l0coTAFTVXVP276FXpL3ZJIVVbWvDcd8qq/86r7jVwF7W/zsw+J3DfqBVbUF2AIwMTExMBGUJEnS9003rF1Sd806yauq/51kT5KfrqqHgXOBh9prI3Bte7+1HbIDuDzJdnqLrDzbEsHbgT/oW2zlPODK2dZLo8fGR5IkSRqeuT4M/T8AH0lyPPAocAm9IaA3J7kUeAK4qJW9DbgQmASea2WpqgNJ3gvc28pdfXARFkmSJEnSsZlTkldV9wMTA3adO6BsAZdN8z1bga1zqYskSZIkae7PyZMkSZIkLSFzHa4pSZKkJcA57pIOMsnTgrHxkSRJkuafwzUlSZIkqUNM8iRJkiSpQxyuKUmSJI2A6aa+PH7tmxa4JlrqTPIkdZINoSRJGlcO15QkSZKkDjHJkyRJkqQOMcmTJEmSpA4xyZMkaciSrE5yZ5LdSR5M8o4WPzHJziSPtPdlLZ4kH0gymeSBJKf3fdfGVv6RJBsX65wkSaPDJE+SpOF7AXh3Vb0WOAu4LMmpwGbgjqpaB9zRtgEuANa11ybgBuglhcBVwJnAGcBVBxNDSZKm4+qakiQNWVXtA/a1z99MshtYCWwAzm7FtgF3AVe0+E1VVcDdSU5IsqKV3VlVBwCS7ATWAx9dsJPRkjPd6sGSdJBJnqSxMug/Rz5WQfMpyVrgDcA9wMktAaSq9iU5qRVbCezpO2yqxaaLH/4zNtHrAWTNmjXDPQFJ0sgxydPQeYdRknqSvBL4OPDOqvpGkmmLDojVDPFDA1VbgC0AExMTR+yXJI0X5+RJkjQPkryMXoL3kar6RAs/2YZh0t6favEpYHXf4auAvTPEJUmalkmeJElDll6X3Y3A7qp6f9+uHcDBFTI3Arf2xd/eVtk8C3i2Deu8HTgvybK24Mp5LSZJ0rQcrilJ0vC9Efg14CtJ7m+x3wGuBW5OcinwBHBR23cbcCEwCTwHXAJQVQeSvBe4t5W7+uAiLJIkTcckT5KkIauqzzN4Ph3AuQPKF3DZNN+1Fdg6vNpJkrrOJE+SJEkaYdMteufq0ePLOXmSJEmS1CEmeZIkSZLUISZ5kiRJktQhJnmSJEmS1CEuvCJp7DlhXdJSNN2/TZJ0NCZ5mhMbIEmSJGlpcbimJEmSJHWISZ4kSZIkdYhJniRJkiR1iEmeJEmSJHWISZ4kSZIkdYira0qSJEkdNGgVdB8PNB5M8iRpGj4/T5IkjSKHa0qSJElSh9iTpxfFh55LkjR/bGclDZM9eZIkSZLUIfbkSZIkSWPC+ebjwZ48SZIkSeoQkzxJkiRJ6hCTPEmSJEnqEJM8SZIkSeoQF17RIVzCWTo6J61LkqSlbM49eUmOS/LlJJ9q26ckuSfJI0k+luT4Fn95255s+9f2fceVLf5wkvPnWidJkiRJGlfDGK75DmB33/Z1wPVVtQ54Gri0xS8Fnq6q1wDXt3IkORW4GHgdsB740yTHDaFekiQtiiRbkzyV5Kt9sfck+VqS+9vrwr59A292JlnfYpNJNi/0eUiSRtOchmsmWQW8CbgGeFeSAOcAv9KKbAPeA9wAbGifAW4BPtjKbwC2V9XzwGNJJoEzgC/MpW6SJC2iDwMfBG46LH59Vf1hf+Cwm50/Bvxdkp9quz8E/AIwBdybZEdVPTSfFdf8clqEpIUw1568PwZ+G/hu234V8ExVvdC2p4CV7fNKYA9A2/9sK/+9+IBjDpFkU5JdSXbt379/jlWXJGl+VNXngAMvsvj3bnZW1WPAwZudZwCTVfVoVX0b2N7KSpI0o1n35CV5M/BUVd2X5OyD4QFF6yj7Zjrm0GDVFmALwMTExMAykiQtYZcneTuwC3h3VT1N78bm3X1l+m92Hn4T9MwFqaWkseOiYt0yl568NwK/mORxencXz6HXs3dCkoPJ4ypgb/s8BawGaPt/lN5dzu/FBxwjSVJX3AD8JHAasA/4oxaf801QR7pIkvrNOsmrqiuralVVraU3l+CzVfWrwJ3AW1uxjcCt7fOOtk3b/9mqqha/uK2+eQqwDvjibOslSdJSVFVPVtV3quq7wJ/TG44J09/sfNE3QatqS1VNVNXE8uXLh195SdJImY+HoV9BbxGWSXpz7m5s8RuBV7X4u4DNAFX1IHAz8BDwGeCyqvrOPNRLkqRFk2RF3+YvAQdX3pzuZue9wLr2aKLj6d1Q3bGQdZYkjaahPAy9qu4C7mqfH+X7dyf7y/wTcNE0x19Db4VOLRBX95Kk+ZPko8DZwKuTTAFXAWcnOY3ekMvHgd+E3s3OJAdvdr5A383OJJcDtwPHAVvbjVFJkmY0lCRPkjT45okT1sdTVb1tQPjGAbGD5Qfe7Kyq24Dbhlg1SdIYmI/hmpIkSZKkRWJPniRJ0pA5LUJd4aMVRpM9eZIkSZLUISZ5kiRJktQhJnmSJEmS1CHOyZMkSZol595JWorsyZMkSZKkDjHJkyRJkqQOcbimJM0jl56WJEkLzSRvDDhfQJIkSRofJnmSJEmSjsmgTgRHqSwdzsmTJEmSpA4xyZMkSZKkDjHJkyRJkqQOMcmTJEmSpA4xyZMkSZKkDjHJkyRJkqQO8REKkiRJkuZsumcz+2iFhWeS1yE+9FySJEmSSZ4kLQLvdkqjx5upkkaFc/IkSZIkqUNM8iRJkiSpQ0zyJEmSJKlDTPIkSZIkqUNceEWSpCFLshV4M/BUVb2+xU4EPgasBR4Hfrmqnk4S4E+AC4HngF+vqi+1YzYCv9e+9n1VtW0hz2NcucCKpFFnkidJ0vB9GPggcFNfbDNwR1Vdm2Rz274CuABY115nAjcAZ7ak8CpgAijgviQ7qurpBTsLSRoCV5ReeA7XlCRpyKrqc8CBw8IbgIM9cduAt/TFb6qeu4ETkqwAzgd2VtWBltjtBNbPf+0lSaPOJE+SpIVxclXtA2jvJ7X4SmBPX7mpFpsuLknSjByuOYKcKyB1l0NaxlIGxGqG+JFfkGwCNgGsWbNmeDWTJI0ke/IkSVoYT7ZhmLT3p1p8CljdV24VsHeG+BGqaktVTVTVxPLly4decUnSaDHJkyRpYewANrbPG4Fb++JvT89ZwLNtOOftwHlJliVZBpzXYpIkzcjhmpIkDVmSjwJnA69OMkVvlcxrgZuTXAo8AVzUit9G7/EJk/QeoXAJQFUdSPJe4N5W7uqqOnwxF82B0x8kdZVJniRJQ1ZVb5tm17kDyhZw2TTfsxXYOsSqSdKS4Tz0+eNwTUmSJEnqEJM8SZIkSeoQh2suYc4VkCRJknSs7MmTJEmSpA6xJ0+SRsCgnn0npksvjiNjJI0be/IkSZIkqUNM8iRJkiSpQ2ad5CVZneTOJLuTPJjkHS1+YpKdSR5p78taPEk+kGQyyQNJTu/7ro2t/CNJNs79tCRJkiRpPM2lJ+8F4N1V9VrgLOCyJKcCm4E7qmodcEfbBrgAWNdem4AboJcUAlcBZwJnAFcdTAwlSZIkScdm1kleVe2rqi+1z98EdgMrgQ3AtlZsG/CW9nkDcFP13A2ckGQFcD6ws6oOVNXTwE5g/WzrJUmSJEnjbCirayZZC7wBuAc4uar2QS8RTHJSK7YS2NN32FSLTReXJEk6Jq6kKUlDSPKSvBL4OPDOqvpGkmmLDojVDPFBP2sTvaGerFmz5tgrK0mSJGlJ87FBczenJC/Jy+gleB+pqk+08JNJVrRevBXAUy0+BazuO3wVsLfFzz4sftegn1dVW4AtABMTEwMTwVHkXUdJszHdvx02hJIkjbe5rK4Z4EZgd1W9v2/XDuDgCpkbgVv74m9vq2yeBTzbhnXeDpyXZFlbcOW8FpMkSZIkHaO59OS9Efg14CtJ7m+x3wGuBW5OcinwBHBR23cbcCEwCTwHXAJQVQeSvBe4t5W7uqoOzKFekiRJkjS2Zp3kVdXnGTyfDuDcAeULuGya79oKbJ1tXSRJkiR1l1MUjs1cnpMnSZIkSVpiTPIkSZIkqUNM8iRJkiSpQ0zyJEmSJKlD5vwwdEmSpIXmM2YlgQuyTMckb4HZKEmabzZ4kiSNN4drSpIkSVKHmORJkiRJUoeY5EmStICSPJ7kK0nuT7KrxU5MsjPJI+19WYsnyQeSTCZ5IMnpi1t7SdIocE6eJEkL799V1df7tjcDd1TVtUk2t+0rgAuAde11JnBDex8bzmWXpGNnkidJ0uLbAJzdPm8D7qKX5G0AbqqqAu5OckKSFVW1b1FqKUkjYtwXIXO4piRJC6uAv01yX5JNLXbywcStvZ/U4iuBPX3HTrWYJEnTsidPkqSF9caq2pvkJGBnkn+YoWwGxOqIQr1kcRPAmjVrhlNLSdLIsidPkqQFVFV72/tTwCeBM4Ank6wAaO9PteJTwOq+w1cBewd855aqmqiqieXLl89n9SVJI8CevHniRHFJS824z09YCpL8EPCSqvpm+3wecDWwA9gIXNveb22H7AAuT7Kd3oIrzzofT5J0NCZ5kiQtnJOBTyaBXhv811X1mST3AjcnuRR4Ariolb8NuBCYBJ4DLln4Ki8Mb45K0vCY5EmStECq6lHgZwfE/w9w7oB4AZctQNUkaSwMuqHUxREtzsmTJEmSpA4xyZMkSZKkDjHJkyRJkqQOMcmTJEmSpA5x4RVJkrRgXEVTkuafSd4c2VhJGnXjstKYJEmDdPE5sg7XlCRJkqQOMcmTJEmSpA4xyZMkSZKkDjHJkyRJkqQOceEVSZI0L1ycTJIWhz15kiRJktQh9uRJkiRJ0mFG+dEKJnkvkkNOJI2TUW7YJEkadw7XlCRJkqQOsSdPkiTNiaNdJI2TURjtYk+eJEmSJHWISZ4kSZIkdYhJniRJkiR1iEmeJEmSJHWIC68M4ARySRpsFCaba/7YPkrSaDDJkyRJkqQ5Wko3Qh2uKUmSJEkdYk+eJEk6hMMyJWm0LZmevCTrkzycZDLJ5sWujyRJS4HtoyTpWC2JnrwkxwEfAn4BmALuTbKjqh5a3JpJkrR4bB8lafQNGh0x3/P0lkSSB5wBTFbVowBJtgMbgHltxByOIknDsZQmm3fMorSPkqTRtlSGa64E9vRtT7WYJEnjzPZRknTMlkpPXgbE6ohCySZgU9v8VpKH57VWC+/VwNcXuxJLiNfjSF6TQ3k9DrXkrkeuG8rX/PhQvmU02T4uwd/rJcBrciSvyaG8HkdaUtdkSO0jTNNGLpUkbwpY3be9Cth7eKGq2gJsWahKLbQku6pqYrHrsVR4PY7kNTmU1+NQXo9OGvv20d/rI3lNjuQ1OZTX40jjdk2WynDNe4F1SU5JcjxwMbBjkeskSdJis32UJB2zJdGTV1UvJLkcuB04DthaVQ8ucrUkSVpUto+SpNlYEkkeQFXdBty22PVYZJ0cajMHXo8jeU0O5fU4lNejg2wf/b0ewGtyJK/JobweRxqra5KqI+ZvS5IkSZJG1FKZkydJkiRJGgKTvAWU5L8l+YckDyT5ZJIT+vZdmWQyycNJzu+Lr2+xySSb++KnJLknySNJPtYm5I+UJBcleTDJd5NMHLZv7K7HTKY77y5KsjXJU0m+2hc7McnO9ue7M8myFk+SD7Tr8kCS0/uO2djKP5Jk42Kcy1wlWZ3kziS729+Vd7T4WF4PdZNt46FsG4+N7eP4tge2kUdRVb4W6AWcB7y0fb4OuK59PhX4e+DlwCnAP9KbYH9c+/wTwPGtzKntmJuBi9vnPwN+a7HPbxbX47XATwN3ARN98bG8HjNcp2nPu4sv4OeA04Gv9sX+K7C5fd7c93fnQuDT9J4ldhZwT4ufCDza3pe1z8sW+9xmcS1WAKe3zz8M/M/292Msr4evbr5sG4+4HraNL/5a2Sx8QnAAAAMySURBVD6OcXtgGznzy568BVRVf1tVL7TNu+k97whgA7C9qp6vqseASeCM9pqsqker6tvAdmBDkgDnALe047cBb1mo8xiWqtpdVYMe2DuW12MGA897kes0b6rqc8CBw8Ib6P25wqF/vhuAm6rnbuCEJCuA84GdVXWgqp4GdgLr57/2w1VV+6rqS+3zN4HdwErG9Hqom2wbD2XbeExsH8e4PbCNnJlJ3uL5DXp3E6D3C7mnb99Ui00XfxXwTF+jeDDeFV6PQ0133uPk5KraB71/1IGTWvxYf1dGVpK1wBuAe/B6qLtsG6fn9TiS/7bZHgC2kYMsmUcodEWSvwP+xYBdv1tVt7Yyvwu8AHzk4GEDyheDk/CaofyS82Kux6DDBsQ6cT1mqevnNxfTXZtOXbMkrwQ+Dryzqr7Ru0E/uOiAWOeuh0aPbeOhbBuHZhzOcbbGpj2wjRzMJG/IqurnZ9rfJnO+GTi3qg7+Ak0Bq/uKrQL2ts+D4l+n18X80naHrr/8knK06zGNzl6PWZrpeoyLJ5OsqKp9bWjFUy0+3bWZAs4+LH7XAtRz6JK8jF7j9ZGq+kQLj+310GiybTyUbePQ2D6OeXtgGzk9h2suoCTrgSuAX6yq5/p27QAuTvLyJKcA64AvAvcC69rqWMcDFwM7WgN4J/DWdvxGYLo7f6PI63Gogee9yHVaaDvo/bnCoX++O4C3txWzzgKebUMzbgfOS7Ksrap1XouNlDan5kZgd1W9v2/XWF4PdZNt44vm9TiS7eMYtwe2kUdxrCu1+Jr9i94k6T3A/e31Z337fpfeClEPAxf0xS+kt1rQP9IbxnEw/hP0/nGfBP4GePlin98srscv0bt78jzwJHD7OF+Po1yrgefdxRfwUWAf8M/t9+NSenNL7gAeae8ntrIBPtSuy1c4dCW632i/D5PAJYt9XrO8Fv+W3pCRB/r+3bhwXK+Hr26+bBuPuB62jcd2vWwfx7Q9sI2c+ZV2YpIkSZKkDnC4piRJkiR1iEmeJEmSJHWISZ4kSZIkdYhJniRJkiR1iEmeJEmSJHWISZ4kSZIkdYhJniRJkiR1iEmeJEmSJHXI/weesXbEVJjhYgAAAABJRU5ErkJggg==\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": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "p1x=a['p1x'].values[:]\n",
    "p1y=a['p1y'].values[:]\n",
    "p1z=a['p1z'].values[:]\n",
    "p1E=a['p1E'].values[:]\n",
    "\n",
    "p2x=a['p2x'].values[:]\n",
    "p2y=a['p2y'].values[:]\n",
    "p2z=a['p2z'].values[:]\n",
    "p2E=a['p2E'].values[:]\n",
    "\n",
    "qx =a['qx'].values[:]\n",
    "qy =a['qy'].values[:]\n",
    "qz =a['qz'].values[:]\n",
    "qE =a['qE'].values[:]\n",
    "\n",
    "pKstx =a['pKstx'].values[:]\n",
    "pKsty =a['pKsty'].values[:]\n",
    "pKstz =a['pKstz'].values[:]\n",
    "pKstE =a['pKstE'].values[:]\n",
    "\n",
    "p1=np.array([p1x,p1y,p1z,p1E]).transpose()\n",
    "p2=np.array([p2x,p2y,p2z,p2E]).transpose()\n",
    "q = np.array([qx,qy,qz,qE]).transpose()\n",
    "pKst = np.array([pKstx,pKsty,pKstz,pKstE]).transpose()\n",
    "#q = p1 + p2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "q2_from_sample=scalar_product_4_np(q,q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3kAAAEvCAYAAAD4uAgWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAaiElEQVR4nO3df7Bmd10f8PfHXRMrVn6ErdUk66Yl7cxiKZUl6FgpU2rcGMtqm8gGpwabzko1bZ1Wa2g1YtSZxV9IS9qaktQAYsC0tDuTxchIx844QHcTEVxidE1XcgkDgaShkWJc+PSP5wl9eHju3nNzf+zdc1+vmcye8z3f89zv3W/O8zzv/X7P+VZ3BwAAgHH4krPdAAAAANaPkAcAADAiQh4AAMCICHkAAAAjIuQBAACMiJAHAAAwIjvPdgPmPfvZz+49e/ac7WYAAACcFffcc88nunvXUz1/y4W8PXv25Pjx42e7GQAAAGdFVf3xWs43XRMAAGBEhDwAAIAREfIAAABGRMgDAAAYESEPAABgRIQ8AACAERHyAAAARkTIAwAAGBEhDwAAYESEPAAAgBER8gAAAEZk59luwBjtueGuwXVPHb5yA1sCAABsN0byAAAARkTIAwAAGBEhDwAAYESEPAAAgBHx4JWzzENaAACA9STknUMEQgAAYCWmawIAAIyIkAcAADAiQh4AAMCICHkAAAAjIuQBAACMiJAHAAAwIkIeAADAiAh5AAAAIyLkAQAAjMigkFdV+6vq/qo6WVU3LDj+4qq6t6pOV9VVC45/ZVV9pKresB6NBgAAYLEVQ15V7Uhyc5IrkuxNck1V7Z2r9uEkr0zy1mVe5ieT/NZTbyYAAABDDBnJuyzJye5+oLufSHJHkgOzFbr7VHd/IMnn5k+uqhck+aokv7EO7QUAAOAMhoS8C5M8OLO/NC1bUVV9SZKfT/LDq28aAAAAq7VzQJ1aUNYDX//7kxzt7gerFr3M9AdUHUpyKEl279498KU5kz033DW47qnDV25gSwAAgM00JOQtJbl4Zv+iJA8NfP1vTPLNVfX9Sb4iyXlV9Xh3f8HDW7r7liS3JMm+ffuGBkgAAADmDAl5x5JcWlWXJPlIkoNJXjHkxbv7u5/crqpXJtk3H/AAAABYPyvek9fdp5Ncn+TuJPcleXt3n6iqm6rqZUlSVS+sqqUkVyf5pao6sZGNBgAAYLEhI3np7qNJjs6V3TizfSyTaZxneo1fTvLLq24hAAAAgw1aDB0AAIBzg5AHAAAwIkIeAADAiAy6J4/VrTsHAABwthjJAwAAGBEjeaxqlPLU4Ss3sCUAAMBaGckDAAAYESEPAABgRIQ8AACAEXFPHqvi/j0AANjajOQBAACMiJAHAAAwIkIeAADAiAh5AAAAIyLkAQAAjIiQBwAAMCJCHgAAwIgIeQAAACNiMXQ2zNCF0y2aDgAA68dIHgAAwIgIeQAAACMi5AEAAIyIkAcAADAiQh4AAMCIeLomZ93Qp3AmnsQJAAArMZIHAAAwIoNCXlXtr6r7q+pkVd2w4PiLq+reqjpdVVfNlD+/qt5TVSeq6gNV9fL1bDwAAABfaMWQV1U7ktyc5Ioke5NcU1V756p9OMkrk7x1rvzTSb6nu5+bZH+SX6yqZ6y10QAAACw25J68y5Kc7O4HkqSq7khyIMmHnqzQ3aemxz43e2J3/8HM9kNV9fEku5L87zW3HAAAgC8yZLrmhUkenNlfmpatSlVdluS8JH+04NihqjpeVccffvjh1b40AAAAU0NCXi0o69X8kKr66iRvTvK93f25+ePdfUt37+vufbt27VrNSwMAADBjSMhbSnLxzP5FSR4a+gOq6iuT3JXkR7v7vatrHgAAAKsx5J68Y0kurapLknwkycEkrxjy4lV1XpJ3JHlTd//aU24lTFlTDwAAzmzFkbzuPp3k+iR3J7kvydu7+0RV3VRVL0uSqnphVS0luTrJL1XVienp35XkxUleWVXvn/73/A35TQAAABg0kpfuPprk6FzZjTPbxzKZxjl/3luSvGWNbQQAAGCgQYuhAwAAcG4Q8gAAAEZEyAMAABgRIQ8AAGBEhDwAAIAREfIAAABGRMgDAAAYkUHr5MG5aM8Ndw2ue+rwlRvYEgAA2DxG8gAAAEbESB7EqB8AAONhJA8AAGBEhDwAAIAREfIAAABGRMgDAAAYESEPAABgRDxdE1bJkzgBANjKjOQBAACMiJAHAAAwIkIeAADAiLgnDzaQ+/cAANhsRvIAAABGRMgDAAAYESEPAABgRIQ8AACAERHyAAAARmTQ0zWran+S1yfZkeSN3X147viLk/xikuclOdjdd84cuzbJj053f6q7b1+PhsPYeBInAADrYcWRvKrakeTmJFck2ZvkmqraO1ftw0lemeStc+c+K8mPJ3lRksuS/HhVPXPtzQYAAGCRIdM1L0tysrsf6O4nktyR5MBshe4+1d0fSPK5uXO/Ncm7uvuR7n40ybuS7F+HdgMAALDAkJB3YZIHZ/aXpmVDDDq3qg5V1fGqOv7www8PfGkAAADmDQl5taCsB77+oHO7+5bu3tfd+3bt2jXwpQEAAJg3JOQtJbl4Zv+iJA8NfP21nAsAAMAqDQl5x5JcWlWXVNV5SQ4mOTLw9e9OcnlVPXP6wJXLp2UAAABsgBVDXnefTnJ9JuHsviRv7+4TVXVTVb0sSarqhVW1lOTqJL9UVSem5z6S5CczCYrHktw0LQMAAGADDFonr7uPJjk6V3bjzPaxTKZiLjr3tiS3raGNAAAADDRkuiYAAADnCCEPAABgRIQ8AACAERHyAAAARkTIAwAAGBEhDwAAYESEPAAAgBEZtE4ecG7ac8Ndg+ueOnzlBrYEAIDNYiQPAABgRIQ8AACAETFdE85Bq5mGCQDA9mIkDwAAYESEPAAAgBER8gAAAEZEyAMAABgRIQ8AAGBEhDwAAIARsYQCkGR1yzKcOnzlBrYEAIC1MJIHAAAwIkbygFUz6gcAsHUZyQMAABgRIQ8AAGBETNcENpSpnQAAm8tIHgAAwIgIeQAAACMyKORV1f6qur+qTlbVDQuOn19Vb5sef19V7ZmWf2lV3V5VH6yq+6rq1evbfAAAAGatGPKqakeSm5NckWRvkmuqau9cteuSPNrdz0nyuiSvnZZfneT87v5rSV6Q5PueDIAAAACsvyEjeZclOdndD3T3E0nuSHJgrs6BJLdPt+9M8tKqqiSd5GlVtTPJn0vyRJJPrUvLAQAA+CJDQt6FSR6c2V+ali2s092nkzyW5IJMAt+fJPlokg8n+bnufmT+B1TVoao6XlXHH3744VX/EgAAAEwMCXm1oKwH1rksyWeTfE2SS5L8i6r6S19UsfuW7t7X3ft27do1oEkAAAAsMiTkLSW5eGb/oiQPLVdnOjXz6UkeSfKKJL/e3X/W3R9P8ttJ9q210QAAACw2JOQdS3JpVV1SVeclOZjkyFydI0munW5fleTd3d2ZTNH82zXxtCTfkOT316fpAAAAzFsx5E3vsbs+yd1J7kvy9u4+UVU3VdXLptVuTXJBVZ1M8s+TPLnMws1JviLJ72USFv9Td39gnX8HAAAApnYOqdTdR5McnSu7cWb7M5kslzB/3uOLygEAANgYg0IewGbYc8Ndg+ueOnzlBrYEAODcNeSePAAAAM4RQh4AAMCImK4JnJNM7QQAWMxIHgAAwIgIeQAAACMi5AEAAIyIkAcAADAiQh4AAMCICHkAAAAjYgkFYPQstwAAbCdG8gAAAEbESB7AjKGjfkb8AICtykgeAADAiAh5AAAAIyLkAQAAjIiQBwAAMCJCHgAAwIgIeQAAACMi5AEAAIyIkAcAADAiFkMH2GBDF1hPLLIOAKydkAfwFKwmuAEAbCbTNQEAAEZEyAMAABiRQSGvqvZX1f1VdbKqblhw/Pyqetv0+Puqas/MsedV1Xuq6kRVfbCqvmz9mg8AAMCsFUNeVe1IcnOSK5LsTXJNVe2dq3Zdkke7+zlJXpfktdNzdyZ5S5JXdfdzk7wkyZ+tW+sBAAD4AkNG8i5LcrK7H+juJ5LckeTAXJ0DSW6fbt+Z5KVVVUkuT/KB7v7dJOnuT3b3Z9en6QAAAMwbEvIuTPLgzP7StGxhne4+neSxJBck+StJuqrurqp7q+pfrr3JAAAALGfIEgq1oKwH1tmZ5G8meWGSTyf5zaq6p7t/8wtOrjqU5FCS7N69e0CTAMbJmnoAwFoNGclbSnLxzP5FSR5ars70PrynJ3lkWv5b3f2J7v50kqNJvn7+B3T3Ld29r7v37dq1a/W/BQAAAEmGhbxjSS6tqkuq6rwkB5McmatzJMm10+2rkry7uzvJ3UmeV1VfPg1/fyvJh9an6QAAAMxbcbpmd5+uquszCWw7ktzW3Seq6qYkx7v7SJJbk7y5qk5mMoJ3cHruo1X1C5kExU5ytLuHz0UCAABgVYbck5fuPprJVMvZshtntj+T5Oplzn1LJssoAAAAsMEGLYYOAADAuWHQSB4AW48ncQIAixjJAwAAGBEhDwAAYESEPAAAgBER8gAAAEbEg1cAtgEPaQGA7cNIHgAAwIgIeQAAACMi5AEAAIyIe/IAeMqG3uvnPj8A2DxG8gAAAEZEyAMAABgR0zUB+AKrWW4BANh6jOQBAACMiJAHAAAwIkIeAADAiAh5AAAAI+LBKwBsuNU8zMWaegCwNkIeAFuKQAgAa2O6JgAAwIgIeQAAACMi5AEAAIyIkAcAADAiQh4AAMCIDAp5VbW/qu6vqpNVdcOC4+dX1dumx99XVXvmju+uqser6ofWp9kAAAAssmLIq6odSW5OckWSvUmuqaq9c9WuS/Jodz8nyeuSvHbu+OuSvHPtzQUAAOBMhqyTd1mSk939QJJU1R1JDiT50EydA0leM92+M8kbqqq6u6vqO5I8kORP1q3VABBr6gHAIkOma16Y5MGZ/aVp2cI63X06yWNJLqiqpyX5kSQ/sfamAgAAsJIhI3m1oKwH1vmJJK/r7serFlWZnlx1KMmhJNm9e/eAJgHA6hj1A2C7GBLylpJcPLN/UZKHlqmzVFU7kzw9ySNJXpTkqqr6mSTPSPK5qvpMd79h9uTuviXJLUmyb9+++QAJAADAQENC3rEkl1bVJUk+kuRgklfM1TmS5Nok70lyVZJ3d3cn+eYnK1TVa5I8Ph/wAAAAWD8rhrzuPl1V1ye5O8mOJLd194mquinJ8e4+kuTWJG+uqpOZjOAd3MhGAwAAsNiQkbx099EkR+fKbpzZ/kySq1d4jdc8hfYBAACwCoMWQwcAAODcMGgkDwC2E0/iBOBcZiQPAABgRIQ8AACAETFdEwDWwNROALYaI3kAAAAjIuQBAACMiOmaALAFDZ0GagooAPOM5AEAAIyIkTwA2CSreUgLADxVRvIAAABGRMgDAAAYESEPAABgRNyTBwDnMIuxAzDPSB4AAMCICHkAAAAjIuQBAACMiHvyAIAv4l4/gHOXkTwAAIARMZIHANvEakbnADh3GckDAAAYESN5AMCauH8PYGsxkgcAADAiQh4AAMCICHkAAAAjMuievKran+T1SXYkeWN3H547fn6SNyV5QZJPJnl5d5+qqm9JcjjJeUmeSPLD3f3udWw/AHAOcf8ewMZbMeRV1Y4kNyf5liRLSY5V1ZHu/tBMteuSPNrdz6mqg0lem+TlST6R5O9290NV9XVJ7k5y4Xr/EgDA+AiEAE/NkOmalyU52d0PdPcTSe5IcmCuzoEkt0+370zy0qqq7v6d7n5oWn4iyZdNR/0AAADYAENC3oVJHpzZX8oXj8Z9vk53n07yWJIL5ur8/SS/091/+tSaCgAAwEqG3JNXC8p6NXWq6rmZTOG8fOEPqDqU5FCS7N69e0CTAAAAWGRIyFtKcvHM/kVJHlqmzlJV7Uzy9CSPJElVXZTkHUm+p7v/aNEP6O5bktySJPv27ZsPkAAAZ+T+PYD/b8h0zWNJLq2qS6rqvCQHkxyZq3MkybXT7auSvLu7u6qekeSuJK/u7t9er0YDAACw2Iojed19uqquz+TJmDuS3NbdJ6rqpiTHu/tIkluTvLmqTmYygndwevr1SZ6T5Meq6semZZd398fX+xcBABhi6KifET/gXDVonbzuPprk6FzZjTPbn0ly9YLzfirJT62xjQAAAAw0ZLomAAAA54hBI3kAANvNah7mshqmgQIbzUgeAADAiAh5AAAAIyLkAQAAjIh78gAANpGF24GNJuQBAGxRAiHwVJiuCQAAMCJG8gAARsCoH/AkIQ8AYJsRCGHcTNcEAAAYESEPAABgREzXBABgWUOndprWCVuHkAcAwJq5zw+2DtM1AQAARsRIHgAAm8qoH2wsI3kAAAAjIuQBAACMiOmaAABsWauZ2rkapoEyZkIeAADbjvsCGTPTNQEAAEbESB4AAJyBUT/ONUIeAACsE4GQrUDIAwCAs2BoIBQGWS0hDwAAtjCjg6zWoJBXVfuTvD7JjiRv7O7Dc8fPT/KmJC9I8skkL+/uU9Njr05yXZLPJvmn3X33urUeAAD4PEtOkAx4umZV7Uhyc5IrkuxNck1V7Z2rdl2SR7v7OUlel+S103P3JjmY5LlJ9if5d9PXAwAAYAMMGcm7LMnJ7n4gSarqjiQHknxops6BJK+Zbt+Z5A1VVdPyO7r7T5P8r6o6OX2996xP8wEAgI22USOEq2E0cbghIe/CJA/O7C8ledFydbr7dFU9luSCafl758698Cm3FgAA2JbcmzjckJBXC8p6YJ0h56aqDiU5NN19vKruH9CuzfTsJJ84243g8/TH1qEvthb9sbXoj61Ff2wt+mNrGV1/1GvPdguesif74mvX8iJDQt5Skotn9i9K8tAydZaqameSpyd5ZOC56e5bktwyvNmbq6qOd/e+s90OJvTH1qEvthb9sbXoj61Ff2wt+mNr0R9bx3r1xYoPXklyLMmlVXVJVZ2XyYNUjszVOZLk2un2VUne3d09LT9YVedX1SVJLk3yP9faaAAAABZbcSRveo/d9UnuzmQJhdu6+0RV3ZTkeHcfSXJrkjdPH6zySCZBMNN6b8/kIS2nk/xAd392g34XAACAbW/QOnndfTTJ0bmyG2e2P5Pk6mXO/ekkP72GNm4FW3Yq6TalP7YOfbG16I+tRX9sLfpja9EfW4v+2DrWpS9qMqsSAACAMRhyTx4AAADnCCFvqqr2V9X9VXWyqm5YcPz8qnrb9Pj7qmrP5rdye6iqi6vqv1fVfVV1oqr+2YI6L6mqx6rq/dP/blz0WqyPqjpVVR+c/l0fX3C8qurfTK+PD1TV15+Ndm4HVfVXZ/6/f39VfaqqfnCujutjA1XVbVX18ar6vZmyZ1XVu6rqD6d/PnOZc6+d1vnDqrp2UR1WZ5n++Nmq+v3p+9E7quoZy5x7xvc2Vm+Z/nhNVX1k5j3p25Y594zfxVi9ZfrjbTN9caqq3r/Mua6PdbTc99uN+vwwXTNJVe1I8gdJviWTZR+OJbmmuz80U+f7kzyvu19VVQeTfGd3v/ysNHjkquqrk3x1d99bVX8+yT1JvmOuP16S5Ie6+9vPUjO3lao6lWRfdy9cQ2f6gf1PknxbkhcleX13v2jzWrg9Td+7PpLkRd39xzPlL4nrY8NU1YuTPJ7kTd39ddOyn0nySHcfnn45fWZ3/8jcec9KcjzJvkzWjL0nyQu6+9FN/QVGZpn+uDyTJ32frpqsljXfH9N6p3KG9zZWb5n+eE2Sx7v7585w3orfxVi9Rf0xd/znkzzW3TctOHYqro91s9z32ySvzAZ8fhjJm7gsycnufqC7n0hyR5IDc3UOJLl9un1nkpdW1aLF3lmj7v5od9873f4/Se5LcuHZbRUrOJDJB0h393uTPGP6ZsbGemmSP5oNeGy87v4fmTxJetbsZ8TtmXxwz/vWJO/q7kemH8zvSrJ/wxq6TSzqj+7+je4+Pd19bybr9LIJlrk+hhjyXYxVOlN/TL/HfleSX93URm1TZ/h+uyGfH0LexIVJHpzZX8oXh4rP15l+cDyW5IJNad02VpNpsX8jyfsWHP7GqvrdqnpnVT13Uxu2/XSS36iqe6rq0ILjQ64h1t/BLP/h7PrYXF/V3R9NJh/kSf7Cgjquk7PjHyZ55zLHVnpvY/1cP50+e9sy09FcH5vvm5N8rLv/cJnjro8NMvf9dkM+P4S8iUUjcvPzWIfUYR1V1Vck+c9JfrC7PzV3+N4kX9vdfz3Jv03yXze7fdvMN3X31ye5IskPTKd/zHJ9bLKqOi/Jy5L82oLDro+tyXWyyarqX2eyTu+vLFNlpfc21se/T/KXkzw/yUeT/PyCOq6PzXdNzjyK5/rYACt8v132tAVlZ7w+hLyJpSQXz+xflOSh5epU1c4kT89Tm47AAFX1pZlcAL/S3f9l/nh3f6q7H59uH03ypVX17E1u5rbR3Q9N//x4kndkMq1m1pBriPV1RZJ7u/tj8wdcH2fFx56cojz98+ML6rhONtH0wQTfnuS7e5kHEAx4b2MddPfHuvuz3f25JP8xi/+eXR+baPpd9u8ledtydVwf62+Z77cb8vkh5E0cS3JpVV0y/dfxg0mOzNU5kuTJJ9lclckN3f6FaQNM54jfmuS+7v6FZer8xSfviayqyzL5f/mTm9fK7aOqnja9QThV9bQklyf5vblqR5J8T018QyY3cX90k5u63Sz7L7Cuj7Ni9jPi2iT/bUGdu5NcXlXPnE5Xu3xaxjqrqv1JfiTJy7r708vUGfLexjqYu0f7O7P473nIdzHWz99J8vvdvbTooOtj/Z3h++2GfH7sXHuTz33Tp29dn8lf1o4kt3X3iaq6Kcnx7j6SSae8uapOZjKCd/DstXj0vinJP0jywZnH+v6rJLuTpLv/QyZB+x9X1ekk/zfJQaF7w3xVkndMM8POJG/t7l+vqlcln++Po5k8WfNkkk8n+d6z1NZtoaq+PJMn0H3fTNlsf7g+NlBV/WqSlyR5dlUtJfnxJIeTvL2qrkvy4SRXT+vuS/Kq7v5H3f1IVf1kJl9mk+Sm7jYjZI2W6Y9XJzk/ybum713vnT4d+2uSvLG7vy3LvLedhV9hVJbpj5dU1fMzmV52KtP3rtn+WO672Fn4FUZlUX90961ZcE+362PDLff9dkM+PyyhAAAAMCKmawIAAIyIkAcAADAiQh4AAMCICHkAAAAjIuQBAACMiJAHAAAwIkIeAADAiAh5AAAAI/L/AG8gppJPnGCtAAAAAElFTkSuQmCC\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=70, density=True);\n",
    "fig=plt.gcf()\n",
    "fig.set_size_inches(15,5);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([11025., 11025., 11025., ..., 11025., 11025., 11025.])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scalar_product_4_np(p1,p1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3MAAAEvCAYAAADvmpjfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAUaklEQVR4nO3dXYjl530f8O8vUuyL5EJStVKEXroq7EWUkjphkQ3phVsXvZauU6pgQZuNbdgWJEggF5Htgopdlw2lSZ3WEai1aglcy6KJ44UoVTYixb2Ro5VrHDuK48XZWlsJSY4Ux0XgIvvXi/PfeLSa95lzZp45nw8MM/P8/+fMc+DZs/Od3/NS3R0AAADG8kN73QEAAAC2TpgDAAAYkDAHAAAwIGEOAABgQMIcAADAgIQ5AACAAV261x1Yz5VXXtmHDx/e624AAADsiWeeeeZb3X1otWv7OswdPnw4Z86c2etuAAAA7Imq+t9rXTPNEgAAYEDCHAAAwICEOQAAgAEJcwAAAAMS5gAAAAYkzAEAAAxImAMAABiQMAcAADAgYQ4AAGBAwhwAAMCAhDkAAIABXbrXHQCAvXT4vt9d9/q5k3cuqCcAsDUqcwAAAAMS5gAAAAYkzAEAAAxImAMAABiQDVAAYB0bbZCyHpunADBPwhwAw7MjJQDLSJgD4MDbSXUNAPYra+YAAAAGJMwBAAAMSJgDAAAYkDAHAAAwIGEOAABgQMIcAADAgIQ5AACAAQlzAAAAA3JoOADMyUaHlZ87eeeCegLAQaQyBwAAMCBhDgAAYECmWQLAHllvGqYpmABsRJgDYAgbrT8DgGVjmiUAAMCAhDkAAIABCXMAAAADsmYOgH3BmjgA2BqVOQAAgAEJcwAAAAMyzRKAhTCNEgB2l8ocAADAgIQ5AACAAQlzAAAAAxLmAAAABmQDFADYhzbaMObcyTsX1BMA9ithDoBdY8dKAFgc0ywBAAAGJMwBAAAMSJgDAAAY0IZr5qrq+iSPJPmxJN9P8mB3f6yqrkjymSSHk5xL8nPd/WpVVZKPJbkjyWtJfqG7vzg91/Ek/3J66n/d3Q/v7ssBgOVggxQANlOZez3JL3f3jyd5R5J7quqmJPclebK7jyR5cvo+SW5PcmT6OJHkgSSZwt/9Sd6e5OYk91fV5bv4WgAAAJbGhmGuu1+4UFnr7u8keTbJtUmOJblQWXs4ybunr48leaRnnkpyWVVdk+TWJKe7+5XufjXJ6SS37eqrAQAAWBJbWjNXVYeT/FSSLyS5urtfSGaBL8lV023XJnluxcPOT21rtV/8M05U1ZmqOvPyyy9vpXsAAABLY9Nhrqp+NMlvJfml7v6r9W5dpa3XaX9jQ/eD3X20u48eOnRos90DAABYKpsKc1X1w5kFuU91929PzS9O0yczfX5paj+f5PoVD78uyfPrtAMAALBFG4a5aXfKTyR5trt/bcWlU0mOT18fT/K5Fe0/XzPvSPLtaRrmE0luqarLp41PbpnaAAAA2KINjyZI8jNJ/lmSP66qL01tH0xyMsljVfX+JN9Mctd07fHMjiU4m9nRBO9Nku5+pao+kuTp6b4Pd/cru/IqAAAAlkx1v2nZ2r5x9OjRPnPmzF53A4BN2ujsM/YHZ9ABjKOqnunuo6td29JulgAAAOwPwhwAAMCAhDkAAIABCXMAAAAD2sxulgCQxAYnALCfqMwBAAAMSJgDAAAYkDAHAAAwIGEOAABgQMIcAADAgOxmCcBfs1slAIxDmAOAJbNRaD938s4F9QSAnRDmAIA3WC/sCXoA+4c1cwAAAAMS5gAAAAYkzAEAAAxImAMAABiQMAcAADAgYQ4AAGBAjiYAWDIOBmcnnFEHsH+ozAEAAAxImAMAABiQMAcAADAgYQ4AAGBAwhwAAMCAhDkAAIABCXMAAAADcs4cwAHjHDkAWA4qcwAAAAMS5gAAAAYkzAEAAAxImAMAABiQDVAABmSTEwBAZQ4AAGBAKnMAwMKsV1U+d/LOBfYEYHzCHACwa0wBBlgcYQ5gD2z0C68KBQCwEWvmAAAABiTMAQAADMg0S4B9yLojAGAjKnMAAAADEuYAAAAGJMwBAAAMyJo5AGBfcGQHwNaozAEAAAxowzBXVQ9V1UtV9ZUVbf+qqv5PVX1p+rhjxbUPVNXZqvpaVd26ov22qe1sVd23+y8FAABgeWymMvfJJLet0v7r3f226ePxJKmqm5K8J8lPTI/5zaq6pKouSfLxJLcnuSnJ3dO9AAAAbMOGa+a6+/NVdXiTz3csyaPd/d0kf15VZ5PcPF07293fSJKqenS690+23GMAAAB2tGbu3qr68jQN8/Kp7dokz6245/zUtlY7AAAA27Dd3SwfSPKRJD19/ndJ3pekVrm3s3po7NWeuKpOJDmRJDfccMM2uwcAHDTr7XZpp0tgGW2rMtfdL3b397r7+0n+U34wlfJ8kutX3HpdkufXaV/tuR/s7qPdffTQoUPb6R4AAMCBt60wV1XXrPj2Z5Nc2OnyVJL3VNVbq+rGJEeS/FGSp5Mcqaobq+otmW2Scmr73QYAAFhuG06zrKpPJ3lnkiur6nyS+5O8s6reltlUyXNJ/nmSdPdXq+qxzDY2eT3JPd39vel57k3yRJJLkjzU3V/d9VcDAACwJDazm+XdqzR/Yp37P5rko6u0P57k8S31DgAAgFXtZDdLAAAA9sh2d7MEANg31tvpMrHbJXAwqcwBAAAMSJgDAAAYkDAHAAAwIGvmAOZkozU8AAA7IcwBrGO9QGZDBQBgLwlzAMDS84cbYETWzAEAAAxIZQ5gm6yJAwD2ksocAADAgFTmgKWmugbLwb914CBSmQMAABiQMAcAADAgYQ4AAGBAwhwAAMCAhDkAAIABCXMAAAADEuYAAAAGJMwBAAAMyKHhAADr2OjA8XMn71xQTwDeSJgDDrSNfgkDABiVaZYAAAADEuYAAAAGJMwBAAAMSJgDAAAYkDAHAAAwIGEOAABgQI4mAACYo/WOSHFGHbATKnMAAAADEuYAAAAGJMwBAAAMyJo5AIAdWG9NHMA8qcwBAAAMSJgDAAAYkGmWwPBMcQIAlpEwB+x7whpwUG30/uYcOmA9plkCAAAMSGUO2BdU3wAAtkZlDgAAYEDCHAAAwICEOQAAgAEJcwAAAAMS5gAAAAYkzAEAAAxImAMAABjQhmGuqh6qqpeq6isr2q6oqtNV9fXp8+VTe1XVb1TV2ar6clX99IrHHJ/u/3pVHZ/PywEAAFgOm6nMfTLJbRe13Zfkye4+kuTJ6fskuT3JkenjRJIHkln4S3J/krcnuTnJ/RcCIAAAAFu3YZjr7s8neeWi5mNJHp6+fjjJu1e0P9IzTyW5rKquSXJrktPd/Up3v5rkdN4cEAEAANikS7f5uKu7+4Uk6e4Xquqqqf3aJM+tuO/81LZWOwAAazh83++uee3cyTsX2BNgP9rtDVBqlbZep/3NT1B1oqrOVNWZl19+eVc7BwAAcFBsN8y9OE2fzPT5pan9fJLrV9x3XZLn12l/k+5+sLuPdvfRQ4cObbN7AAAAB9t2p1meSnI8ycnp8+dWtN9bVY9mttnJt6dpmE8k+TcrNj25JckHtt9tYDTrTRUCAGDrNgxzVfXpJO9McmVVnc9sV8qTSR6rqvcn+WaSu6bbH09yR5KzSV5L8t4k6e5XquojSZ6e7vtwd1+8qQoAAACbtGGY6+6717j0rlXu7ST3rPE8DyV5aEu9A4ah8gYAsFjbnWYJAMAe2ukf0eyGCePb7d0sAQAAWABhDgAAYEDCHAAAwICEOQAAgAEJcwAAAAMS5gAAAAbkaAIAgCW03tEGji2AMajMAQAADEiYAwAAGJBplsCmrTclBwCAxRLmYMkIZAAAB4NplgAAAAMS5gAAAAZkmiUAAG+w0ZR8RxfA/qAyBwAAMCBhDgAAYECmWcIBY7dKAIDloDIHAAAwIGEOAABgQKZZAgCwJetN6bfTJSyOyhwAAMCAhDkAAIABCXMAAAADsmYOAIBds9EROdbUwe5RmQMAABiQMAcAADAg0ywBAFgY0zBh96jMAQAADEiYAwAAGJBpljCYjaanAMDI1vt/zhRMeCNhDgCAA0EQZNkIc7APqb4BALARYQ4AgCH4Yye8kTAHc2DbZQAA5s1ulgAAAAMS5gAAAAZkmiXsAXP+AQDYKZU5AACAAQlzAAAAAxLmAAAABiTMAQAADMgGKLBNNjEBgHE4A5aDSGUOAABgQCpzsAaVNwAA9rMdhbmqOpfkO0m+l+T17j5aVVck+UySw0nOJfm57n61qirJx5LckeS1JL/Q3V/cyc8HAIDdYBomI9qNaZZ/r7vf1t1Hp+/vS/Jkdx9J8uT0fZLcnuTI9HEiyQO78LMBAACW0jzWzB1L8vD09cNJ3r2i/ZGeeSrJZVV1zRx+PgAAwIG30zDXSX6/qp6pqhNT29Xd/UKSTJ+vmtqvTfLciseen9oAAADYop1ugPIz3f18VV2V5HRV/ek699Yqbf2mm2ah8ESS3HDDDTvsHgAAwMG0o8pcdz8/fX4pyWeT3JzkxQvTJ6fPL023n09y/YqHX5fk+VWe88HuPtrdRw8dOrST7gEAABxY267MVdWPJPmh7v7O9PUtST6c5FSS40lOTp8/Nz3kVJJ7q+rRJG9P8u0L0zFhrzh+AACAUe1kmuXVST47O3Eglyb5r93936vq6SSPVdX7k3wzyV3T/Y9ndizB2cyOJnjvDn42AADAUtt2mOvubyT5O6u0/0WSd63S3knu2e7PAwCAvbLebB5n0LFX5nE0AQAAAHO2090sAQBgqW20Bl/ljnkR5jjQbHACAMBBZZolAADAgFTmGJrKGwAAy0qYY98T2AAA4M2EOQAAmCPHGjAv1swBAAAMSJgDAAAYkDAHAAAwIGEOAABgQDZAAQCAPbLRrt02SGE9KnMAAAADUpljzzlHDgBgdSp3rEdlDgAAYEDCHAAAwIBMs2QhTKUEAIDdpTIHAAAwIJU5doXKGwDA4q33O5jNUQ4+lTkAAIABCXMAAAADEuYAAAAGZM0cm2JNHAAA7C/CHAAALKGN/lhvA5X9zzRLAACAAQlzAAAAAzLNEgAADiB7Hhx8whx/zT94AAAYh2mWAAAAA1KZWyIqbwAAcHAIcwAAwJvspBDgWIPFEOYOGNU3AABYDtbMAQAADEiYAwAAGJBplgAAwK5ab+mP9XS7R5gbjDVxAABAYpolAADAkIQ5AACAAZlmCQAALMxGy4asqds8YW6fsSYOAADYDNMsAQAABqQytwdU3wAAYHWmYW6eyhwAAMCAVOYAAIBhOJD8B4S5OTCNEgAAmDfTLAEAAAa08MpcVd2W5GNJLknyn7v75KL7AAAAHDzLtnnKQitzVXVJko8nuT3JTUnurqqbFtkHAACAg2DRlbmbk5zt7m8kSVU9muRYkj9ZcD92xJo4AAAYz0HbPGXRa+auTfLciu/PT20AAABswaIrc7VKW7/hhqoTSU5M3/7fqvra3HvFolyZ5Ft73QkOJGOLeTG2mAfjinkxtnagfnWve7Cmv7nWhUWHufNJrl/x/XVJnl95Q3c/mOTBRXaKxaiqM919dK/7wcFjbDEvxhbzYFwxL8bW8ln0NMunkxypqhur6i1J3pPk1IL7AAAAMLyFVua6+/WqujfJE5kdTfBQd391kX0AAAA4CBZ+zlx3P57k8UX/XPYF02eZF2OLeTG2mAfjinkxtpZMdffGdwEAALCvLHrNHAAAALtAmGNXVNW/rao/raovV9Vnq+qyFdc+UFVnq+prVXXrivbbprazVXXfivYbq+oLVfX1qvrMtFkOS6qq7qqqr1bV96vq6EXXjC3mYq0xBGupqoeq6qWq+sqKtiuq6vT0nnO6qi6f2quqfmMaX1+uqp9e8Zjj0/1fr6rje/Fa2D+q6vqq+sOqenb6v/AXp3ZjiyTCHLvndJK/3d0/meTPknwgSarqpsx2Lf2JJLcl+c2quqSqLkny8SS3J7kpyd3TvUnyq0l+vbuPJHk1yfsX+krYb76S5B8n+fzKRmOLedlgDMFaPpnZe9FK9yV5cnrPeXL6PpmNrSPTx4kkDySzX9CT3J/k7UluTnL/hV/SWVqvJ/nl7v7xJO9Ics/0fmRskUSYY5d09+939+vTt09ldoZgkhxL8mh3f7e7/zzJ2czeRG5Ocra7v9Hd/y/Jo0mOVVUl+ftJ/tv0+IeTvHtRr4P9p7uf7e6vrXLJ2GJeVh1De9wn9rnu/nySVy5qPpbZe03yxvecY0ke6ZmnklxWVdckuTXJ6e5+pbtfzewPpRcHRJZId7/Q3V+cvv5OkmeTXBtji4kwxzy8L8nvTV9fm+S5FdfOT21rtf+NJH+5IhheaIeLGVvMy1pjCLbq6u5+IZn9Up7kqql9q+9fkKo6nOSnknwhxhaThR9NwLiq6g+S/Ngqlz7U3Z+b7vlQZlMCPnXhYavc31n9Dwm9zv0cYJsZW6s9bJU2Y4vdYKwwb2uNMWOPVVXVjyb5rSS/1N1/NZtssvqtq7QZWweYMMemdfc/WO/6tJj2HyZ5V//gzIvzSa5fcdt1SZ6fvl6t/VuZTQm4dKqgrLyfA2qjsbUGY4t5WW9swVa8WFXXdPcL01S3l6b2tcbY+STvvKj9fyygn+xjVfXDmQW5T3X3b0/NxhZJTLNkl1TVbUl+Jck/6u7XVlw6leQ9VfXWqroxswW5f5Tk6SRHpt0F35LZRhanphD4h0n+yfT440nWqsyw3Iwt5mXVMbTHfWJMpzJ7r0ne+J5zKsnPTzsPviPJt6epck8kuaWqLp82p7hlamNJTeu9P5Hk2e7+tRWXjC2SODScXVJVZ5O8NclfTE1Pdfe/mK59KLN1dK9nNj3g96b2O5L8+ySXJHmouz86tf+tzDYcuCLJ/0ryT7v7uwt8OewjVfWzSf5DkkNJ/jLJl7r71umascVcrDWGYC1V9enMKh9XJnkxs50DfyfJY0luSPLNJHd19yvTL+j/MbMNKF5L8t7uPjM9z/uSfHB62o92939Z5Otgf6mqv5vkfyb54yTfn5o/mNm6OWMLYQ4AAGBEplkCAAAMSJgDAAAYkDAHAAAwIGEOAABgQMIcAADAgIQ5AACAAQlzAAAAAxLmAAAABvT/AaxsM9MYXFaJAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(p1[:,1],bins=100);\n",
    "fig = plt.gcf()\n",
    "fig.set_size_inches(15,5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3MAAAEvCAYAAADvmpjfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAYU0lEQVR4nO3dfYyl1X0f8O8vEJM6UQI2Y4fsLl3SbNMQq5XRCqOmqiyTYDCW161MZffFWxtpFRW3Tp0oXuw/qBK5WitVSFy7SMQQg2SBXccpq0LqbIgtt1IgrF+CwYSywluYQMy6YJIWxQ71r3/cZ+1hmZ3ZnTtvz9zPRxrN85zn3HvPZY8u873nPOdUdwcAAIBx+b6NbgAAAACnT5gDAAAYIWEOAABghIQ5AACAERLmAAAARkiYAwAAGKEzN7oBSzn33HN7586dG90MAACADfGFL3zhG909t9i1TR3mdu7cmcOHD290MwAAADZEVf2vk10zzRIAAGCEhDkAAIAREuYAAABGSJgDAAAYIWEOAABghIQ5AACAERLmAAAARkiYAwAAGCFhDgAAYISWDXNVdXNVPVVVDyxy7Zeqqqvq3OG8qupDVXWkqu6vqosW1N1bVY8MP3tX920AAADMllMZmftYkstPLKyqHUl+LsljC4qvSLJr+NmX5Iah7suSXJfkNUkuTnJdVZ0zTcMBAABm2ZnLVejuz1fVzkUuXZ/kl5PcsaBsT5Jbu7uT3FNVZ1fVeUlem+RQdz+dJFV1KJOAeNtUrQeYUTv333nSa0cPXLmOLQEANsqK7pmrqjcl+bPu/pMTLm1L8viC8/mh7GTlAAAArMCyI3MnqqqXJnl/kssWu7xIWS9Rvtjz78tkimbOP//8020eAADATFjJyNzfSnJBkj+pqqNJtif5YlX9aCYjbjsW1N2e5Iklyl+ku2/s7t3dvXtubm4FzQMAANj6TjvMdfdXuvsV3b2zu3dmEtQu6u4/T3IwyduHVS0vSfJsdz+Z5DNJLquqc4aFTy4bygAAAFiBU9ma4LYkf5TkJ6tqvqquXqL6XUkeTXIkyW8l+VdJMix88qtJ7ht+fuX4YigAAACcvlNZzfJty1zfueC4k1xzkno3J7n5NNsHADNpqRVLE6uWArDC1SwBAADYWMIcAADACAlzAAAAI3Ta+8wBsLm512rz8G8BwFoyMgcAADBCRuYAOGVGml5suf8mALBWhDkARmGp0DSLIRIAhDkANoW1HOESBNfPRo3eGjUGZpEwBwBsCgIZwOkR5gCYae55A2CshDlg9EyhAwBmkTAHwHcZpWI50/aRzdrHfCkEjJEwB2xp7sGZDZs1IKwlfRsAYQ5gjfimHwBYS8IcALDlzeLoLbD1fd9GNwAAAIDTZ2QOYANs5P1ORig2D/8WAEzDyBwAAMAIGZkDgC3IAjwAW58wB7CEjZoGtxWn323F97RV+bcCGAdhDoBVIwSMw1j/ncbaboC14p45AACAERLmAAAARsg0SwCADbKR25QA4yfMAZuClfcAAE6PMAfMNAsqAMvZrKNnm7VdwPoR5gAApiBUARtFmAMAWEPTzAAwewBYyrKrWVbVzVX1VFU9sKDs16rqT6vq/qr63ao6e8G1a6vqSFU9XFWvX1B++VB2pKr2r/5bAQAAmB2nMjL3sSQfTnLrgrJDSa7t7uer6oNJrk3y3qq6MMlbk/x0kh9L8gdV9beHx3wkyc8lmU9yX1Ud7O6vrs7bALYy30wDrC9TR2Eclh2Z6+7PJ3n6hLLf7+7nh9N7kmwfjvckub27v9XdX0tyJMnFw8+R7n60u7+d5PahLgAAACuwGvfMvTPJJ4bjbZmEu+Pmh7IkefyE8tcs9mRVtS/JviQ5//zzV6F5AACzx+gabH3Ljswtparen+T5JB8/XrRItV6i/MWF3Td29+7u3j03NzdN8wAAALasFY/MVdXeJG9Mcml3Hw9m80l2LKi2PckTw/HJygEAADhNKxqZq6rLk7w3yZu6+7kFlw4meWtVnVVVFyTZleSPk9yXZFdVXVBVL8lkkZSD0zUdAABgdi07MldVtyV5bZJzq2o+yXWZrF55VpJDVZUk93T3z3f3g1X1ySRfzWT65TXd/f+G53lXks8kOSPJzd394Bq8H2CTsiIlAMDqWjbMdffbFim+aYn6H0jygUXK70py12m1DgAAgEWtxmqWwBZh5TMAgPEQ5gAAZpDp7zB+U21NAAAAwMYQ5gAAAEbINEvglC01Jcf9dAAA60uYAwDgtExzv50v/2D1mGYJAAAwQsIcAADACJlmCVuMpaYBAGaDMAeb0BgXGhEiAQDWl2mWAAAAIyTMAQAAjJAwBwAAMELCHAAAwAgJcwAAACMkzAEAAIyQrQkAAFg3y21ls1m34IHNSJiDDTDNnmz2cwMAIDHNEgAAYJSEOQAAgBES5gAAAEZImAMAABghYQ4AAGCEhDkAAIAREuYAAABGyD5zAACMgg3H4YWMzAEAAIzQsmGuqm6uqqeq6oEFZS+rqkNV9cjw+5yhvKrqQ1V1pKrur6qLFjxm71D/karauzZvBwAAYDacysjcx5JcfkLZ/iR3d/euJHcP50lyRZJdw8++JDckk/CX5Lokr0lycZLrjgdAAAAATt+yYa67P5/k6ROK9yS5ZTi+JcmbF5Tf2hP3JDm7qs5L8vokh7r76e5+JsmhvDggAgAAcIpWes/cK7v7ySQZfr9iKN+W5PEF9eaHspOVAwAAsAKrvQBKLVLWS5S/+Amq9lXV4ao6fOzYsVVtHAAAwFax0jD39WH6ZIbfTw3l80l2LKi3PckTS5S/SHff2N27u3v33NzcCpsHAACwta10n7mDSfYmOTD8vmNB+buq6vZMFjt5trufrKrPJPn3CxY9uSzJtStvNmxuy+2DAwAA01o2zFXVbUlem+TcqprPZFXKA0k+WVVXJ3ksyVVD9buSvCHJkSTPJXlHknT301X1q0nuG+r9SnefuKgKjIrABgDARlo2zHX3205y6dJF6naSa07yPDcnufm0WgcAwEzxZSmcutVeAAUAAIB1sNJ75gAAYFNZalTv6IEr17ElsD6MzAEAAIyQMAcAADBCplkCALDlLbewimmYjJEwBydhNS0AADYz0ywBAABGSJgDAAAYIdMsmWmmUgIAMFbCHAAALMMedmxGplkCAACMkJE5AABmnlsvGCMjcwAAACMkzAEAAIyQMAcAADBCwhwAAMAICXMAAAAjZDVLAABYQ/aoY60IcwAAMAXbGrBRTLMEAAAYIWEOAABghIQ5AACAERLmAAAARsgCKGxpbkgGAGCrMjIHAAAwQsIcAADACAlzAAAAIyTMAQAAjNBUYa6q/m1VPVhVD1TVbVX1A1V1QVXdW1WPVNUnquolQ92zhvMjw/Wdq/EGAAAAZtGKw1xVbUvyb5Ls7u5XJTkjyVuTfDDJ9d29K8kzSa4eHnJ1kme6+yeSXD/UAwAAYAWm3ZrgzCR/o6r+OslLkzyZ5HVJ/ulw/ZYk/y7JDUn2DMdJ8qkkH66q6u6esg3MMFsPAABjttzfMkcPXLlOLWGMVjwy191/luQ/JHkskxD3bJIvJPlmdz8/VJtPsm043pbk8eGxzw/1X77S1wcAAJhl00yzPCeT0bYLkvxYkh9McsUiVY+PvNUS1xY+776qOlxVh48dO7bS5gEAAGxp0yyA8rNJvtbdx7r7r5N8OsnfT3J2VR2fvrk9yRPD8XySHUkyXP+RJE+f+KTdfWN37+7u3XNzc1M0DwAAYOuaJsw9luSSqnppVVWSS5N8Nclnk7xlqLM3yR3D8cHhPMP1P3S/HAAAwMpMc8/cvZksZPLFJF8ZnuvGJO9N8p6qOpLJPXE3DQ+5KcnLh/L3JNk/RbsBAABm2lSrWXb3dUmuO6H40SQXL1L3r5JcNc3rAQAAMDHVpuEAAABsDGEOAABghKbdNBzWnI3BAQDgxYzMAQAAjJAwBwAAMELCHAAAwAgJcwAAACNkARQAANikllsI7uiBK9epJWxGRuYAAABGyMgc68L2AgAAsLqMzAEAAIyQMAcAADBCwhwAAMAICXMAAAAjJMwBAACMkNUsAQBgpJZaMdwedFufkTkAAIAREuYAAABGSJgDAAAYIWEOAABghIQ5AACAERLmAAAARkiYAwAAGCH7zAEAwAxaao+6xD51Y2BkDgAAYISMzLEqlvtmBwAAWF1G5gAAAEZImAMAABihqaZZVtXZST6a5FVJOsk7kzyc5BNJdiY5muSfdPczVVVJfjPJG5I8l+RfdvcXp3l9AABgcW6D2fqmHZn7zST/rbv/TpK/l+ShJPuT3N3du5LcPZwnyRVJdg0/+5LcMOVrAwAAzKwVh7mq+uEk/zDJTUnS3d/u7m8m2ZPklqHaLUnePBzvSXJrT9yT5OyqOm/FLQcAAJhh04zM/XiSY0l+u6q+VFUfraofTPLK7n4ySYbfrxjqb0vy+ILHzw9lAAAAnKZpwtyZSS5KckN3vzrJ/833plQuphYp6xdVqtpXVYer6vCxY8emaB4AAMDWNU2Ym08y3933DuefyiTcff349Mnh91ML6u9Y8PjtSZ448Um7+8bu3t3du+fm5qZoHgAAwNa14jDX3X+e5PGq+smh6NIkX01yMMneoWxvkjuG44NJ3l4TlyR59vh0TAAAAE7PVFsTJPnXST5eVS9J8miSd2QSED9ZVVcneSzJVUPduzLZluBIJlsTvGPK1wYAAJhZU4W57v5ykt2LXLp0kbqd5JppXg8AAICJafeZAwAAYAMIcwAAACM07T1zzIid++/c6CYAAAALGJkDAAAYIWEOAABghIQ5AACAERLmAAAARkiYAwAAGCFhDgAAYISEOQAAgBES5gAAAEbIpuF8l43BAQA4bqm/DY8euHIdW8LJGJkDAAAYIWEOAABghIQ5AACAERLmAAAARkiYAwAAGCFhDgAAYISEOQAAgBES5gAAAEbIpuEAAMBpWWpD8cSm4uvFyBwAAMAIGZkDAABW1VIjd0btVo+ROQAAgBES5gAAAEZImAMAABghYQ4AAGCELIAyQ5ZbQhYAABiPqUfmquqMqvpSVf3X4fyCqrq3qh6pqk9U1UuG8rOG8yPD9Z3TvjYAAMCsWo1plu9O8tCC8w8mub67dyV5JsnVQ/nVSZ7p7p9Icv1QDwAAgBWYKsxV1fYkVyb56HBeSV6X5FNDlVuSvHk43jOcZ7h+6VAfAACA0zTtyNxvJPnlJN8Zzl+e5Jvd/fxwPp9k23C8LcnjSTJcf3ao/wJVta+qDlfV4WPHjk3ZPAAAgK1pxWGuqt6Y5Knu/sLC4kWq9ilc+15B943dvbu7d8/Nza20eQAAAFvaNKtZ/kySN1XVG5L8QJIfzmSk7uyqOnMYfdue5Imh/nySHUnmq+rMJD+S5OkpXh8AAGBmrTjMdfe1Sa5Nkqp6bZJf6u5/VlX/OclbktyeZG+SO4aHHBzO/2i4/ofd/aKROaZj+wEAAJgNa7Fp+HuTvKeqjmRyT9xNQ/lNSV4+lL8nyf41eG0AAICZsCqbhnf355J8bjh+NMnFi9T5qyRXrcbrAQAAzLq1GJkDAABgja3KyBwAAMCpWG6Nh6MHrlynloyfkTkAAIARMjI3MlarBAAAEiNzAAAAoyTMAQAAjJAwBwAAMELCHAAAwAgJcwAAACMkzAEAAIyQMAcAADBCwhwAAMAICXMAAAAjJMwBAACMkDAHAAAwQsIcAADACAlzAAAAIyTMAQAAjNCZG90AAACA43buv3PJ60cPXLlOLdn8hLlNZrnOCwAAkJhmCQAAMErCHAAAwAgJcwAAACMkzAEAAIyQMAcAADBCwhwAAMAICXMAAAAjtOIwV1U7quqzVfVQVT1YVe8eyl9WVYeq6pHh9zlDeVXVh6rqSFXdX1UXrdabAAAAmDXTjMw9n+QXu/unklyS5JqqujDJ/iR3d/euJHcP50lyRZJdw8++JDdM8doAAAAzbcVhrruf7O4vDsd/meShJNuS7Elyy1DtliRvHo73JLm1J+5JcnZVnbfilgMAAMywVblnrqp2Jnl1knuTvLK7n0wmgS/JK4Zq25I8vuBh80MZAAAAp+nMaZ+gqn4oye8k+YXu/ouqOmnVRcp6kefbl8k0zJx//vnTNg8AANhCdu6/86TXjh64ch1bsvGmGpmrqu/PJMh9vLs/PRR//fj0yeH3U0P5fJIdCx6+PckTJz5nd9/Y3bu7e/fc3Nw0zQMAANiyVjwyV5MhuJuSPNTdv77g0sEke5McGH7fsaD8XVV1e5LXJHn2+HTMWbPUtwkAAACnYppplj+T5F8k+UpVfXkoe18mIe6TVXV1kseSXDVcuyvJG5IcSfJckndM8doAAAAzbcVhrrv/Rxa/Dy5JLl2kfie5ZqWvBwAAwPesymqWAAAArC9hDgAAYISEOQAAgBES5gAAAEZImAMAABghYQ4AAGCEhDkAAIAREuYAAABGSJgDAAAYIWEOAABghIQ5AACAETpzoxuwFe3cf+dGNwEAAGbOcn+HHz1w5Tq1ZH0YmQMAABghYQ4AAGCEhDkAAIAREuYAAABGSJgDAAAYIWEOAABghIQ5AACAERLmAAAARkiYAwAAGCFhDgAAYISEOQAAgBES5gAAAEbozI1uwBjt3H/nRjcBAACYccIcAAAwE5YalDl64Mp1bMnqMM0SAABghIQ5AACAEVr3MFdVl1fVw1V1pKr2r/frAwAAbAXrGuaq6owkH0lyRZILk7ytqi5czzYAAABsBes9MndxkiPd/Wh3fzvJ7Un2rHMbAAAARm+9w9y2JI8vOJ8fygAAADgN6701QS1S1i+oULUvyb7h9P9U1cNr3io20rlJvrHRjWBL0adYbfoUq02fYrXpU6ugPrjRLTipv3myC+sd5uaT7Fhwvj3JEwsrdPeNSW5cz0axcarqcHfv3uh2sHXoU6w2fYrVpk+x2vSp2bXe0yzvS7Krqi6oqpckeWuSg+vcBgAAgNFb15G57n6+qt6V5DNJzkhyc3c/uJ5tAAAA2ArWe5pluvuuJHet9+uyaZlSy2rTp1ht+hSrTZ9itelTM6q6e/laAAAAbCrrfc8cAAAAq0CYY81U1a9V1Z9W1f1V9btVdfaCa9dW1ZGqeriqXr+g/PKh7EhV7V9QfkFV3VtVj1TVJ4YFdJgxVXVVVT1YVd+pqt0nXNOnWFUn6ztwoqq6uaqeqqoHFpS9rKoODZ8xh6rqnKG8qupDQ7+6v6ouWvCYvUP9R6pq70a8FzaHqtpRVZ+tqoeG/++9eyjXr3gBYY61dCjJq7r77yb5n0muTZKqujCTlUx/OsnlSf5TVZ1RVWck+UiSK5JcmORtQ90k+WCS67t7V5Jnkly9ru+EzeKBJP84yecXFupTrLZl+g6c6GOZfPYstD/J3cNnzN3DeTLpU7uGn31Jbkgmf6QnuS7Ja5JcnOS643+oM5OeT/KL3f1TSS5Jcs3wGaRf8QLCHGumu3+/u58fTu/JZF/BJNmT5Pbu/lZ3fy3JkUw+YC5OcqS7H+3ubye5Pcmeqqokr0vyqeHxtyR583q9DzaP7n6oux9e5JI+xWpbtO9scJvYpLr780mePqF4TyafLckLP2P2JLm1J+5JcnZVnZfk9UkOdffT3f1MJl+InhgQmRHd/WR3f3E4/sskDyXZFv2KEwhzrJd3Jvm94XhbkscXXJsfyk5W/vIk31wQDI+Xw3H6FKvtZH0HTtUru/vJZPKHeZJXDOWn+3nFjKuqnUleneTe6FecYN23JmBrqao/SPKji1x6f3ffMdR5fybTBT5+/GGL1O8s/uVCL1GfLehU+tRiD1ukTJ9iGvoIa+VkfUuf40Wq6oeS/E6SX+juv5hMLFm86iJl+tUMEOaYSnf/7FLXhxtt35jk0v7ePhjzSXYsqLY9yRPD8WLl38hkusCZw0jKwvpsMcv1qZPQp1htS/UpOBVfr6rzuvvJYbrbU0P5yfrWfJLXnlD+uXVoJ5tUVX1/JkHu49396aFYv+IFTLNkzVTV5Unem+RN3f3cgksHk7y1qs6qqgsyuVn3j5Pcl2TXsMrgSzJZ0OLgEAI/m+Qtw+P3JjnZCA2zSZ9itS3adza4TYzLwUw+W5IXfsYcTPL2YfXBS5I8O0yX+0ySy6rqnGGBisuGMmbQcG/3TUke6u5fX3BJv+IFbBrOmqmqI0nOSvK/h6J7uvvnh2vvz+Q+uuczmTrwe0P5G5L8RpIzktzc3R8Yyn88kwUIXpbkS0n+eXd/ax3fDptAVf2jJP8xyVySbyb5cne/frimT7GqTtZ34ERVdVsmox/nJvl6JqsH/pckn0xyfpLHklzV3U8Pf6R/OJNFKJ5L8o7uPjw8zzuTvG942g9092+v5/tg86iqf5Dkvyf5SpLvDMXvy+S+Of2K7xLmAAAARsg0SwAAgBES5gAAAEZImAMAABghYQ4AAGCEhDkAAIAREuYAAABGSJgDAAAYIWEOAABghP4/dmthEAo3GbcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(pKst[:,1],bins=100);\n",
    "fig = plt.gcf()\n",
    "fig.set_size_inches(15,5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "p1_q = lorentz_boost_np(p1, p1+p2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "105.0"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sqrt(scalar_product_4_np(p1_q,p1_q)).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "z=np.array([[0.,0.,1.,0] for i in range(n_events)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAASmElEQVR4nO3dfaxkdX3H8fdXtkDV6C6w0HUXXIhbldhUyA1STXxgjQIalqZg19S60m02WrS2tClLNbGxT9A0pRob7FbQpTUIrhq2FUtWHmKaCHVRBGGLe0EL113Za3mwloig3/4xvyuHu3P3zp0zM/fh934lm3vO7/zOOd/5zZnPnDnzsJGZSJLq8Jz5LkCSNDqGviRVxNCXpIoY+pJUEUNfkiqybL4LOJRjjjkm165dO99lSNKicscdd/wgM1d2W7agQ3/t2rXs3r17vsuQpEUlIv57pmVe3pGkihj6klQRQ1+SKmLoS1JFDH1JqoihL0kVMfQlqSKGviRVxNCXpIos6G/kSrNZu/WLP5/+7qVvmcdKNCre5+14pi9JFTH0JakiXt4ZMl+KSkvfYnqcG/oN83XHzecBM+x9L6YHg0an9uNiPm+/oa++jfLAHda+ag8f1cfQV1dLNQybt6tGi+l+XUy1LiaGvmblg290hjHWtT/RzdVSP96XdOgv9TtvLhbjWMwUVqMIscU4XsPiWCwtSzr0tXAthSAZ1G1YCGfi02tYrPfJbJbCcdfWrKEfEVcBbwUOZOYrSttRwLXAWuC7wNsy89GICOAjwNnAE8C7MvPrZZ1NwAfLZv8iM7cP9qZI86fGT35pcerlTP9TwMeAqxttW4GbMvPSiNha5i8GzgLWlX+vAq4AXlWeJD4EjAEJ3BEROzPz0UHdkIVk2GduC/GBvhBr0uLgsTNas4Z+Zn4lItZOa94AvL5MbwdupRP6G4CrMzOB2yJieUSsKn13ZeYjABGxCzgTuKb1LVgg5usl+kz7HdaDZyFciljohnGfLMQ3eA3rxanfa/rHZeZ+gMzcHxHHlvbVwEONfhOlbab2g0TEFmALwAknnNBneb0zxObfMK6NL/TP8nvcab4M+o3c6NKWh2g/uDFzG7ANYGxsrGuftnp5wPXyQG/zwPVBr8VorsftUnqzG5bGq5t+Q//hiFhVzvJXAQdK+wRwfKPfGmBfaX/9tPZb+9z3yC21A3e+jPL2j+Ksf5Tran4sxfus39DfCWwCLi1/r2+0vzciPkPnjdzHyxPDjcBfRcSK0u9NwCX9l62laik+yFSXhf5qoJePbF5D5yz9mIiYoPMpnEuB6yJiM/AgcH7pfgOdj2uO0/nI5gUAmflIRPw58LXS78NTb+rq2Rb6ASMtBnN9HNV0stHLp3fePsOi9V36JnDhDNu5CrhqTtUtQDUdHNIgDeOxM+z33tqaa32jONGr5hu5hrUWG4/Z+oziCaCa0Nczhv2bNobVcM3102ejsBAuS7b5ZFFNDP1FotYDVPNjvj6l5HE+fIa+lrxagqSW26l2DP0FzAexpEEz9KXCJ1nVwNCXpD4s1pMEQ3+JWqwHpKThMvQlaQQWyonYc+a7AEnS6Bj6klQRQ1+SKmLoS1JFDH1JqoihL0kVMfQlqSKGviRVxNCXpIoY+pJUEUNfkipi6EtSRQx9SaqIoS9JFTH0Jakihr4kVcTQl6SKGPqSVBFDX5IqYuhLUkUMfUmqSKvQj4g/jIh7IuJbEXFNRBwZESdGxO0RsTciro2Iw0vfI8r8eFm+dhA3QJLUu75DPyJWA78PjGXmK4DDgI3AZcDlmbkOeBTYXFbZDDyamS8BLi/9JEkj1PbyzjLgFyNiGfBcYD9wBrCjLN8OnFumN5R5yvL1EREt9y9JmoO+Qz8zvwf8LfAgnbB/HLgDeCwzny7dJoDVZXo18FBZ9+nS/+jp242ILRGxOyJ2T05O9lueJKmLNpd3VtA5ez8ReBHwPOCsLl1zapVDLHumIXNbZo5l5tjKlSv7LU+S1EWbyztvBL6TmZOZ+RTweeDVwPJyuQdgDbCvTE8AxwOU5S8EHmmxf0nSHLUJ/QeB0yPiueXa/HrgXuAW4LzSZxNwfZneWeYpy2/OzIPO9CVJw9Pmmv7tdN6Q/Tpwd9nWNuBi4KKIGKdzzf7KssqVwNGl/SJga4u6JUl9WDZ7l5ll5oeAD01rfgA4rUvfHwPnt9mfJKkdv5ErSRUx9CWpIoa+JFXE0Jekihj6klQRQ1+SKmLoS1JFDH1JqoihL0kVMfQlqSKGviRVxNCXpIoY+pJUEUNfkipi6EtSRQx9SaqIoS9JFTH0Jakihr4kVcTQl6SKGPqSVBFDX5IqYuhLUkUMfUmqiKEvSRUx9CWpIoa+JFXE0Jekihj6klSRVqEfEcsjYkdE/FdE7ImIX4uIoyJiV0TsLX9XlL4RER+NiPGIuCsiTh3MTZAk9artmf5HgH/PzJcBvwrsAbYCN2XmOuCmMg9wFrCu/NsCXNFy35KkOeo79CPiBcBrgSsBMvMnmfkYsAHYXrptB84t0xuAq7PjNmB5RKzqu3JJ0py1OdM/CZgEPhkR34iIT0TE84DjMnM/QPl7bOm/Gniosf5EaZMkjUib0F8GnApckZmnAP/HM5dyuokubXlQp4gtEbE7InZPTk62KE+SNF2b0J8AJjLz9jK/g86TwMNTl23K3wON/sc31l8D7Ju+0czclpljmTm2cuXKFuVJkqbrO/Qz8/vAQxHx0tK0HrgX2AlsKm2bgOvL9E7gneVTPKcDj09dBpIkjcayluu/D/h0RBwOPABcQOeJ5LqI2Aw8CJxf+t4AnA2MA0+UvpKkEWoV+pl5JzDWZdH6Ln0TuLDN/iRJ7fiNXEmqiKEvSRUx9CWpIoa+JFXE0Jekihj6klQRQ1+SKmLoS1JFDH1JqoihL0kVMfQlqSKGviRVxNCXpIoY+pJUEUNfkipi6EtSRQx9SaqIoS9JFTH0Jakihr4kVcTQl6SKGPqSVBFDX5IqYuhLUkUMfUmqiKEvSRUx9CWpIoa+JFXE0Jekihj6klSR1qEfEYdFxDci4t/K/IkRcXtE7I2IayPi8NJ+RJkfL8vXtt23JGluBnGm/35gT2P+MuDyzFwHPApsLu2bgUcz8yXA5aWfJGmEWoV+RKwB3gJ8oswHcAawo3TZDpxbpjeUecry9aW/JGlE2p7p/z3wJ8DPyvzRwGOZ+XSZnwBWl+nVwEMAZfnjpf+zRMSWiNgdEbsnJydblidJauo79CPircCBzLyj2dyla/aw7JmGzG2ZOZaZYytXruy3PElSF8tarPsa4JyIOBs4EngBnTP/5RGxrJzNrwH2lf4TwPHAREQsA14IPNJi/5KkOer7TD8zL8nMNZm5FtgI3JyZvwXcApxXum0Cri/TO8s8ZfnNmXnQmb4kaXiG8Tn9i4GLImKczjX7K0v7lcDRpf0iYOsQ9i1JOoQ2l3d+LjNvBW4t0w8Ap3Xp82Pg/EHsT5LUH7+RK0kVMfQlqSKGviRVxNCXpIoY+pJUEUNfkipi6EtSRQx9SaqIoS9JFTH0Jakihr4kVcTQl6SKGPqSVBFDX5IqYuhLUkUMfUmqiKEvSRUx9CWpIoa+JFXE0Jekihj6klQRQ1+SKmLoS1JFDH1JqoihL0kVMfQlqSKGviRVxNCXpIoY+pJUEUNfkirSd+hHxPERcUtE7ImIeyLi/aX9qIjYFRF7y98VpT0i4qMRMR4Rd0XEqYO6EZKk3rQ5038a+KPMfDlwOnBhRJwMbAVuysx1wE1lHuAsYF35twW4osW+JUl96Dv0M3N/Zn69TP8vsAdYDWwAtpdu24Fzy/QG4OrsuA1YHhGr+q5ckjRnA7mmHxFrgVOA24HjMnM/dJ4YgGNLt9XAQ43VJkrb9G1tiYjdEbF7cnJyEOVJkorWoR8Rzwc+B/xBZv7wUF27tOVBDZnbMnMsM8dWrlzZtjxJUkOr0I+IX6AT+J/OzM+X5oenLtuUvwdK+wRwfGP1NcC+NvuXJM1Nm0/vBHAlsCcz/66xaCewqUxvAq5vtL+zfIrndODxqctAkqTRWNZi3dcAvw3cHRF3lrY/BS4FrouIzcCDwPll2Q3A2cA48ARwQYt9S5L60HfoZ+Z/0P06PcD6Lv0TuLDf/UmS2vMbuZJUEUNfkipi6EtSRQx9SaqIoS9JFTH0Jakihr4kVcTQl6SKGPqSVBFDX5IqYuhLUkUMfUmqiKEvSRUx9CWpIoa+JFXE0Jekihj6klQRQ1+SKmLoS1JFDH1JqoihL0kVMfQlqSKGviRVxNCXpIoY+pJUEUNfkipi6EtSRQx9SaqIoS9JFTH0JakiIw/9iDgzIu6LiPGI2Drq/UtSzUYa+hFxGPAPwFnAycDbI+LkUdYgSTUb9Zn+acB4Zj6QmT8BPgNsGHENklStZSPe32rgocb8BPCqZoeI2AJsKbM/ioj7WuzvGOAHLdYfFuuaG+uaG+uamwVZV1zWqq4Xz7Rg1KEfXdryWTOZ24BtA9lZxO7MHBvEtgbJuubGuubGuuamtrpGfXlnAji+Mb8G2DfiGiSpWqMO/a8B6yLixIg4HNgI7BxxDZJUrZFe3snMpyPivcCNwGHAVZl5zxB3OZDLRENgXXNjXXNjXXNTVV2RmbP3kiQtCX4jV5IqYuhLUkUWdehHxPkRcU9E/CwiZvxo00w//VDeUL49IvZGxLXlzeVB1HVUROwq290VESu69HlDRNzZ+PfjiDi3LPtURHynseyVo6qr9PtpY987G+3zOV6vjIivlvv7roj4zcaygY7XbD8VEhFHlNs/XsZjbWPZJaX9voh4c5s6+qjrooi4t4zPTRHx4sayrvfpiOp6V0RMNvb/u41lm8r9vjciNo24rssbNX07Ih5rLBvmeF0VEQci4lszLI+I+Gip+66IOLWxrP14Zeai/Qe8HHgpcCswNkOfw4D7gZOAw4FvAieXZdcBG8v0x4H3DKiuvwG2lumtwGWz9D8KeAR4bpn/FHDeEMarp7qAH83QPm/jBfwysK5MvwjYDywf9Hgd6nhp9Pk94ONleiNwbZk+ufQ/AjixbOewEdb1hsYx9J6pug51n46orncBH+uy7lHAA+XvijK9YlR1Tev/PjofLBnqeJVtvxY4FfjWDMvPBr5E53tNpwO3D3K8FvWZfmbuyczZvrHb9acfIiKAM4Adpd924NwBlbahbK/X7Z4HfCkznxjQ/mcy17p+br7HKzO/nZl7y/Q+4ACwckD7b+rlp0Ka9e4A1pfx2QB8JjOfzMzvAONleyOpKzNvaRxDt9H5HsywtflplTcDuzLzkcx8FNgFnDlPdb0duGZA+z6kzPwKnZO8mWwArs6O24DlEbGKAY3Xog79HnX76YfVwNHAY5n59LT2QTguM/cDlL/HztJ/IwcfcH9ZXtpdHhFHjLiuIyNid0TcNnXJiQU0XhFxGp2zt/sbzYMar5mOl659yng8Tmd8ell3mHU1baZztjil2306yrp+o9w/OyJi6guaC2K8ymWwE4GbG83DGq9ezFT7QMZr1D/DMGcR8WXgl7os+kBmXt/LJrq05SHaW9fV6zbKdlYBv0LnuwtTLgG+TyfYtgEXAx8eYV0nZOa+iDgJuDki7gZ+2KXffI3XPwObMvNnpbnv8eq2iy5t02/nUI6pWfS87Yh4BzAGvK7RfNB9mpn3d1t/CHX9K3BNZj4ZEe+m8yrpjB7XHWZdUzYCOzLzp422YY1XL4Z6fC340M/MN7bcxEw//fADOi+blpWztTn9JMSh6oqIhyNiVWbuLyF14BCbehvwhcx8qrHt/WXyyYj4JPDHo6yrXD4hMx+IiFuBU4DPMc/jFREvAL4IfLC87J3adt/j1UUvPxUy1WciIpYBL6Tzcn2YPzPS07Yj4o10nkhfl5lPTrXPcJ8OIsRmrSsz/6cx+0/AZY11Xz9t3VsHUFNPdTVsBC5sNgxxvHoxU+0DGa8aLu90/emH7Lwzcgud6+kAm4BeXjn0YmfZXi/bPehaYgm+qevo5wJd3+UfRl0RsWLq8khEHAO8Brh3vser3HdfoHOt87PTlg1yvHr5qZBmvecBN5fx2QlsjM6ne04E1gH/2aKWOdUVEacA/wick5kHGu1d79MR1rWqMXsOsKdM3wi8qdS3AngTz37FO9S6Sm0vpfOm6FcbbcMcr17sBN5ZPsVzOvB4ObEZzHgN6x3qUfwDfp3Os9+TwMPAjaX9RcANjX5nA9+m80z9gUb7SXQelOPAZ4EjBlTX0cBNwN7y96jSPgZ8otFvLfA94DnT1r8ZuJtOeP0L8PxR1QW8uuz7m+Xv5oUwXsA7gKeAOxv/XjmM8ep2vNC5XHROmT6y3P7xMh4nNdb9QFnvPuCsAR/vs9X15fI4mBqfnbPdpyOq66+Be8r+bwFe1lj3d8o4jgMXjLKuMv9nwKXT1hv2eF1D59NnT9HJr83Au4F3l+VB5z+bur/sf6yxbuvx8mcYJKkiNVzekSQVhr4kVcTQl6SKGPqSVBFDX5IqYuhLUkUMfUmqyP8DGSgy3uwtFv0AAAAASUVORK5CYII=\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,bins=100);"
   ]
  },
  {
   "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}