{
 "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 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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "zfit.run.check_numerics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING:tensorflow:From /Users/davide/miniconda3/envs/zfit_env/lib/python3.6/site-packages/tensorflow/python/ops/resource_variable_ops.py:435: colocate_with (from tensorflow.python.framework.ops) is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "Colocations handled automatically by placer.\n"
     ]
    }
   ],
   "source": [
    "mmu_mass = 105.0\n",
    "mB_mass = 5280.0\n",
    "mKst_mass = 892.0\n",
    "\n",
    "minq=2*mmu_mass\n",
    "maxq=(mB_mass-mKst_mass)\n",
    "\n",
    "el1 = Particle('l1', mmu_mass)\n",
    "el2 = Particle('l2', mmu_mass)\n",
    "\n",
    "def modq(min_mass, max_mass, n_events):\n",
    "    \n",
    "        min_mass = tf.cast(min_mass, tf.float64)\n",
    "        max_mass = tf.cast(max_mass, tf.float64)\n",
    "\n",
    "        min_mass = tf.broadcast_to(min_mass, shape=(n_events,))\n",
    "        \n",
    "        modq_mass = tfp.distributions.Uniform(low=min_mass, high=max_mass).sample()\n",
    "        \n",
    "        return modq_mass\n",
    "    \n",
    "q=Particle('q', modq).set_children(el1, el2)\n",
    "Kst = Particle('Kst', mKst_mass)\n",
    "\n",
    "B = Particle('B', mB_mass).set_children(Kst, q)\n",
    "\n",
    "#weights, particles = B.generate(n_events=10000)\n",
    "weights_np, particles_np = B.generate(n_events=1000000)\n",
    "weights, particles = B.generate_tensor(n_events=1000000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "105.00049473122846 5279.995436140764\n"
     ]
    }
   ],
   "source": [
    "Kst=particles_np['Kst']\n",
    "l1 =particles_np['l1']\n",
    "l2 =particles_np['l2']\n",
    "q_np=particles_np['q']\n",
    "q2_array=scalar_product_4_np(q_np,q_np)\n",
    "q2_range=np.arange(q2_array.min(), q2_array.max(),1000)\n",
    "print(np.sqrt(q2_array.min()/4), np.sqrt(q2_array.max())+mKst_mass)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3EAAAE6CAYAAAC4dNvKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd5iV1bn38e89haGjVKmCgAgWLIgNe8MG0WiCGqNRgz0aU16SY0yOJznGFHsLicYaS6wYiMYOKipgQ0QFFAVFQVGwIQLr/WM2nnEcZICZeeaZ+X6uay72fp619/5NNA73rLXuFSklJEmSJEn5UJR1AEmSJElS9VnESZIkSVKOWMRJkiRJUo5YxEmSJElSjljESZIkSVKOWMRJkiRJUo5kWsRFxNURMT8iXqzG2Asi4rnC16sR8WFdZJQkSZKk+iSyPCcuInYBPgauSylttgavOw3YKqV0bK2FkyRJkqR6KNOZuJTSeGBhxWsR0Tsi7o2IKRExISI2qeKlhwM31UlISZIkSapHSrIOUIXRwIkppRkRsR1wObDHypsRsSHQC3goo3ySJEmSlJl6VcRFREtgR+CfEbHyclmlYSOA21JKy+symyRJkiTVB/WqiKN8eeeHKaUtv2HMCOCUOsojSZIkSfVKvTpiIKW0GHg9Ig4DiHIDV96PiH7A+sDEjCJKkiRJUqayPmLgJsoLsn4RMTcijgOOBI6LiOeBacDwCi85HLg5ZdlSU5IkSZIylOkRA5IkSZKkNVOvllNKkiRJkr6ZRZwkSZIk5Uhm3Snbt2+fevbsmdXHS5Lq0JQpU95LKXXIOkde+DNSkhqHtf35mFkR17NnTyZPnpzVx0uS6lBEvJF1hjzxZ6QkNQ5r+/PR5ZSSJEmSlCMWcZIkSZKUIxZxkiRJkpQjFnGSJEmSlCMWcZIkSZKUIxZxkiRJkpQjFnGSJEmSlCMWcZIkSZKUIxZxkiRJkpQjJVkHqC3vLFrC83M/5KMly2hWWky39ZvRb4NWNC0tzjqaJEmSJK21BlfEzV+8hF/d/SL/eeldUvrqvSbFRWzZfT1236QjB2zemR7tmmcTUpIkSZLWUoMq4t5ZtIRvX/EE73/yOafu3oc9NulI+5ZlfLRkGW8u/IRn3vyQJ2a9x3n3vsx5977MwG5tGDG4B8MGdqFFWYP6n0KSpGrpOWrsV57P/v0B1bonScpOg6lcUkr85J/P8eGnS7ntxB3ZrGubr9wf0KU1QzfrDMCchZ8ybuo87nz2LX5xx1R+N3Y6B2/VlaN33JA+HVtlEV+SJEmSqqXBFHGPvrqAx2e+zznDN/1aAVdZ97bNOWHX3ozcZSOeefMDbnzyTW6ZPIfrn3yDvfp35IRdezNow/WJiDpKL0lS/VdxZs5ZOUnKToMp4m548g3atyzj8ME9qv2aiGCbDduyzYZtOevAAVw3cTbXPjGbw66cyFY91uOEXXqzz4BOFBVZzEmSJEmqHxpEEffRki94+JUFHD+kF6XFa3dqQtsWTThjr405YZfe3DZlDqMnvMaJN0xh404t+dGefdl/s84Wc5KkBq/yPri1eZ2zdJJUu1ZbxEVEd+A6YANgBTA6pXRRpTG7AXcDrxcu3ZFSOqdmo67axFnvs3xFYrd+Hdf5vZo1KeaoHXpy+OAejJ06j4sfnMGp/3iWjTvNsJiTJDUIa1uo1fR7SJLWTnVm4pYBP0kpPRMRrYApEXF/SumlSuMmpJQOrPmIq/fU6wtpWlrE1huuV2PvWVJcxPAtu3LgFl0s5iRJkiTVG6st4lJK84B5hccfRcR0oCtQuYjLzPR5i+m3QWvKSmr+IO/iomDYwC4csHnnrxVzZ+7dj3037WQDFEmSKvBoAkmqXWu0Jy4iegJbAU9VcXuHiHgeeBv4aUppWhWvHwmMBOjRo/oNSL5JSomX3/mIvft3qpH3W5XKxdyF97/KiTdMYWD39fj5vv3YqU/7Wv18SVJ2ImIocBFQDPwtpfT7SvfLKN96sA3wPvDdlNLsiCgF/gZsTfnP3OtSSufWaXhc+ihJDU21u4BEREvgduCMlNLiSrefATZMKQ0ELgHuquo9UkqjU0qDUkqDOnTosLaZv+K9j5ey8JOl9Nugbs53W1nM/efHu3Detzdn/uIlHPm3pzjyb0/y3JwP6ySDJKnuREQxcBmwHzAAODwiBlQadhzwQUqpD3ABcF7h+mFAWUppc8oLvBMKvxCVJGmtVWsmrvCbxNuBG1NKd1S+X7GoSymNi4jLI6J9Sum9motatbkffApAz/bNa/ujvqKkuIjvbtuD4Vt25can3uSyh2fyrcseZ99NO/HTffrRt5OHhktSAzEYmJlSeg0gIm4GhvPVbQXDgd8UHt8GXBrla+0T0CIiSoBmwFKg8i9CG7zqdq60w6UkVc9qZ+IKP4SuAqanlM5fxZgNCuOIiMGF932/JoOuyrxFSwDYoHWzuvi4r2laWsxxQ3ox/ue78+O9Nubxme+z74XjOfPW55iz8NNMMkmSalRXYE6F53ML16ock1JaBiwC2lFe0H1C+d7yN4E/pZQW1nZgSVLDVp2ZuJ2Ao4CpEfFc4dovgR4AKaUrgUOBkyJiGfAZMCKllGoh79esLOI6t2laFx+3Si3LSjh9r74ctcOGXPHITK6d+Ab3PP8239t+Q07boy9tWzTJNJ8kaa1V1b2q8s+4VY0ZDCwHugDrAxMi4oGVs3pfeYNa2DcuSWqYqtOd8jGq/uFUccylwKU1FWpNvLPoM8pKiliveWkWH/81bVs04b8OGMCxQ3px8YMzuPaJ2dw2eS4n7tabY3fqRbMmNd9BU5JUq+YC3Ss870Z5E6+qxswtLJ1sAywEjgDuTSl9AcyPiMeBQcDXiriU0mhgNMCgQYPq5BehkqR8WqPulPXRO4s/Z4M2Tetdm//ObZpx7iFbcOxOvTjv3lf4432vcP3ENzhzn4359tbdKPaMOUnKi0lA34joBbwFjKC8OKtoDHA0MJHy1SkPpZRSRLwJ7BERNwDNge2BC+sseT1kp0xJWnfV7k5ZX3346VLWa15/lyr27dSKvx09iFtGbk+nNk35+W0vsP9FE3j4lfnU0YpTSdI6KOxxOxW4D5gO3JpSmhYR50TEsMKwq4B2ETETOBMYVbh+GdASeJHyYvDvKaUX6vQbkCQ1OLmfiVv82Re0qcdF3ErbbdSOu07ekXFT3+EP973MD/4+iR17t+MX+/Vn825tso4nSfoGKaVxwLhK186u8HgJ5ccJVH7dx1VdlyRpXeR+Jm7RZ1/Qpln92A+3OhHBAVt05v4f78p/D9uUl9/5iIMufYwf3fSsnSwlSZIkVUsDKeLyNaHYpKSIo3fsyaM/241Td+/Df156hz3+/Ajn3PMSH3yyNOt4kiRJkuqxXBdxKSUWL1mWm5m4ylo1LeWn+/bjkZ/uzre37sY1T7zOLn98mCsemcWSL5ZnHU+SJElSPZSvKaxKPlm6nOUrUm6LuJU2aNOU3397C44d0ovz/v0y5937Mjc8+QY/H9qPg7boQpGdLCVJjcw3dbGc/fsD6jCJJNU/uZ6JW/zZFwC0bprvIm6ljTu14qpjtuUfP9yO9ZqXcvrNz3HwFU8wefbCrKNJkiRJqidyXcR9urR8yWHzslxPKH7Njr3bc8+pQ/jTYQN5Z9FnHHrlRE66YQpvvP9J1tEkSZIkZSzX1c/KfWNNS3Jdi1apqCg4dJtu7L/5Bvx1/Ov8ZfwsHpj+Lkfv0JPT9uhLm+YNY/ZRkqQ15VJLSY1drqufL4u40uKMk9Se5k1KOH2vvjzy0904ZKtuXPX46+z6p4e5+rHXWbpsRdbxJEmSJNWxnBdx5UVMQy7iVurYuinnHboFY0/bmc26tOGcf73EvheO575p75BSyjqeJEmSpDqS8yJu5Uxcrr+NNTKgS2uuP24wfz9mW4qLghOun8KI0U8yde6irKNJkiRJqgO5rn6WLGv4yymrEhHsvklH7j19Z/7nW5sxY/7HHHTpY5x5y3PMW/RZ1vEkSZIk1aKcNzYpX07ZrJEVcSuVFBdx1PYbMnzLLlz+8Cyufvx1xr04jx/uvBEn7tqbFg2sa6ckSatTuemJjU4kNUS5/lv+Z4XllGWNaDllVVo3LWXUfptw5HY9+ON9r3DJQzO5edIcfrL3xhw2qDvFHhYuSWqk7GQpqSHKdfXzeSPoTrkmurdtzsWHb8WdJ+9Ij7bNGXXHVA64eAITZizIOpokSZKkGpLrIu7/zomziKtoqx7rc9uJO3D5kVvzydJlHHXV0xx3zSRmLfg462iSJEmS1lHOi7gVFAWUFrtcsLKIYP/NO/PAmbvyi/024anXF7LvBeM5556XWPTpF1nHkyRJkrSWcr0nbunyFTQpKSLCIm5VykqKOWHX3hyydTfOv/9V/v7E69z57Fx+vPfGHDG4ByXFua7jJUmqdTZLkVTf5LqI+2L5CkqLLEKqo0OrMs49ZHOO2n5D/udfL3H23dO4fuIb/OrAAeyycYes40mSVOe+qemJJNVnua6Alq9IlLiUco0M6NKaf/xwO0YftQ1Ll6/g+1c/zbHXTGLmfPfLSZIkSXmQ6yLui+WJYmfi1lhEsM+mG/CfH+/CL/ffhEmvL2ToheP573um8eGnS7OOJ0mSJOkb5LoCWr5ihU1N1kFZSTEjd+nNwz/bje9s251rn5jNbn96hGufmM0Xy1dkHU+SJElSFXJdxC1bnjzIuga0b1nG/x68OWN/tDObdmnNr8dMY7+LJvDIK/OzjiZJkiSpknwXcSsSpXZXrDH9O7fmhuO246/fH8Sy5Ss45u+TOObvTzNz/kdZR5OkTEXE0Ih4JSJmRsSoKu6XRcQthftPRUTPwvUjI+K5Cl8rImLLus4vSWpYct2dctmKFc7E1bCIYO8Bndh14w5cN3E2Fz04g30vnMBR22/IGXv1Zb3mTbKOKEl1KiKKgcuAvYG5wKSIGJNSeqnCsOOAD1JKfSJiBHAe8N2U0o3AjYX32Ry4O6X0XN1+B6ppFbtaetyApCzkehpr2fJEiUVcrWhSUsTxO2/EIz/djRHbdue6ibPZ9Y+PcM3jr7tfTlJjMxiYmVJ6LaW0FLgZGF5pzHDg2sLj24A94+uHmB4O3FSrSSVJjULOZ+I8YqC2tWtZxu8O3pyjdtiQ3/5rOr+55yWuf/INzjpwALv365h1PEmqC12BORWezwW2W9WYlNKyiFgEtAPeqzDmu3y9+FMD4qHgkupKvmfiViRKPGKgTmyyQWuuP24wVx09iBUJfvD3SRx99dPMeNf9cpIavKp+W5jWZExEbAd8mlJ6cZUfEjEyIiZHxOQFCxasXVJJUqOQ6wpo2fIVLqesQxHBnv07cd8Zu3DWAf155s0PGHrRBH5994t88Inny0lqsOYC3Ss87wa8vaoxEVECtAEWVrg/gtUspUwpjU4pDUopDerQocM6h5YkNVz5LuJcTpmJlfvlHv3Z7hwxuAfXP/kGu//5Ea6bOJtl7peT1PBMAvpGRK+IaEJ5QTam0pgxwNGFx4cCD6WUEkBEFAGHUb6XTpKkdZbvPXHLV9C8Sa6/hVxr26IJ//OtzThy+x6cc89LnH33NG588k1+fdAAduzTPut4klQjCnvcTgXuA4qBq1NK0yLiHGBySmkMcBVwfUTMpHwGbkSFt9gFmJtSeq2us6v2Vd4HJ0l1IdcV0HJn4uqFTTZozY3Hb8d9097ht2Onc8TfnmLophvwXwf0p3vb5lnHk6R1llIaB4yrdO3sCo+XUD7bVtVrHwG2r818kqTGJddF3BceMVBvRARDN+vMbv068rcJr3HZw7N46JX5jNx5I07evbczppKkRuebzpPzrDlJ6yLXf7NeviJ52Hc907S0mFP36Muh23Tn9/+ezqUPz+S2KXMZtd8mDN+yC18/NkmSVBtc5le/+M9DUk3KeWOTFZQU5/pbaLA2aNOUC0dsxe0n7UDH1mWccctzHHrlRF6Y+2HW0SRJkqRcy3UFVH5OnDM79dk2G7blrpN34g+HbsEb73/K8Mse5+e3Pc/8j5ZkHU2SJEnKpdUup4yI7sB1wAbACmB0SumiSmMCuAjYH/gUOCal9EzNx/2qZctdTpkHRUXBdwZ1Z7/NNuCSh2by98dfZ9zUd/jRnn04ZsdeNCnJ9e8SJElaJ5WXWrpHTtLqVOdvz8uAn6SU+lPeXeuUiBhQacx+QN/C10jgihpNuQorkjNxedKqaSm/3L8/952xC4N7teV/x73MvheO56GX3806miRJkpQbqy3iUkrzVs6qpZQ+AqYDXSsNGw5cl8o9CawXEZ1rPG0lK1IisIjLm406tOTqY7bl7z/Ylgg49prJHPP3p5k5/+Oso0mSJEn13hqtY4uInsBWwFOVbnUF5lR4PpevF3o1LiUociVebu3eryP3nr4LZx3QnymzP2DoheP57b9eYvGSL7KOJkmSJNVb1S6BIqIlcDtwRkppceXbVbwkVfEeIyNickRMXrBgwZolrcKKhC3rc65JSRHH77wRD/9sNw4b1I2rHn+d3f/4CDc//SbLV3ztXyFJkiSp0atWERcRpZQXcDemlO6oYshcoHuF592AtysPSimNTikNSikN6tChw9rkrfx+LqZsINq3LOPcQ7bgnlOHsFGHFoy6YyrDL3uMSbMXZh1NkqRc6Dlq7Jdfkhq26nSnDOAqYHpK6fxVDBsDnBoRNwPbAYtSSvNqLmbVElDkTFyDslnXNtx6wg7c88I8zh03ncOunMiwgV0Ytd8mdFmvWdbxJEmqdd9UhNm5UhJUo4gDdgKOAqZGxHOFa78EegCklK4ExlF+vMBMyo8Y+EHNR/26FSlhc8qGJyIYNrALe/XvyJWPvsZfHp3F/S+9y0m79WbkLhvRtLQ464iSJElSZlZbxKWUHqPqPW8VxyTglJoKVV0rViT3xDVgzZuUcObeG3PYNt0499/TOf/+V7l18hzOPnAAew/o5D97SVKjU3GWzlk5qfGqzkxcvZUS+Pf4hq972+ZcfuQ2PDHrPX4zZhojr5/CLht34NcHDaB3h5ZZx5MkKRPufZMar1w36HdPXOOyY+/2jP3Rzpx94ACefbP8SIJzx03n48+XZR1NkiRJqjO5nolzT1zjU1pcxLFDejFsyy784d6X+cv417jz2bf4xf6b8K0tu7rEUpKkSirP2LkMU8q/XM/ErUjuiWus2rcs4w+HDuTOk3ekc5um/PiW5znsyolMe3tR1tEkSZKkWpXrIs49cdqqx/rcefJOnPftzXn9vU846JLHOOuuqXzwydKso0mSJEm1IvdFnHviVFQUfHfbHjz00934/g49uenpOez+50e44ck3WL4iZR1PkiRJqlG5LuLcE6eK2jQr5TfDNmXsj4bQr1MrzrrrRYZd+hiTZy/MOpokSfVGz1Fjv/ySlE+5b2wS33yEnRqhTTZozc0jt+dfL8zjf8dN59ArJ3LIVl0Ztd8mdGzdNOt4kiTVOgs0qWHL9Uxc+REDWadQfRQRHDSwCw/+ZFdO2b03/3phHnv8+VFGj5/F0mUrso4nSZIkrbXcFnEppUJjE6s4rVrzJiX8bN9N+M+Pd2Fwr7b877iXGXrReMa/uiDraJJyJCKGRsQrETEzIkZVcb8sIm4p3H8qInpWuLdFREyMiGkRMTUiXBIgSVonOS7iyv+0hlN19GzfgquP2ZarjxnE8hWJ71/9NCdcP5k5Cz/NOpqkei4iioHLgP2AAcDhETGg0rDjgA9SSn2AC4DzCq8tAW4ATkwpbQrsBnxRR9ElSQ1UbvfErew5aHdKrYk9NunEjr3bc9Vjr3PpQzPZ65VHOXHX3py0W2+alhZnHU9S/TQYmJlSeg0gIm4GhgMvVRgzHPhN4fFtwKVRvlRkH+CFlNLzACml9+sqtFQd33QQuIeES/VXbmfiVhSm4twTpzXVtLSYU3bvw4M/2ZW9B3TiogdnsOefH+XeF+eRkkcSSPqarsCcCs/nFq5VOSaltAxYBLQDNgZSRNwXEc9ExM/rIK8kqYHLfRHnnjitrS7rNePSI7bmph9uT8uyEk684Rm+f/XTzFrwcdbRJNUvVf2gqfwbn1WNKQGGAEcW/jw4Ivas8kMiRkbE5IiYvGCB+3YlSauW2yLOPXGqKTv0bsfYHw3h1wcN4Lk5HzL0wvH84d6X+XTpsqyjSaof5gLdKzzvBry9qjGFfXBtgIWF64+mlN5LKX0KjAO2rupDUkqjU0qDUkqDOnToUMPfgiSpIcnvnrhCEeeeONWEkuIifrBTLw7cogvn/ns6lz8yi7uefYuzDxrAvptu4Iyv1LhNAvpGRC/gLWAEcESlMWOAo4GJwKHAQymlFBH3AT+PiObAUmBXyhufSPWS58tJ+ZDbmTj3xKk2dGhVxvnf2ZJbT9iB1s1KOfGGZzj675N4/b1Pso4mKSOFPW6nAvcB04FbU0rTIuKciBhWGHYV0C4iZgJnAqMKr/0AOJ/yQvA54JmUkn9LliStk9zOxH25J67KbQjSuhncqy3/Om0I1018g/Pvf5V9LxjPCbtuxMm79aFZE7tYSo1NSmkc5UshK147u8LjJcBhq3jtDZQfMyDlWsVZOjtVStnKbRG3cke5q9xUW0qKizh2SC8O3KIz5/77ZS55aCZ3PPMWvz5oAHsP6OQSS0lSo+XxA1K2crucMq0o/9M9captHVs35YLvbsnNI7enRVkxI6+fwrHXTOKN911iKUmSpLqX2yLu/44YyDiIGo3tN2rH2B/tzFkH9Ofp1xey9wXjOf/+V1nyxfKso0mSJKkRyf1ySmfiVJdKi4s4fueNOGhgF343djoXPziDO5+dy28O2pQ9+3fKOp4kSZlwv5xUt3I/E2d3SmWhU+umXHz4Vvzjh9tRVlLMcddO5vhrJ/Hm+59mHU2SJEkNXO6LOJtLKEs79m7Pv0/fmV/uvwlPzHqfvS94lIsemOESS0mSJNWa/C6nLKyntIZT1kqLixi5S+8vl1he8MCr3P7MXP572KbsvknHrONJklSn7Fwp1b7czsStLOLcE6f6onObZlx6xNbcePx2lBYHP7hmEj+8bjJzFrrEUpLUePUcNfbLL0k1I7dFnHviVF/t1Kc9/z59F/7f0E14bMZ77HX+o1zy4Aw+X+YSS0mSJK273BdxgVWc6p8mJUWctFtvHvzJruzZvyN/vv9V9rtwAo/PfC/raJIkScq53BZx7olTHnRZrxmXH7kN1x07mOUpceTfnuJHNz3L/I+WZB1NkiRJOdUAGptYxan+22XjDtx3xi5c8cgsrnhkFg+/PJ+fDe3HkdttSLFrgiVJjYRNT6SakduZuJX866/yomlpMT/ee2PuPWNnBnZfj7PvnsbBlz/OC3M/zDqaJEmSciT3RZyUNxt1aMn1xw3m4sO3Yt6iJQy/7HHOvvtFFi/5IutokiTVKTtXSmsnt8sppTyLCIYN7MJu/Trw5/te4bon3+DfL77DWQf0Z9jALi4TliQ1OhULOZdZSt8stzNxiZR1BGmdtW5ayn8P34y7T9mJDVo35fSbn+N7Vz3Faws+zjqaJEmS6qncz8Q5YaGGYItu63HXKTtx41Nv8Md7X2HohRM4cbfenLxbb5qWFmcdT5KkOmUDFOmb5b6IkxqK4qLg+zv0ZOhmG/C7sdO5+MEZ3P3cW5wzfDN23bhD1vEkScqMSy2lr8rtckqpoerYqikXjdiKG4/fjuIIjr76aU75xzO8u9iz5SRJkpTjIi65JU4N3E592vPvM3bmzL035v6X3mXPPz/K1Y+9zrLlK7KOJkmSpAzlfjmle+LUkJWVFPOjPfsyfMsunH33NM7510vc/sxcfvutzdiqx/pZx5Mkqc65X06qxkxcRFwdEfMj4sVV3N8tIhZFxHOFr7NrPqbUuG3YrgXX/GBbLjtia977+HMOueIJzrprKos+82w5SZKkxqY6yymvAYauZsyElNKWha9z1j2WpMoiggO26MwDZ+7KMTv25B9Pvcle5z/K2BfmkVxfLNWqiBgaEa9ExMyIGFXF/bKIuKVw/6mI6Fm43jMiPqvwi84r6zq71NB5YLgao9UWcSml8cDCOsgiqRpaNS3l1wdtyl2n7ETHVmWc8o9nOO7aycz94NOso0kNUkQUA5cB+wEDgMMjYkClYccBH6SU+gAXAOdVuDerwi86T6yT0JKkBq2mGpvsEBHPR8S/I2LTGnrPb+S8gxq7Lbqtx92n7MRZB/Rn4qz32eeC8fxtwms2PpFq3mBgZkrptZTSUuBmYHilMcOBawuPbwP2jHDXtiSpdtREEfcMsGFKaSBwCXDXqgZGxMiImBwRkxcsWFADHw2BPyPVeJUUF3H8zhtx/5m7sP1G7fjt2Ol86/LHmTp3UdbRpIakKzCnwvO5hWtVjkkpLQMWAe0K93pFxLMR8WhE7FzbYSVJDd86F3EppcUppY8Lj8cBpRHRfhVjR6eUBqWUBnXo4OHFUk3ptn5zrjp6EJcdsTXvLv6c4Zc9xjn3vMQnny/LOprUEFT128LKC0JWNWYe0COltBVwJvCPiGhd5YfUwi86JUkN0zofMRARGwDvppRSRAymvDB8f52TSVojKxufDOnbnj/c+zJXP/469744j3OGb8ZeAzplHU/Ks7lA9wrPuwFvr2LM3IgoAdoAC1N516HPAVJKUyJiFrAxMLnyh6SURgOjAQYNGuSuAWktePyAGovVFnERcROwG9A+IuYCvwZKAVJKVwKHAidFxDLgM2BEqoNWeXbjk6rWplkpvzt4cw7Zuiu/uGMqx183mf0224DfDNuUTq2bZh1PyqNJQN+I6AW8BYwAjqg0ZgxwNDCR8p+LDxV+udmB8mJueURsBPQFXqu76FLjVrGos6BTQ7LaIi6ldPhq7l8KXFpjidaQ28alqm2zYVv+ddrO/HXCa1z84AwmzHiPnw/tx5HbbUhxkf/HkaorpbQsIk4F7gOKgatTStMi4hxgckppDHAVcH1EzKS8o/OIwst3Ac4p/KJzOXBiSsmOz1I9YIGnPFvn5ZSS6q8mJUWcsnsfDti8M2fd9SJn3z2NO555i3MP2Zz+navcliOpCoU93+MqXTu7wuMlwGFVvO524PZaDyhJalQs4qRGoGf7Flx/3E32hJMAACAASURBVGDueu4t/udf0znoksc4fueNOH3PvjRrUpx1PEmSap2Hgashqalz4uqcO+KkNRMRHLxVNx48c1cO3qorVz46i30vHM+EGXbBkyRJypPcFnGS1s76LZrwx8MGctMPt6ekKDjqqqf5ya3P88EnS7OOJkmSpGpwOaXUSO3Qux3jTt+ZSx+ayZWPzuLRV+fz64M25cAtOhN2DJIkNSIeTaC8cSZOasSalhbz0337MebUIXRZrxmn3fQsP7xuMvMWfZZ1NEmSJK2CRZwkBnRpzR0n7ch/7d+fx2a+xz7nj+fGp95gxQp3n0qSGp+eo8Z++SXVR7ldTulZ31LNKiku4oe7bMQ+m3biF3dM5b/ufJG7n3ub3x+yORt1aJl1PEmSMuF5cqqPcj8T594dqWZt2K4FNx6/Hed9e3Omz1vM0IsmcPkjM/li+Yqso0mSJIkGUMRJqnkRwXe37cGDZ+7Knpt05A/3vsLwSx9n6txFWUeTJElq9CziJK1Sx9ZNueJ723Dl97Zhwcef863LH+fccdP5bOnyrKNJkiQ1Wjku4twUJ9WVoZttwANn7sph23TjL+NfY+hF43li1ntZx5IkSWqUctvYZCV3xEl1o02zUn7/7S0YNrALo+6YyhF/fYrDB3dn1H79adOsNOt4kiTVum/qVmnTE9WlHM/EScrCjn3ac98ZuzByl424ZdIc9j7/UR546d2sY0mSJDUaFnGS1lizJsX8cv/+3HXKTqzfvAnHXzeZM295jg8/XZp1NEmSpAbPIk7SWtui23rcc9oQfrRHH8Y8/zZ7XzCe/0x7J+tYkiRJDVpu98R52LdUPzQpKeLMffqxz6Yb8LPbXmDk9VMYNrALvxm2KW1bNMk6niRJdaLyfjn3yKk25X4mzrO+pfphs65tuPuUnfjxXhszbuo89rngUe59cV7WsSRJkhqc3BdxkuqPJiVFnL5XX8acOoROrZty4g3PcOo/nuH9jz/POpokSVKDkdvllJLqrwFdWnPXKTtx5SOzuPihGUyc9T7nDN+MA7bonHU0SZLqRMXllS6tVE3L7UycW+Kk+q20uIjT9uzLv07bmS7rNeOUfzzDSTdMYcFHzspJkiSti9zPxIXHfUv1Wr8NWnHnyTsyesJrXHj/DJ587VF+M2xThg3sQripVZLUCNj0RDUt90WcpPqvpLiIk3frw979O/HT217g9JufY+wL8/jtwZvRsVXTrONJklSnXGqpdZXb5ZSS8qdvp1bcfuIOjNpvEx55dQH7XjCecVPtYKn6LyKGRsQrETEzIkZVcb8sIm4p3H8qInpWut8jIj6OiJ/WVWZJUsOV2yLOc+KkfCopLuLEXXsz9rQhdFu/OSff+Ayn3/wsiz79IutoUpUiohi4DNgPGAAcHhEDKg07DvggpdQHuAA4r9L9C4B/13ZWSVLjkNsibiW31Ej51LdTK+44eUfO2KsvY1+Yxz4XPsojr8zPOpZUlcHAzJTSaymlpcDNwPBKY4YD1xYe3wbsGYVNnxHxLeA1YFod5ZUkNXC5L+Ik5VdpcRFn7LUxd568E62blnLM3yfxyzun8snny7KOJlXUFZhT4fncwrUqx6SUlgGLgHYR0QL4f8B/r+5DImJkREyOiMkLFiyokeCSpIbJxiaSMrd5tzbcc9oQzr//Vf464TUmzFjAnw/bksG92mYdTQKqbINceVH/qsb8N3BBSunj1XVjTSmNBkYDDBo0yE0DUiNh50qtDYs4SfVC09Jifrl/f/bq34mf/vN5vjt6IscP6cVP9ulH09LirOOpcZsLdK/wvBvw9irGzI2IEqANsBDYDjg0Iv4ArAesiIglKaVLaz+2pDyyc6WqI7fLKZPHfUsN0uBebfn36TtzxOAe/HXC6xx0yWNMnbso61hq3CYBfSOiV0Q0AUYAYyqNGQMcXXh8KPBQKrdzSqlnSqkncCHwvxZwkqR1ldsibiX7mkgNT4uyEn538OZce+xgFi/5goMvf5wLH3iVL5avyDqaGqHCHrdTgfuA6cCtKaVpEXFORAwrDLuK8j1wM4Ezga8dQyBJUk1xOaWkemvXjTvwnzN25ddjXuTCB2bw4PT5XPDdLenTsWXW0dTIpJTGAeMqXTu7wuMlwGGreY/f1Eo4SQ2W++W0KrmfiZPUsLVpXsqFI7bi8iO3Zu4Hn3LgJRO4/sk3SB4WKUmSGqnczsT59zepcdl/885ss+H6/PSfz/Oru17kkZfnc96hW9C+ZVnW0SRJqhM2PdFKuZ+J87BvqfHo1Lop1/5gML8+aAATZr7H0AvH89DL72YdS5IkqU7lvoiT1LgUFQU/2KkX95w6hPYtyzj2msmcdddUPlu6POtokiRJdSK3yyklNW79NmjF3afuxB/vfYW/PfY6E2e9z0UjtmKzrm2yjiZJUq2z6UnjltuZOPfESSorKeasAwdw4/Hb8cnnyzn48se54pFZLF/hfyAkSVLDldsi7v+4KU5q7Hbq0557z9iZvQd04rx7X+aIvz7JWx9+lnUsSZIy0XPU2C+/1DCttoiLiKsjYn5EvLiK+xERF0fEzIh4ISK2rvmYkvTN1mvehMuO2Jo/HroFL761iKEXjuee59/OOpYkSVKNq86euGuAS4HrVnF/P6Bv4Ws74IrCn5JUpyKCwwZ1Z7te7Tjjlmc57aZneWzGe/x62ACaN3ELsCSp4XLWrXFZ7UxcSmk8sPAbhgwHrkvlngTWi4jONRVQktZUj3bNueWEHThl997cOmUOB13yGNPnLc46liRJUo2oiT1xXYE5FZ7PLVyrVQkbF0hatdLiIn627ybccNx2LF6yjOGXPc71E2eT7IokSZJyriaKuKo6i1T5t6SIGBkRkyNi8oIFC2rgoz3sW9I326lPe/59+s7s2Lsdv7p7GifeMIUPP12adSxJkqS1VhObROYC3Ss87wZU2U0gpTQaGA0waNAgfx0uqU60b1nG1Udvy1WPvc4f7nuZ/S+awMWHb8Wgnm2zjiZJUq3yPLmGqSZm4sYA3y90qdweWJRSmlcD7ytJNaaoKPjhLhtx+0k7UlpSxHdHP8klD87wTDlJUqPi8QMNQ3WOGLgJmAj0i4i5EXFcRJwYEScWhowDXgNmAn8FTq61tBW4rUXS2tii23r867QhHLhFZ/58/6t8729P8e7iJVnHkiRJqrbVLqdMKR2+mvsJOKXGEq0ht8RJWlOtmpZy4Xe3ZEif9px99zT2u2gCF43Ykp37dsg6miRJ0mrVxHJKScqdlWfK3XPaENq3bML3r36aCx941eWVkiSp3vP0W0mNWp+OLbnrlJ04664XufCBGUx54wMu/O6WtGtZlnU0SZJqlU1P8suZOEmNXvMmJfz5sIH8/pDNeer1hRxw8WNMmr0w61iSJElVyn0RFx4UJ6kGRAQjBvfgzpN3pKy0iBGjn2T0+FkeDi5JajTsXJkfuS/iJKkmbdqlDfecNoS9+3fif8e9zMjrp7Dosy+yjiVJkvQlizhJqqR101Ku+N7W/OrAATz88nwOvGQCU+cuyjqWJEkSYBEnSVWKCI4b0otbTtiBZcsT377yCW6dPCfrWMpIRAyNiFciYmZEjKrifllE3FK4/1RE9CxcHxwRzxW+no+Ig+s6uyStjYpLK11eWf/ktohzm4qkurDNhusz9kc7s23P9fn5bS/wq7teZOmyFVnHUh2KiGLgMmA/YABweEQMqDTsOOCDlFIf4ALgvML1F4FBKaUtgaHAXyLCztCSpHWS2yJuJduaSKptbVs04dofDOaEXTbi+iff4Ii/Psn8xUuyjqW6MxiYmVJ6LaW0FLgZGF5pzHDg2sLj24A9IyJSSp+mlJYVrjcF/BWkJGmd5b6Ik6S6UFJcxC/2788lh2/FtLcXc+AljzHljQ+yjqW60RWouJZ2buFalWMKRdsioB1ARGwXEdOAqcCJFYq6r4iIkRExOSImL1iwoIa/BUlaNy6trF8s4iRpDRw0sAt3nLwjTUuLGTF6Ijc+9YbHEDR8VS36qPwPfZVjUkpPpZQ2BbYFfhERTav6kJTS6JTSoJTSoA4dOqxTYElSw5bbIi65IkVSRvp3bs2YU3dix97t+a87X2TU7VNZ8sXyrGOp9swFuld43g14e1VjCnve2gBfOTE+pTQd+ATYrNaSSpIahdxvrvasb0lZWK95E64+ZlvOv/8VLnt4Fi+/+xGjj9qGTq2rnGRRvk0C+kZEL+AtYARwRKUxY4CjgYnAocBDKaVUeM2clNKyiNgQ6AfMrrPkklQLKi+pnP37AzJK0njldiZOkrJWXBT8bN9NuPJ7WzPj3Y8YduljvDD3w6xjqYYV9rCdCtwHTAduTSlNi4hzImJYYdhVQLuImAmcCaw8hmAI8HxEPAfcCZycUnqvbr8DSVJDk/uZOEnK2tDNOrNhuxYcf+1kvvOXifzpsIEcuEWXrGOpBqWUxgHjKl07u8LjJcBhVbzueuD6Wg8oSWpUcjsTZx8BSfVJ/86tueuUndi0SxtO/cezXPjAqzY8kSRJtSL3M3HuiZNUX3RoVcaNx2/HL++YyoUPzGDm/I/502EDaVpanHU0SZJqzTcdO+B+udqR+yJOkuqTpqXF/Pk7A+nbqRV/uO9l3lz4KX/9/iAbnkiSpBqT2+WUklRfRQQn7dab0UcNYub8jxl26WO8+NairGNJkqQGwiJOkmrJ3gM6cftJO1JSVMR3/jKRh1+Zn3UkSZLqVM9RY7/8Us3JbRFnuwBJedC/c2vuPHlHerUv715509NvZh1JkiTlXG6LuJUCO5tIqt86tm7KLSfswJA+7fnFHVP5032v2LlSkiSttdwXcZKUBy3LSvjb0YMYsW13Ln14Jmfe+jxLl63IOpYkScohu1NKUh0pLS7i3EM2p9v6zfjTf17l3cVLuPKobWjdtDTraJIk1bqK++I8emDd5HYmzqVIkvIoIjh1j76c/52BPP36Qg67YiLvLl6SdSxJkpQjuS3ivuSWOEk5dMjW3bj22MHM/eBTvn3FE8x+75OsI0mSpJzIfxEnSTm1U5/23DRyez75fBmHXjmR6fMWZx1JkiTlgEWcJGVoi27r8c8Td6C0OPjOXyYyefbCrCNJklTrKp4f5xlyay63RZw74iQ1FH06tuKfJ+5A+5ZlfO+qp3jEQ8ElSdI3yG0Rt5Jb4iQ1BN3Wb84/T9yBjdq35PhrJ3PP829nHUmSJNVTHjEgSfVE+5Zl3HzC9hx/zWROv/lZvli+gkO27pZ1LEmSap3HD6yZ3M/ESVJD0rppKdceO5gderfjJ/98nn9OnpN1JEmSVM9YxElSPdOsSTFXHb0tQ/q05+e3v8Atk97MOpIkSapHcruc0rO+JTVkTUuL+ev3B3HC9VP4f7dPZfkKOGK7HlnHkiSp1lXuVunyyq/L/UxchK1NJDVMTUuL+ctR27DHJh355Z1TuX7i7KwjSZKkeiD3RZwkNWRNS4u54ntbs1f/Tvzq7mkurZQkSRZxklTflZUUc/mRW7Prxh0YdcdU7n7urawjNToRMTQiXomImRExqor7ZRFxS+H+UxHRs3B974iYEhFTC3/uUdfZJSnvPBT863JcxLkpTlLj0aSkiCu/tw3b9mzLmbc+z/0vvZt1pEYjIoqBy4D9gAHA4RExoNKw44APUkp9gAuA8wrX3wMOSiltDhwNXF83qSVJDVmOi7hy7oiT1Fg0a1LM1cdsy2Zd23DKjc/w2Iz3so7UWAwGZqaUXkspLQVuBoZXGjMcuLbw+DZgz4iIlNKzKaWVJ7dPA5pGRFmdpJYkNVjVKuKqsYzkmIhYEBHPFb6Or/mokqSWZSVc+4Nt2ahDC3543WQmz16YdaTGoCtQ8cC+uYVrVY5JKS0DFgHtKo35NvBsSunzWsopSWokVnvEQIVlJHtT/oNrUkSMSSm9VGnoLSmlU2shoySpgvWaN+H647bjO3+ZyHHXTub2k3agT8dWWcdqyKpa9FF5Tf83jomITSlfYrnPKj8kYiQwEqBHD4+TkKSqePxAuerMxFVnGUmd85w4SY1Zh1ZlXHfsYEqLizj66knMX7wk60gN2Vyge4Xn3YC3VzUmIkqANsDCwvNuwJ3A91NKs1b1ISml0SmlQSmlQR06dKjB+JKkhqY6RVx1lpEAfDsiXoiI2yKiexX3a4XHxElqrLq3bc7fj9mWDz5dyg+umcTHny/LOlJDNQnoGxG9IqIJMAIYU2nMGMoblwAcCjyUUkoRsR4wFvhFSunxOkssSWrQqlPEVWcZyT1Az5TSFsAD/N/m7q++UcTIiJgcEZMXLFiwZkklSV+zebc2XHbk1rz8zkecdMMUvli+IutIDU5hj9upwH3AdODWlNK0iDgnIoYVhl0FtIuImcCZwMr946cCfYBfVdg33rGOvwVJarAa6/ED1SniVruMJKX0foWN2n8FtqnqjVwqIkk1b/d+HTn34M2ZMOM9zrrzRZLrzWtcSmlcSmnjlFLvlNLvCtfOTimNKTxeklI6LKXUJ6U0OKX0WuH6b1NKLVJKW1b4mp/l9yJJyr/qFHGrXUYSEZ0rPB1G+W8qJUl15Dvbdue0Pfpwy+Q5XPPE7KzjSJKkWrTa7pQppWURsXIZSTFw9cplJMDkwm8hf1RYUrKM8o3cx9Ri5vJctf0BkpQzP95rY1555yN+O3Y6fTu2Ykjf9llHkiSpzjSmzpXVOieuGstIfpFS2jSlNDCltHtK6eXaDF1ReNy3JAFQVBSc/90t6d2hBaf84xlmv/dJ1pEkSVItqFYRJ0nKh5ZlJfzt+9sSAT+8brIdKyVJaoAs4iSpgenRrjmXH7E1sxZ8zC/vmGqjE0mSGpjcFnH+nUSSVm3HPu35yT79GPP829z09JzVv0CSJOXGahub1Hce9i1JVTtp19489fpCfnPPNLbsvh4DurTOOpIkSXWmYqOThtbkJLczcZKkb1ZUFFzwnYGs37yUU/7xjPvjJElqICziJKkBa9eyjEsO35o33v+Ec+6ZlnUcSZJUAyziJKmBG9yrLSft1ptbJ8/lgZfezTqOJElaR7ndE2e3NUmqvtP33JgHp89n1B1T+c+G69O2RZOsI0mSVGca2kHguZ+Js6+JJK1ek5Iizv/Oliz6bCm/uvvFrONIkqR1kPsiTpJUPQO6tOaMvTZm7AvzuG/aO1nHkSRJa8kiTpIakZG7bES/Tq347zHT+MRulZIk5VJuizh3xEnSmistLuJ3B2/G24uWcPGDM7KOI0mS1kJuG5t8yU1xkrRGBvVsy3cGdeOqx17nkK270W+DVllHkiSpTuX9IPDczsRJktbeqP3607JpCb8d+1LWUSRJ0hqyiJOkRqhtiyacunsfJsx4jwkzFmQdR5IkrYHcFnEeEydJ6+aoHTak2/rNOHfcy6xY4X9UJUnKi9wWcSuFm+Ikaa2UlRTzs3378dK8xdz9/FtZx5EkSdWU/8YmkqS1dtAWXbjy0de45KGZDBvYleIifzEmSWpcKjY5gXw0Osn9TJwkae0VFQUn79ab1xZ8wn88AFySpFywiJOkRm7/zTvTs11zLn9kFskNx1WKiKER8UpEzIyIUVXcL4uIWwr3n4qInoXr7SLi4Yj4OCIurevckqSGKbdFXPK4b0mqEcVFwYm79mbqW4t4Ytb7WcepdyKiGLgM2A8YABweEQMqDTsO+CCl1Ae4ADivcH0J8Cvgp3UUV5LUCOS2iFsp3L4hSevsW1t1Zf3mpdzw5BtZR6mPBgMzU0qvpZSWAjcDwyuNGQ5cW3h8G7BnRERK6ZOU0mOUF3OSJNWI3BdxkqR117S0mMMGdec/L73Lu4utNyrpCsyp8Hxu4VqVY1JKy4BFQLs6SSdJanQs4iRJABwxuAfLVyRumTRn9YMbl6rWfFRe01+dMd/8IREjI2JyRExesMAD2CUpKz1Hjf3yq77KbxHnljhJqlE927dgx97tuOOZuTY4+aq5QPcKz7sBb69qTESUAG2AhWvyISml0SmlQSmlQR06dFiHuJKkhi6/RVyBW+IkqeYMG9iF2e9/yrS3F2cdpT6ZBPSNiF4R0QQYAYypNGYMcHTh8aHAQ8lKWJJUS3JfxEmSas6+m25ASVEwduq8rKPUG4U9bqcC9wHTgVtTStMi4pyIGFYYdhXQLiJmAmcCXx5DEBGzgfOBYyJibhWdLSVJWiMlWQeQJNUf67dowk592jNu6jz+39BNso5Tb6SUxgHjKl07u8LjJcBhq3htz1oNJ0lqdHI7E9e6WSmDe7WlVdPSrKNIUoOyZ/+OvPH+p7z5/qdZR5EkSVXIbRG3Wdc23HrCDgzo0jrrKJLUoOzUpz0AE2baIVGSpPrI5ZSSpK/YqH0LurRpymMz3uPI7TbMOo4kSZmpeMzA7N8fkGGSr8rtTJwkqXZEBNtt1I4pb3yQdRRJklQFizhJ0tds3rUN8z/6nPmLl2QdRZIkVWIRJ0n6ms27tQFg6luLMk4iSZIqs4iTJH3NgM7lTaOmz/PQb0mS6huLOEnS17QoK6FDqzLeXOgxA5Ik1TcWcZKkKvVo25w5Cz/LOoYkSarEIk6SVKXu6zdjzgfOxEmSVN9YxEmSqtS2RRkffLI06xiSJKmSahVxETE0Il6JiJkRMaqK+2URcUvh/lMR0bOmg0qS6labZqV8snQ5XyxfkXUUSZJUwWqLuIgoBi4D9gMGAIdHxIBKw44DPkgp9QEuAM6r6aCSpLrVplkJAB8tWZZxEkmSVFFJNcYMBmamlF4DiIibgeHASxXGDAd+U3h8G3BpRERKKdVgVklSHWpaWgzA58uWZ5xEkqTs9Rw19ivPZ//+gIySVK+I6wrMqfB8LrDdqsaklJZFxCKgHfBexUERMRIYCdCjR4+1jCxJqgt7D+jEZl3b0K5FWdZRJElSBdXZExdVXKs8w1adMaSURqeUBqWUBnXo0KE6+SRJGWnXsozNurahSYk9sCRJqk+q85N5LtC9wvNuwNurGhMRJUAbYGFNBJQkSZIk/Z/qFHGTgL4R0SsimgAjgDGVxowBji48PhR4yP1wkiRJklTzVrsnrrDH7VTgPqAYuDqlNC0izgEmp5TGAFcB10fETMpn4EbUZmhJkiRJaqyq09iElNI4YFyla2dXeLwEOKxmo0mSJEmSKnO3uiRJkiTliEWcJEmSJOWIRZwkSasREUMj4pWImBkRo6q4XxYRtxTuPxURPSvc+0Xh+isRsW9d5pYkNUwWcZIkfYOIKAYuA/YDBgCHR8SASsOOAz5IKfUBLgDOK7x2AOXNvjYFhgKXF95PkqS1ZhEnSdI3GwzMTCm9llJaCtwMDK80ZjhwbeHxbcCeERGF6zenlD5PKb0OzCy8nyRJa80iTpKkb9YVmFPh+dzCtSrHpJSWAYuAdtV8rSRJa6RaRwzUhilTprwXEW+s49u0B96riTx1yMx1w8x1w8x1oyFk3jCrIDUgqriWqjmmOq8tf4OIkcDIwtOPI+KVaif8uobw70wemLlumLlumHkNxXlr9bIa+fmYWRGXUuqwru8REZNTSoNqIk9dMXPdMHPdMHPdMHPm5gLdKzzvBry9ijFzI6IEaAMsrOZrAUgpjQZG10TgPP7vb+a6Yea6Yea60Zgzu5xSkqRvNgnoGxG9IqIJ5Y1KxlQaMwY4uvD4UOChlFIqXB9R6F7ZC+gLPF1HuSVJDVRmM3GSJOVBSmlZRJwK3AcUA1enlKZFxDnw/9u7/9C76jqO488XzB9Ush+JuUy0iQT2RzpEtB8iGGuO2CoyFkIjgxAV9I+ggTAk8A+L/MOIBHVkIrnSVl9iow0T/GsrG9/Nra3tqwyczg0qtiLEX+/+OJ8bp/M95+62vt/P+Zyvrwdc7rnnfM53r33u537enHvPPZcXI2IKeBx4UtIM1Sdw69O+ByT9Avgz8A5wV0S828t/xMzMFoyhH8TNyWknmTlzHs6chzPn4cw9i4htwLbGuk215TeBWzv2fQB4YF4DzjbE/nfmPJw5D2fO432bWdXZHmZmZmZmZjYE/k6cmZmZmZnZgBR/ECdptaS/SJqRtLFl+3mStqTtuyVdnj/lf+W5VNLzkg5KOiDpnpY2N0k6JWk63Ta1/a3cJB2V9FLK9GLLdkl6OPX1Pkkr+8hZy/OJWh9OSzot6d5Gm977WtJmSScl7a+tWyZpp6Qj6X5px74bUpsjkja0tcmY+QeSDqXnfqukJR37jh1HmTPfL+m12vO/pmPfsfNM5sxbanmPSpru2Dd7P3fNb6WP54XMNTIP18d5y+n62F9m18e5z5y3RkZEsTeqL5C/DKwAzgX2Alc12twJPJKW1wNbes68HFiZli8ADrdkvgn4bd/925L9KHDhmO1rgO1Uv3t0PbC778yNsfIGcFlpfQ3cCKwE9tfWfR/YmJY3Ag+27LcMeCXdL03LS3vMvApYlJYfbMs8yTjKnPl+4DsTjJ2x80zOzI3tPwQ2ldLPXfNb6eN5od5cI7Pmdn2cn2yuj/1ldn2c+8xZa2Tpn8RdB8xExCsR8RbwNLCu0WYd8ERafga4WVLbj6tmERHHI2JPWv4HcBC4pK88c2wd8LOo7AKWSFred6jkZuDliPh/f0B+zkXEC1RXq6urj9sngC+17PoFYGdE/C0i/g7sBFbPW9CatswRsSMi3kkPd1H93lUxOvp5EpPMM/NiXOY0j30N+HmOLJMYM78VPZ4XMNfIcrg+ngXXxzxcH/PIXSNLP4i7BHi19vgYsyf7/7RJL6BTwIezpDuDdNrKNcDuls03SNorabukT2YN1i2AHZL+JOnbLdsneT76sp7uF3OJff2RiDgO1YseuKilTcn9fTvVu85tzjSOcrs7neKyueMUhlL7+XPAiYg40rG9135uzG9DH89D5RqZj+tjPkOfT1wf51/R9RHy1MjSD+La3i1sXk5zkjbZSfoQ8Cxwb0ScbmzeQ3Vaw6eAHwG/zp2vw2ciYiVwC3CXpBsb20vt63OBj31cFwAAAsdJREFUtcAvWzaX2teTKLW/76P6vaunOpqcaRzl9BPgCuBq4DjV6RdNRfYz8HXGv8vYWz+fYX7r3K1lXQn9PGSukfm4Ppal1P52fcyj2PoI+Wpk6Qdxx4BLa48/Brze1UbSImAxZ/eR8ZyRdA7Vk/dURPyquT0iTkfEP9PyNuAcSRdmjjlLRLye7k8CW6k+Rq+b5Pnowy3Anog40dxQal8DJ0an2qT7ky1tiuvv9EXbLwK3RTqJu2mCcZRNRJyIiHcj4j3g0Y4sJfbzIuArwJauNn31c8f8NsjxvAC4Rmbi+pjVIOcT18c8Sq6PkLdGln4Q90fgSkkfT+8mrQemGm2mgNEVXL4K/L7rxZNDOk/3ceBgRDzU0ebi0XcSJF1H9Tz8NV/K1kwflHTBaJnqS7r7G82mgG+ocj1wavTxcM8635Epsa+T+rjdAPympc3vgFWSlqbTHFaldb2QtBr4LrA2Iv7V0WaScZRN4zspX+7IMsk8k9vngUMRcaxtY1/9PGZ+G9x4XiBcIzNwfcxucPOJ62NWRdbH9O/lrZGR+cot/+uN6opPh6mujnNfWvc9qhcKwPlUpwnMAH8AVvSc97NUH3/uA6bTbQ1wB3BHanM3cIDqKj+7gE8X0M8rUp69Kduor+u5Bfw4PRcvAdcWkPsDVEVncW1dUX1NVUCPA29TvdPyLarvpDwHHEn3y1Lba4HHavvensb2DPDNnjPPUJ2vPRrXoyvefRTYNm4c9Zj5yTRW91FNosubmdPjWfNMX5nT+p+OxnCtbe/9PGZ+K3o8L+Rb29jFNXKuM7s+zl9G18f+Mrs+zn3mrDVSaSczMzMzMzMbgNJPpzQzMzMzM7MaH8SZmZmZmZkNiA/izMzMzMzMBsQHcWZmZmZmZgPigzgzMzMzM7MB8UGcmZmZmZnZgPggzszMzMzMbEB8EGdmZmZmZjYg/wbmRIscYdWrIAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1080x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "val=np.sqrt(lambda_function_np(mB_mass, mKst_mass, q2_range))*np.sqrt(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": 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, qvec_B/qE)\n",
    "        \n",
    "        costheta_l= tf.expand_dims(get_costheta_l(p1q, z), axis=1)\n",
    "        \n",
    "        #pdf=costheta_l*beta_l(q2, ml)\n",
    "        \n",
    "        #pdf = 1./(1.-beta_l(q2, ml)*costheta_l)        \n",
    "        #pdf = (1/(q2))*(1+ztf.square(costheta_l))*(1/(q2+21e6))#*ztf.exp(-ztf.square(q2-9e6))\n",
    "        #pdf = ztf.sqrt(lambda_function(mB, mKst, q2))*beta_l(q2,ml)\n",
    "        #pdf = tf.ones_like(q2) \n",
    "        #pdf = q2/(1.-beta_l(q2, ml)*costheta_l)\n",
    "        pdf = 1/(1.-costheta_l*beta_l(q2,mmu_mass))\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",
    "#lower = ((lower_Kstx, lower_Ksty,  lower_Kstz, lower_KstE, lower_l1x, lower_l1y,  lower_l1z, lower_l1E, lower_l2x, lower_l2y,  lower_l2z, lower_l2E,  ),)\n",
    "#upper = ((upper_Kstx, upper_Ksty,  upper_Kstz, upper_KstE, upper_l1x, upper_l1y,  upper_l1z, upper_l1E, upper_l2x, upper_l2y,  upper_l2z, upper_l2E,  ),)\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))\n",
    "\n"
   ]
  },
  {
   "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=100000)"
   ]
  },
  {
   "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+17YcXAAAgAElEQVR4nO3df7DddX3n8eer/LC2tiUo2hiSDV3TjuhMlb0D7LrTYcXyy05DZ2Qau1Mzlmm6HWh1t90V7B+0Wjpxt9XVrcVNS2roWCNFHTIWSyPKOM6UH0EpguhyC6zckkLcIOo4pQ197x/nc/EkOeeS3Hvuufd8z/Mxc+Z+z/v7+Z77+X4S8uF9Pj++qSokSZIkSd3wfStdAUmSJEnS6JjkSZIkSVKHmORJkiRJUoeY5EmSJElSh5jkSZIkSVKHmORJkiRJUoecuNIVWKyXvOQltXHjxpWuhiRpmd1zzz3fqKrTVroek8L+UZKmx7A+cmKTvI0bN7Jv376VroYkaZkl+b8rXYdJYv8oSdNjWB/pdE1JkiRJ6hCTPEmSJEnqEJM8SZIkSeoQkzxJkiRJ6hCTPEmSJEnqEJM8SZIkSeoQkzxJkiRJ6hCTPEmSJEnqEJM8SZIkSeoQkzxJkiRJ6hCTPEmSJEnqkBNXugKSjrbxqr88rvKPbn/jMtVEkqTVZVAfOawfHNaf2m+q6xzJkyRpxJJ8f5K7kvxtkgeS/E6LfzjJI0nuba/XtHiSfCDJbJL7kpzV91lbkzzUXltX6p4kSZPDkTxJkkbvGeD1VfWdJCcBX0jy6Xbuv1bVTUeUvxjY1F7nANcB5yQ5FbgGmAEKuCfJnqp6aix3IUmaSI7kSZI0YtXznfb2pPaqBS7ZDNzQrrsDOCXJWuBCYG9VHWyJ3V7gouWsuyRp8j1vkpdkZ5Ink9zfFzs1yd42dWRvkjUtftzTTZL8myRfbtd8IElGfZOSJI1bkhOS3As8SS9Ru7Odurb1ke9L8oIWWwc81nf5XIsNix/5u7Yl2Zdk34EDB0Z+L5KkyXIs0zU/DPwhcENf7CrgtqranuSq9v4dLG66yXXANuAO4BZ631B+GmkKHO8GK5ImR1U9C7wmySnAJ5O8Grga+AfgZGAHvb7zXcCgLzhrgfiRv2tH+zxmZmYWGjGUJE2B503yqurzSTYeEd4MnNeOdwG30+uonptuAtyRZH66yXm06SYASfYCFyW5HfjhqvqbFr8BuBSTPOm4HM9OY5LGq6q+2fq7i6rq91v4mSR/Cvxmez8HrO+77HTg8RY/74j47ctZX0nS5FvsmryXVdV+gPbzpS1+vNNN1rXjI+OSJE2sJKe1ETySvBB4A/DV9sUnbWnCpcD8Uog9wFvasodzgadb/3orcEGSNW1pxAUtJknSUKPeXfN4p5sc0zSU5z482UZvaicbNmxYTP0kSRqHtcCuJCfQ+0L1xqr6VJLPJjmNXv93L/CfWvlbgEuAWeC7wFsBqupgkncDd7dy75qfFSN1nUsapMVbbJL3RJK1VbW/fSv5ZIsf73STuXZ8ZPmBXHMgSZoEVXUf8NoB8dcPKV/AFUPO7QR2jrSC0pTzIenqusUmeXuArcD29vPmvviVSXbT23jl6ZYI3gr83vwunPSmm1zdvqH8dpuacifwFuB/LbJOkiRJmkKO+kmHe94kL8lH6Y3CvSTJHL1dMrcDNya5HPg6cFkrvpjpJr9KbwfPF9LbcMVNVyRJkiRpkY5ld803Dzl1/oCyxz3dpKr2Aa9+vnpIk85vGSVJkjQOi91dU5IkSZK0Co16d01Jq4SLyiVJkqaTSZ4kSZJWjMsZpNEzyZMkSZIYnHA6A0aTyDV5kiRJktQhJnmSJEmS1CEmeZIkSZLUIa7Jk0bMBeSSJElaSY7kSZIkSVKHmORJkiRJUoc4XVOaMj4kXZIkqdtM8iRJkrTsXLMujY/TNSVJkiSpQxzJkyRJ0kg5aietLEfyJEmSJKlDTPIkSZIkqUNM8iRJkiSpQ1yTJy2S6w0kSZK0GjmSJ0mSJEkdYpInSZIkSR3idE1JkiRpiGHLMx7d/sYx10Q6do7kSZIkSVKHOJInCfCbSkmSpK5wJE+SJEmSOsQkT5KkEUvy/UnuSvK3SR5I8jstfkaSO5M8lORjSU5u8Re097Pt/Ma+z7q6xb+W5MKVuSNJ0iRxuqYkSaP3DPD6qvpOkpOALyT5NPBfgPdV1e4kHwIuB65rP5+qqlck2QK8B/j5JGcCW4BXAS8HPpPkx6vq2ZW4KUnf4zIHrWaO5EmSNGLV85329qT2KuD1wE0tvgu4tB1vbu9p589PkhbfXVXPVNUjwCxw9hhuQZI0wRzJk57HsG/qJGkhSU4A7gFeAXwQ+Dvgm1V1qBWZA9a143XAYwBVdSjJ08CLW/yOvo/tv0aSpIFM8iQtaFCS61QU6fm1KZWvSXIK8EnglYOKtZ8Zcm5Y/DBJtgHbADZs2LCo+kqSusPpmpIkLaOq+iZwO3AucEqS+S9YTwceb8dzwHqAdv5HgIP98QHX9P+OHVU1U1Uzp5122nLchiRpgpjkSZI0YklOayN4JHkh8AbgQeBzwJtasa3Aze14T3tPO//ZqqoW39J23zwD2ATcNZ67kCRNKqdrSpI0emuBXW1d3vcBN1bVp5J8Bdid5HeBLwHXt/LXA3+WZJbeCN4WgKp6IMmNwFeAQ8AV7qyp1cR169LqZJInSdKIVdV9wGsHxB9mwO6YVfWPwGVDPuta4NpR11GS1F1O15QkSZKkDjHJkyRJkqQOMcmTJEmSpA5xTZ7UuHhckiRJXeBIniRJkiR1iCN5ko7bsFHPR7e/ccw1kSRJ0pEcyZMkSZKkDjHJkyRJkqQOWVKSl+Q/J3kgyf1JPprk+5OckeTOJA8l+ViSk1vZF7T3s+38xr7PubrFv5bkwqXdkiRJkiRNr0UneUnWAb8OzFTVq4ETgC3Ae4D3VdUm4Cng8nbJ5cBTVfUK4H2tHEnObNe9CrgI+KMkJyy2XpIkSZI0zZa68cqJwAuT/DPwA8B+4PXAL7Tzu4DfBq4DNrdjgJuAP0ySFt9dVc8AjySZBc4G/maJdZMkSZLGys3JtBoseiSvqv4e+H3g6/SSu6eBe4BvVtWhVmwOWNeO1wGPtWsPtfIv7o8PuEaSJEmSdByWMl1zDb1RuDOAlwM/CFw8oGjNXzLk3LD4oN+5Lcm+JPsOHDhw/JWWJEmSpI5bynTNNwCPVNUBgCSfAP4dcEqSE9to3enA4638HLAemEtyIvAjwMG++Lz+aw5TVTuAHQAzMzMDE0HpWAybSiFJkiRNuqXsrvl14NwkP9DW1p0PfAX4HPCmVmYrcHM73tPe085/tqqqxbe03TfPADYBdy2hXpIkSZI0tRY9kldVdya5CfgicAj4Er1Rtr8Edif53Ra7vl1yPfBnbWOVg/R21KSqHkhyI70E8RBwRVU9u9h6SZIkabScASNNliXtrllV1wDXHBF+mN7umEeW/UfgsiGfcy1w7VLqImnluaOYJE02kzmpG5b0MHRJkiRJ0uqy1OfkSZIkSXoeg0ZJnemi5eJIniRJkiR1iEmeJEmSJHWISZ4kSZIkdYhJniRJkiR1iBuvqNPcCnp18NEKkiRJ4+NIniRJkiR1iEmeJEmSJHWISZ4kSZIkdYhr8iRJGrEk64EbgB8F/gXYUVXvT/LbwC8DB1rRd1bVLe2aq4HLgWeBX6+qW1v8IuD9wAnAn1TV9nHei6Tl45p1LReTPEmSRu8Q8BtV9cUkPwTck2RvO/e+qvr9/sJJzgS2AK8CXg58JsmPt9MfBH4amAPuTrKnqr4ylruQJE0kkzxJkkasqvYD+9vxt5M8CKxb4JLNwO6qegZ4JMkscHY7N1tVDwMk2d3KmuRJkoZyTZ4kScsoyUbgtcCdLXRlkvuS7EyypsXWAY/1XTbXYsPikiQNZZInSdIySfIi4OPA26vqW8B1wL8GXkNvpO8P5osOuLwWiB/5e7Yl2Zdk34EDBwZcIkmaJk7XVCf40HNJq02Sk+gleB+pqk8AVNUTfef/GPhUezsHrO+7/HTg8XY8LP6cqtoB7ACYmZk5KgmUJE0XkzxJkkYsSYDrgQer6r198bVtvR7AzwH3t+M9wJ8neS+9jVc2AXfRG8nblOQM4O/pbc7yC+O5C3WZX45K3WaSJ0nS6L0O+EXgy0nubbF3Am9O8hp6Uy4fBX4FoKoeSHIjvQ1VDgFXVNWzAEmuBG6l9wiFnVX1wDhvRJI0eUzyJK2YQd8k+2wgdUFVfYHB6+luWeCaa4FrB8RvWeg6SZKO5MYrkiRJktQhJnmSJEmS1CEmeZIkSZLUIa7J00RxNzBJkiRpYY7kSZIkSVKHOJInSZLUYc6CkaaPSZ6kVWXY/4z4aAVJkqRj43RNSZIkSeoQR/IkSZKkVcRZLVoqR/IkSZIkqUNM8iRJkiSpQ0zyJEmSJKlDXJOnVcstnyVJkqTj50ieJEmSJHWISZ4kSZIkdYjTNSVJkqQJ4KMVdKwcyZMkSZKkDjHJkyRJkqQOcbqmpIngFBVJkqRj40ieJEmSJHWISZ4kSZIkdYhJniRJkiR1iEmeJEmSJHXIkpK8JKckuSnJV5M8mOTfJjk1yd4kD7Wfa1rZJPlAktkk9yU5q+9ztrbyDyXZutSbkiRJkqRptdTdNd8P/FVVvSnJycAPAO8Ebquq7UmuAq4C3gFcDGxqr3OA64BzkpwKXAPMAAXck2RPVT21xLpJkiRNjWG7EKv73IFaR1r0SF6SHwZ+CrgeoKr+qaq+CWwGdrViu4BL2/Fm4IbquQM4Jcla4EJgb1UdbIndXuCixdZLkiRJkqbZUkbyfgw4APxpkp8E7gHeBrysqvYDVNX+JC9t5dcBj/VdP9diw+JHSbIN2AawYcOGJVRdq4nfPGop/PZSkiTpcEtJ8k4EzgJ+raruTPJ+elMzh8mAWC0QPzpYtQPYATAzMzOwjCRJUpf55aik57OUjVfmgLmqurO9v4le0vdEm4ZJ+/lkX/n1fdefDjy+QFySpImUZH2Sz7VNyR5I8rYWd3MySdKyW3SSV1X/ADyW5Cda6HzgK8AeYL4T2grc3I73AG9pHdm5wNNtWuetwAVJ1rTO7oIWkyRpUh0CfqOqXgmcC1yR5Ex6M15uq6pNwG18bwZM/+Zk2+htTkbf5mTnAGcD18wnhpIkDbPU3TV/DfhI21nzYeCt9BLHG5NcDnwduKyVvQW4BJgFvtvKUlUHk7wbuLuVe1dVHVxivSRJWjHtS8z59enfTvIgvfXmm4HzWrFdwO30dqB+bnMy4I72iKK1reze+X4xyfzmZB8d281IkibOkpK8qrqX3qMPjnT+gLIFXDHkc3YCO5dSF0mSVqMkG4HXAneyTJuTuTGZJKnfkh6GLkmShkvyIuDjwNur6lsLFR0QO+bNyapqR1XNVNXMaaedtrjKSpI6wyRPkqRlkOQkegneR6rqEy3s5mSSpGW31DV50jFzy2dJ0yJJgOuBB6vqvX2n5jcn287Rm5NdmWQ3vU1Wnm7TOW8Ffq9vs5ULgKvHcQ+SpMllkidJ0ui9DvhF4MtJ7m2xd9JL7tycTJK0rEzyJHXSsJHjR7e/ccw10TSqqi8weD0duDmZJGmZuSZPkiRJkjrEJE+SJEmSOsQkT5IkSZI6xDV5kiRJUgcNWp/u2vTp4EieJEmSJHWISZ4kSZIkdYjTNTVyPvRckqSlsz+VtFiO5EmSJElSh5jkSZIkSVKHOF1TkiRJmhLDpgG762a3mORJmipuJy1JkrrO6ZqSJEmS1CEmeZIkSZLUISZ5kiRJktQhJnmSJEmS1CFuvCJJkrTCfPC5pFFyJE+SJEmSOsQkT5IkSZI6xOmaWhKnl6gLfDCsJEnqEpM8SZIkacr5hWe3OF1TkiRJkjrEkTxJkqQxcZmDpHFwJE+SJEmSOsQkT5IkSZI6xCRPkiRJkjrEJE+SJEmSOsSNV3RMXCguSZIkTQaTPEkawmcGabGS7AR+Bniyql7dYr8N/DJwoBV7Z1Xd0s5dDVwOPAv8elXd2uIXAe8HTgD+pKq2j/M+JEmTyemakiSN3oeBiwbE31dVr2mv+QTvTGAL8Kp2zR8lOSHJCcAHgYuBM4E3t7KSJC3IkTxJkkasqj6fZOMxFt8M7K6qZ4BHkswCZ7dzs1X1MECS3a3sV0ZcXUlSxziSJ0nS+FyZ5L4kO5OsabF1wGN9ZeZabFhckqQFOZInSdJ4XAe8G6j28w+AXwIyoGwx+IvYGvTBSbYB2wA2bNgwirpKEuD69EllkidJ0hhU1RPzx0n+GPhUezsHrO8rejrweDseFj/ys3cAOwBmZmYGJoIaL3ellrSSnK4pSdIYJFnb9/bngPvb8R5gS5IXJDkD2ATcBdwNbEpyRpKT6W3OsmecdZYkTSZH8iRJGrEkHwXOA16SZA64BjgvyWvoTbl8FPgVgKp6IMmN9DZUOQRcUVXPts+5EriV3iMUdlbVA2O+FUnSBDLJ02GcXiJJS1dVbx4Qvn6B8tcC1w6I3wLcMsKqSZKmwJKna7Zn+Xwpyafa+zOS3JnkoSQfa1NMaNNQPpZktp3f2PcZV7f415JcuNQ6SZIkSdK0GsVI3tuAB4Efbu/fQ+9hr7uTfAi4nN6OYpcDT1XVK5JsaeV+/oiHwL4c+EySH5+fqiJJq407jUmSpNVsSSN5SU4H3gj8SXsf4PXATa3ILuDSdry5vaedP7+Vf+4hsFX1CND/EFhJkiRJ0nFY6nTN/wn8N+Bf2vsXA9+sqkPtff+DW597qGs7/3Qr78NeJUmSJGlEFp3kJfkZ4Mmquqc/PKBoPc+5ha458nduS7Ivyb4DBw4cV30lSZIkaRosZSTvdcDPJnkU2E1vmub/BE5JMr/Wr//Brc897LWd/xHgIAs/BPYwVbWjqmaqaua0005bQtUlSZIkqZsWvfFKVV0NXA2Q5DzgN6vqPyb5C+BN9BK/rcDN7ZI97f3ftPOfrapKsgf48yTvpbfxyvxDYCVJklY1Hz0kaTVajufkvQPYneR3gS/xvecCXQ/8WZJZeiN4W2Dhh8BKkiRJko7PSJK8qroduL0dP8yA3TGr6h+By4ZcP/AhsFo+fvMoSZKkxRr0/5I+Smj1WI6RPEmaSnZ4kiRpNVjqIxQkSZIkSauISZ4kSZIkdYhJniRJkiR1iEmeJEmSJHWIG69IkiQdA3emlhY27L8RNyEbP0fyJEmSJKlDTPIkSZIkqUOcrjkFnF4iSZIkTQ9H8iRJkiSpQ0zyJEmSJKlDnK4pScvIncYkSdK4OZInSZIkSR3iSJ4kSVIfNyyTNOkcyZMkSZKkDjHJkyRJkqQOcbqmJEmSpGXjJmTj50ieJEmSJHWISZ4kSSOWZGeSJ5Pc3xc7NcneJA+1n2taPEk+kGQ2yX1Jzuq7Zmsr/1CSrStxL5KkyeN0zQ5xNzBJWjU+DPwhcENf7CrgtqranuSq9v4dwMXApvY6B7gOOCfJqcA1wAxQwD1J9lTVU2O7C0nSRDLJkyRpxKrq80k2HhHeDJzXjncBt9NL8jYDN1RVAXckOSXJ2lZ2b1UdBEiyF7gI+OgyV1+SxsK1esvHJE+SVoAd21R6WVXtB6iq/Ule2uLrgMf6ys212LC4JEkLck2eJEkrKwNitUD86A9ItiXZl2TfgQMHRlo5SdLkcSRPkqTxeCLJ2jaKtxZ4ssXngPV95U4HHm/x846I3z7og6tqB7ADYGZmZmAiqKO5ll1SVzmSJ0nSeOwB5nfI3Arc3Bd/S9tl81zg6Tat81bggiRr2k6cF7SYJEkLciRPkqQRS/JReqNwL0kyR2+XzO3AjUkuB74OXNaK3wJcAswC3wXeClBVB5O8G7i7lXvX/CYskiQtxCRPkqQRq6o3Dzl1/oCyBVwx5HN2AjtHWDVJ0hRwuqYkSZIkdYgjeRPIheJSd/loBUmStFSO5EmSJElSh5jkSZIkSVKHOF1TkiR1msscJE0bkzxJkiRJq8agL2Zcm358nK4pSZIkSR1ikidJkiRJHWKSJ0mSJEkdYpInSZIkSR3ixiurmLuBSZIkSTpeJnmSJEmSVrVhgx/uujmYSZ4kTQC3k5YkScfKJE+SJHWGSx0kaQkbryRZn+RzSR5M8kCSt7X4qUn2Jnmo/VzT4knygSSzSe5LclbfZ21t5R9KsnXptyVJkiRJ02kpu2seAn6jql4JnAtckeRM4CrgtqraBNzW3gNcDGxqr23AddBLCoFrgHOAs4Fr5hNDSZIkSdLxWXSSV1X7q+qL7fjbwIPAOmAzsKsV2wVc2o43AzdUzx3AKUnWAhcCe6vqYFU9BewFLlpsvSRJkiRpmo3kOXlJNgKvBe4EXlZV+6GXCAIvbcXWAY/1XTbXYsPikiRJkqTjtOSNV5K8CPg48Paq+laSoUUHxGqB+KDftY3eVE82bNhw/JVdpVwkLkmSJGlUlpTkJTmJXoL3kar6RAs/kWRtVe1v0zGfbPE5YH3f5acDj7f4eUfEbx/0+6pqB7ADYGZmZmAiKEmSus8vSCVpuEUneekN2V0PPFhV7+07tQfYCmxvP2/ui1+ZZDe9TVaebongrcDv9W22cgFw9WLrJUnTwgfDSpKmnX3hYEsZyXsd8IvAl5Pc22LvpJfc3ZjkcuDrwGXt3C3AJcAs8F3grQBVdTDJu4G7W7l3VdXBJdRLkiRJkqbWopO8qvoCg9fTAZw/oHwBVwz5rJ3AzsXWRZIkSZLUM5LdNSVJkiRJq4NJniRJkiR1iEmeJEmSJHXIkp+TJ0mSJEmrybTvummSN2Y+10eSJEnScnK6piRJY5Tk0SRfTnJvkn0tdmqSvUkeaj/XtHiSfCDJbJL7kpy1srWXJE0CR/IkqWOmfYrKhPgPVfWNvvdXAbdV1fYkV7X37wAuBja11znAde2nJElDOZInSdLK2wzsase7gEv74jdUzx3AKUnWrkQFJUmTw5E8SZLGq4C/TlLA/66qHcDLqmo/QFXtT/LSVnYd8FjftXMttr//A5NsA7YBbNiwYZmrP16uZZek42eSJ0nSeL2uqh5vidzeJF9doGwGxOqoQC9R3AEwMzNz1HlJ0nRxuqYkSWNUVY+3n08CnwTOBp6Yn4bZfj7Zis8B6/suPx14fHy1lSRNIpM8SZLGJMkPJvmh+WPgAuB+YA+wtRXbCtzcjvcAb2m7bJ4LPD0/rVOSpGGcrrlMXEMgSRrgZcAnk0CvD/7zqvqrJHcDNya5HPg6cFkrfwtwCTALfBd46/irPB72m5LGYdC/NV3cfdokT5KkMamqh4GfHBD/f8D5A+IFXDGGqkmSOsQkT5KmhM/PkyRpOrgmT5IkSZI6xCRPkiRJkjrEJE+SJEmSOsQkT5IkSZI6xI1XJEnS2PioBElafo7kSZIkSVKHOJK3RH4jKUmSJE2uLj5iyJE8SZIkSeoQR/IkacoN+gZzkr+9lCRp2jmSJ0mSJEkdYpInSZIkSR1ikidJkiRJHeKaPEmStCzcgVqSVoYjeZIkSZLUIY7kHSO/jZQkSZKmxyQ/P88kT5J0lEnu2CRJmnZO15QkSZKkDjHJkyRJkqQOcbqmJElaEtetS5omk7CkwZE8SZIkSeoQkzxJkiRJ6hCnaw7gtBNJkiRJk8okT5IkSZKWaDWt1TPJkyQds9XUgWn8nOkiSZPBNXmSJEmS1CEmeZIkSZLUIasmyUtyUZKvJZlNctVK10eSpNXA/lGSdLxWxZq8JCcAHwR+GpgD7k6yp6q+spy/17UFkqTVbKX6R0nSZFsVSR5wNjBbVQ8DJNkNbAbsxCRpArghy7JZkf7RL0ElaXQG/Zu63P3japmuuQ54rO/9XItJkjTN7B8lScdttYzkZUCsjiqUbAO2tbffSfK1Za3V+L0E+MZKV2IVsT0OZ3sczTY53Kprj7xnJB/zr0byKZPJ/vF7Vt3f7xVmexzNNjmc7XG0VdMmI+ofYUgfuVqSvDlgfd/704HHjyxUVTuAHeOq1Lgl2VdVMytdj9XC9jic7XE02+Rwtkcn2T82/v0+nO1xNNvkcLbH0aapTVbLdM27gU1JzkhyMrAF2LPCdZIkaaXZP0qSjtuqGMmrqkNJrgRuBU4AdlbVAytcLUmSVpT9oyRpMVZFkgdQVbcAt6x0PVZYp6faLILtcTjb42i2yeFsjw6yf3yOf78PZ3sczTY5nO1xtKlpk1QdtX5bkiRJkjShVsuaPEmSJEnSCJjkjUmS/5Hkq0nuS/LJJKf0nbs6yWySryW5sC9+UYvNJrmqL35GkjuTPJTkY20x/kRJclmSB5L8S5KZI85NXXs8n2H33jVJdiZ5Msn9fbFTk+xtf757k6xp8ST5QGuT+5Kc1XfN1lb+oSRbV+JeRiHJ+iSfS/Jg++/lbS0+tW2ibrKPPJr95LGblj4S7CePZD+5gKryNYYXcAFwYjt+D/Cednwm8LfAC4AzgL+jt7j+hHb8Y8DJrcyZ7ZobgS3t+EPAr670/S2iPV4J/ARwOzDTF5/K9niethp67117AT8FnAXc3xf778BV7fiqvv92LgE+Te85YucCd7b4qcDD7eeadrxmpe9tke2xFjirHf8Q8H/afyNT2ya+uvmyjxzYJvaTx9ZOU9NHtvu1nzy8Pewnh7wcyRuTqvrrqjrU3t5B71lHAJuB3VX1TFU9AswCZ7fXbFU9XFX/BOwGNicJ8Hrgpnb9LuDScd3HqFTVg1U16GG9U9kez2Pgva9wnZZFVX0eOHhEeDO9P1c4/M93M3BD9dwBnJJkLXAhsLeqDlbVU8Be4KLlr/3oVdX+qvpiO/428CCwjiluE3WTfeTR7CeP2dT0kWA/eST7yeFM8lbGL9H7FgF6fxEf6zs312LD4i8GvtnXGc7Hu8L2ONqwe58WL6uq/dD7xxx4aYsf79+ViZZkI/Ba4E5sE3WbfeTCbJPD+e+bfQJgP3mkVfMIhS5I8hngRwec+q2qurmV+S3gEPCR+csGlC8GJ+C1QPlV51jaY9BlA2KdaI8lmIZ7XIxh7dK59kryIuDjwNur6lu9L+YHFx0Q62SbaPLYRx7NfnIkun5/SzE1fYL95NFM8kaoqt6w0Pm2iPNngPOrav4vzhywvq/Y6cDj7XhQ/IuRgx0AAAHHSURBVBv0hpZPbN/K9ZdfVZ6vPYbobHsswUJtMg2eSLK2qva3KRVPtviwdpkDzjsifvsY6rkskpxEr+P6SFV9ooWnuk00mewjj2Y/ORLT3kfClPcJ9pODOV1zTJJcBLwD+Nmq+m7fqT3AliQvSHIGsAm4C7gb2NR2xDoZ2ALsaR3f54A3teu3AsO+7ZtEtsfRBt77CtdpnPbQ+3OFw/989wBvaTtlnQs83aZk3ApckGRN203rghabOG0tzfXAg1X13r5TU9sm6ib7yONimxxu2vtImOI+wX5yAce7U4uvxb3oLYx+DLi3vT7Ud+636O0M9TXg4r74JfR2Cfo7elM35uM/Ru8f9FngL4AXrPT9LaI9fo7etybPAE8At05zexxDew289669gI8C+4F/bn8/Lqe3nuQ24KH289RWNsAHW5t8mcN3n/ul9vdhFnjrSt/XEtrj39ObLnJf378dl0xzm/jq5ss+cmCb2E8ee1tNRR/Z7tV+8vD2sJ8c8kq7KUmSJElSBzhdU5IkSZI6xCRPkiRJkjrEJE+SJEmSOsQkT5IkSZI6xCRPkiRJkjrEJE+SJEmSOsQkT5IkSZI6xCRPkiRJkjrk/wPQXKl/yWK9ugAAAABJRU5ErkJggg==\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",
    "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",
    "#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": "iVBORw0KGgoAAAANSUhEUgAAA3MAAAEvCAYAAADvmpjfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAcIklEQVR4nO3dfZBd510f8O8PCTsFWgcclYJtRaY2zChA02Sx6fDSTN0YmVAErV3LYcBt3REpuFOm5cVpi5u6dGpTwE0bt1SN3TqGYlPTUM0gMBncVya4lkOAyMFFUQXeOJPYkcepCcYo/PrHvWFubnats9ld7R7dz2dGo/PyO7vP6ujce7/7POc81d0BAABgXD5rqxsAAADA2glzAAAAIyTMAQAAjJAwBwAAMELCHAAAwAgJcwAAACO0c6sbMO8Vr3hF79mzZ6ubAQAAsCUee+yxZ7p715nqtl2Y27NnT44ePbrVzQAAANgSVfXbQ+oMswQAABghYQ4AAGCEhDkAAIAREuYAAABGSJgDAAAYIWEOAABghIQ5AACAERLmAAAARkiYAwAAGCFhDgAAYISEOQAAgBHaudUNYHPsueXnBteevP0Nm9gSAABgM+iZAwAAGCE9c1tMDxoAAPCZ0DMHAAAwQnrmRmQtvXgAAMC5Tc8cAADACOmZw317AAAwQnrmAAAARkiYAwAAGCFhDgAAYISEOQAAgBEa9ACUqtqX5K1JdiR5e3ffPrf//CTvSPLaJB9Ncn13n6yqz07y9iSvmX6vd3T3P9vA9rNNeagKAABsrjP2zFXVjiR3Jbkmyd4kN1TV3rmym5I8292XJbkzyR3T7dclOb+7vyKToPedVbVnY5oOAACwuIb0zF2R5Hh3n0iSqro/yf4kj8/U7E/ylunyg0neVlWVpJN8blXtTPLHkryY5GMb03QWkR4/AACYGHLP3EVJnpxZX55uW7Gmu08neS7JhZkEu99N8qEkv5PkR7r71DrbDAAAsPCGhLlaYVsPrLkiySeSfHGSS5P8var6kk/7BlUHq+poVR19+umnBzQJAABgsQ0ZZrmc5JKZ9YuTPLVKzfJ0SOUFSU4leWOSX+juP0jykar65SRLSU7MHtzdh5IcSpKlpaX5oMg2spZhjlvNkEwAAM5lQ8Lco0kur6pLk3wwyYFMQtqsw0luTPLuJNcmebi7u6p+J8lfqKqfSPI5Sb46yb/YqMZzbhhTQAQAgO3ijMMsp/fA3ZzkoSTvT/LT3X2sqm6rqm+elt2d5MKqOp7k7ya5Zbr9riSfl+R9mYTCf9/dv77BPwMAAMDCGTTPXHcfSXJkbtutM8svZDINwfxxz6+0HQAAgPUZ8gAUAAAAthlhDgAAYIQGDbNkbTzQAwAA2Gx65gAAAEZImAMAABghwywhJhgHAGB8hDlYI8EPAIDtQJiDkREmAQBI3DMHAAAwSsIcAADACAlzAAAAIyTMAQAAjJAHoMAm8rASAAA2izAH28Ragh8AABhmCQAAMELCHAAAwAgJcwAAACMkzAEAAIyQMAcAADBCwhwAAMAICXMAAAAjNGieuaral+StSXYkeXt33z63//wk70jy2iQfTXJ9d5+sqm9L8n0zpV+Z5DXd/d6NaDywcUxwDgAwLmfsmauqHUnuSnJNkr1JbqiqvXNlNyV5trsvS3JnkjuSpLt/srtf3d2vTvLtSU4KcgAAAOs3ZJjlFUmOd/eJ7n4xyf1J9s/V7E9y73T5wSRXVVXN1dyQ5KfW01gAAAAmhgyzvCjJkzPry0muXK2mu09X1XNJLkzyzEzN9fn0EAiMkCGZAABbb0iYm+9hS5JeS01VXZnk4939vhW/QdXBJAeTZPfu3QOaBAyxltC1Hdog+AEADDdkmOVykktm1i9O8tRqNVW1M8kFSU7N7D+Qlxhi2d2Hunupu5d27do1pN0AAAALbUiYezTJ5VV1aVWdl0kwOzxXczjJjdPla5M83N2dJFX1WUmuy+ReOwAAADbAGYdZTu+BuznJQ5lMTXBPdx+rqtuSHO3uw0nuTnJfVR3PpEfuwMyX+Poky919YuObDwAAsJgGzTPX3UeSHJnbduvM8guZ9L6tdOx/S/LVn3kTAQAAmDdkmCUAAADbjDAHAAAwQoOGWQLwqUy5AABsNWEO2DYEJACA4YQ5YJQEPwBg0blnDgAAYIT0zA20ll4AAACAzaZnDgAAYISEOQAAgBES5gAAAEbIPXPAOW/oPa+eegkAjImeOQAAgBHSMwcwtR2eWmv+PABgKD1zAAAAIyTMAQAAjJAwBwAAMELumQPYZNvhXjwA4NyjZw4AAGCEhDkAAIAREuYAAABGSJgDAAAYIQ9AARgpE4wDwGIb1DNXVfuq6omqOl5Vt6yw//yqemC6/5Gq2jOz7yur6t1VdayqfqOqXrZxzQcAAFhMZwxzVbUjyV1JrkmyN8kNVbV3ruymJM9292VJ7kxyx/TYnUl+IsmbuvtVSV6X5A82rPUAAAALasgwyyuSHO/uE0lSVfcn2Z/k8Zma/UneMl1+MMnbqqqSXJ3k17v715Kkuz+6Qe0GYItt1vx5hoQCwDBDhllelOTJmfXl6bYVa7r7dJLnklyY5EuTdFU9VFXvqarvX3+TAQAAGNIzVyts64E1O5N8bZKvSvLxJL9UVY919y99ysFVB5McTJLdu3cPaBIAAMBiGxLmlpNcMrN+cZKnVqlZnt4nd0GSU9Pt/727n0mSqjqS5DVJPiXMdfehJIeSZGlpaT4oArBOnnwJAOeeIWHu0SSXV9WlST6Y5ECSN87VHE5yY5J3J7k2ycPd3VX1UJLvr6rPSfJikj+fyQNSANimNuteOABgY50xzHX36aq6OclDSXYkuae7j1XVbUmOdvfhJHcnua+qjmfSI3dgeuyzVfVjmQTCTnKku31KAAAAWKdBk4Z395EkR+a23Tqz/EKS61Y59icymZ4AAACADTJo0nAAAAC2F2EOAABghIQ5AACAERp0zxwAbEemXABgkemZAwAAGCFhDgAAYISEOQAAgBES5gAAAEbIA1AA2FbW8lATAFhkeuYAAABGSJgDAAAYIWEOAABghIQ5AACAERLmAAAARkiYAwAAGCFhDgAAYISEOQAAgBEyaTgAC2Etk5GfvP0Nm9gSANgYwhwArMNaQuJQwiQAQxhmCQAAMEJ65gBgmzEkFIAhBvXMVdW+qnqiqo5X1S0r7D+/qh6Y7n+kqvZMt++pqt+rqvdO//z4xjYfAABgMZ2xZ66qdiS5K8nrkywnebSqDnf34zNlNyV5trsvq6oDSe5Icv103we6+9Ub3G4AAICFNqRn7ookx7v7RHe/mOT+JPvnavYnuXe6/GCSq6qqNq6ZAAAAzBpyz9xFSZ6cWV9OcuVqNd19uqqeS3LhdN+lVfWrST6W5B929/9cX5MBYHNtxhMqAWCjDQlzK/Ww9cCaDyXZ3d0frarXJvnZqnpVd3/sUw6uOpjkYJLs3r17QJMAgLXyYBWAc8uQMLec5JKZ9YuTPLVKzXJV7UxyQZJT3d1Jfj9JuvuxqvpAki9NcnT24O4+lORQkiwtLc0HRQBgFXoRARbXkHvmHk1yeVVdWlXnJTmQ5PBczeEkN06Xr03ycHd3Ve2aPkAlVfUlSS5PcmJjmg4AALC4ztgzN70H7uYkDyXZkeSe7j5WVbclOdrdh5PcneS+qjqe5FQmgS9Jvj7JbVV1Osknkrypu09txg8CAACwSAZNGt7dR5Icmdt268zyC0muW+G4n0nyM+tsIwAAAHMGTRoOAADA9iLMAQAAjJAwBwAAMELCHAAAwAgJcwAAACM06GmWAMBiWctk5Cdvf8OWf12ARaRnDgAAYISEOQAAgBEyzBIAWJe1DJ0EYOPomQMAABghYQ4AAGCEhDkAAIAREuYAAABGSJgDAAAYIWEOAABghExNAAAsjLVMo3Dy9jdsYksA1k/PHAAAwAjpmQMAtiWTkQO8ND1zAAAAIyTMAQAAjJAwBwAAMELCHAAAwAgNCnNVta+qnqiq41V1ywr7z6+qB6b7H6mqPXP7d1fV81X1vRvTbAAAgMV2xjBXVTuS3JXkmiR7k9xQVXvnym5K8mx3X5bkziR3zO2/M8nPr7+5AAAAJMN65q5Icry7T3T3i0nuT7J/rmZ/knunyw8muaqqKkmq6luSnEhybGOaDAAAwJB55i5K8uTM+nKSK1er6e7TVfVckgur6veS/ECS1ycxxBIAOCetZU68k7e/YRNbAiySIT1ztcK2Hljzj5Pc2d3Pv+Q3qDpYVUer6ujTTz89oEkAAACLbUjP3HKSS2bWL07y1Co1y1W1M8kFSU5l0oN3bVX9cJKXJ/nDqnqhu982e3B3H0pyKEmWlpbmgyIAAABzhoS5R5NcXlWXJvlgkgNJ3jhXczjJjUneneTaJA93dyf5uk8WVNVbkjw/H+QAALajtQydBNgKZwxz03vgbk7yUJIdSe7p7mNVdVuSo919OMndSe6rquOZ9Mgd2MxGAwAALLqadKBtH0tLS3306NGtbsan8ds5AOBs87AUWExV9Vh3L52pbsgwSwAAtoCnZAIvZcjTLAEAANhmhDkAAIAREuYAAABGyD1zAADnAPfXweLRMwcAADBCwhwAAMAIGWYJALBghg7JNBwTtjc9cwAAACMkzAEAAIyQYZYAAKzIEzJhe9MzBwAAMELCHAAAwAgJcwAAACMkzAEAAIyQMAcAADBCwhwAAMAICXMAAAAjJMwBAACMkEnDAQA4q0xGDhtDzxwAAMAIDQpzVbWvqp6oquNVdcsK+8+vqgem+x+pqj3T7VdU1Xunf36tqr51Y5sPAACwmM4Y5qpqR5K7klyTZG+SG6pq71zZTUme7e7LktyZ5I7p9vclWeruVyfZl+TfVpWhnQAAAOs0JFhdkeR4d59Ikqq6P8n+JI/P1OxP8pbp8oNJ3lZV1d0fn6l5WZJed4sBANh21nIfHLAxhoS5i5I8ObO+nOTK1Wq6+3RVPZfkwiTPVNWVSe5J8sok397dp9fdagAAFsJmhUQPVuFcMOSeuVph23wP26o13f1Id78qyVcleXNVvezTvkHVwao6WlVHn3766QFNAgAAWGxDwtxykktm1i9O8tRqNdN74i5Icmq2oLvfn+R3k3z5/Dfo7kPdvdTdS7t27RreegAAgAU1JMw9muTyqrq0qs5LciDJ4bmaw0lunC5fm+Th7u7pMTuTpKpemeTLkpzckJYDAAAssDPeMze9B+7mJA8l2ZHknu4+VlW3JTna3YeT3J3kvqo6nkmP3IHp4V+b5Jaq+oMkf5jku7r7mc34QQAAABbJoGkCuvtIkiNz226dWX4hyXUrHHdfkvvW2UYAAADmDJo0HAAAgO1FmAMAABghYQ4AAGCEBt0zBwAA55K1TEZugnG2Kz1zAAAAIyTMAQAAjJAwBwAAMELCHAAAwAgJcwAAACMkzAEAAIyQqQkAAGCDDJ3ywHQHbAQ9cwAAACOkZw4AAF7CWiYYh7NJzxwAAMAICXMAAAAjJMwBAACMkDAHAAAwQsIcAADACHmaJQAAnGVreUKmOelYjZ45AACAERLmAAAARkiYAwAAGKFB98xV1b4kb02yI8nbu/v2uf3nJ3lHktcm+WiS67v7ZFW9PsntSc5L8mKS7+vuhzew/QAAcE5zfx2rOWPPXFXtSHJXkmuS7E1yQ1XtnSu7Kcmz3X1ZkjuT3DHd/kySv9TdX5HkxiT3bVTDAQAAFtmQYZZXJDne3Se6+8Uk9yfZP1ezP8m90+UHk1xVVdXdv9rdT023H0vysmkvHgAAAOswJMxdlOTJmfXl6bYVa7r7dJLnklw4V/NXkvxqd//+Z9ZUAAAAPmnIPXO1wrZeS01VvSqToZdXr/gNqg4mOZgku3fvHtAkAACAxTakZ245ySUz6xcneWq1mqrameSCJKem6xcneWeS7+juD6z0Dbr7UHcvdffSrl271vYTAAAALKAhPXOPJrm8qi5N8sEkB5K8ca7mcCYPOHl3kmuTPNzdXVUvT/JzSd7c3b+8cc0GAADmefLlYjljmOvu01V1c5KHMpma4J7uPlZVtyU52t2Hk9yd5L6qOp5Jj9yB6eE3J7ksyQ9W1Q9Ot13d3R/Z6B8EAAAYTvAbv0HzzHX3kSRH5rbdOrP8QpLrVjjuh5L80DrbCAAAwJwh98wBAACwzQhzAAAAIyTMAQAAjJAwBwAAMELCHAAAwAgNepolAACwuExjsD3pmQMAABghYQ4AAGCEhDkAAIAREuYAAABGSJgDAAAYIWEOAABghExNAAAAbBjTGJw9euYAAABGSJgDAAAYIWEOAABghIQ5AACAERLmAAAARkiYAwAAGCFhDgAAYISEOQAAgBEaNGl4Ve1L8tYkO5K8vbtvn9t/fpJ3JHltko8mub67T1bVhUkeTPJVSf5Dd9+8kY0HAADGywTj63PGnrmq2pHkriTXJNmb5Iaq2jtXdlOSZ7v7siR3Jrljuv2FJD+Y5Hs3rMUAAAAMGmZ5RZLj3X2iu19Mcn+S/XM1+5PcO11+MMlVVVXd/bvd/b8yCXUAAABskCFh7qIkT86sL0+3rVjT3aeTPJfkwo1oIAAAAJ9uSJirFbb1Z1Cz+jeoOlhVR6vq6NNPPz30MAAAgIU1JMwtJ7lkZv3iJE+tVlNVO5NckOTU0EZ096HuXurupV27dg09DAAAYGENCXOPJrm8qi6tqvOSHEhyeK7mcJIbp8vXJnm4uwf3zAEAALA2Z5yaoLtPV9XNSR7KZGqCe7r7WFXdluRodx9OcneS+6rqeCY9cgc+eXxVnUzyJ5KcV1XfkuTq7n58438UAACAxTFonrnuPpLkyNy2W2eWX0hy3SrH7llH+wAAAFjBoDAHAACwlYZOML5Ik4sPuWcOAACAbUaYAwAAGCFhDgAAYISEOQAAgBES5gAAAEZImAMAABghYQ4AAGCEhDkAAIARMmk4AABwzhg6uXgy/gnG9cwBAACMkDAHAAAwQsIcAADACAlzAAAAIyTMAQAAjJAwBwAAMELCHAAAwAgJcwAAACMkzAEAAIyQMAcAADBCwhwAAMAICXMAAAAjNCjMVdW+qnqiqo5X1S0r7D+/qh6Y7n+kqvbM7HvzdPsTVfUNG9d0AACAxXXGMFdVO5LcleSaJHuT3FBVe+fKbkrybHdfluTOJHdMj92b5ECSVyXZl+RfT78eAAAA6zCkZ+6KJMe7+0R3v5jk/iT752r2J7l3uvxgkquqqqbb7+/u3+/u/5vk+PTrAQAAsA5DwtxFSZ6cWV+ebluxprtPJ3kuyYUDjwUAAGCNdg6oqRW29cCaIcemqg4mOThdfb6qnhjQrrPpFUme2epGkMS52E6ci+3Dudg+nIvtw7nYPpyL7cO5mFN3bNm3PtO5eOWQLzIkzC0nuWRm/eIkT61Ss1xVO5NckOTUwGPT3YeSHBrS4K1QVUe7e2mr24FzsZ04F9uHc7F9OBfbh3OxfTgX24dzsX1s1LkYMszy0SSXV9WlVXVeJg80OTxXczjJjdPla5M83N093X5g+rTLS5NcnuR/r7fRAAAAi+6MPXPdfbqqbk7yUJIdSe7p7mNVdVuSo919OMndSe6rquOZ9MgdmB57rKp+OsnjSU4n+e7u/sQm/SwAAAALY8gwy3T3kSRH5rbdOrP8QpLrVjn2nyb5p+to43awbYeALiDnYvtwLrYP52L7cC62D+di+3Autg/nYvvYkHNRk9GQAAAAjMmQe+YAAADYZoS5qaraV1VPVNXxqrplhf3nV9UD0/2PVNWes9/Kc19VXVJV/7Wq3l9Vx6rq76xQ87qqeq6q3jv9c+tKX4uNUVUnq+o3pv/WR1fYX1X1L6fXxq9X1Wu2op3nuqr6spn/8++tqo9V1ffM1bg2NklV3VNVH6mq981s+4KqeldV/db0789f5dgbpzW/VVU3rlTDcKuci39eVb85fQ16Z1W9fJVjX/L1jLVZ5Vy8pao+OPM69I2rHPuSn7tYm1XOxQMz5+FkVb13lWNdFxtotc+ym/WeYZhlkqrakeT/JHl9JtMpPJrkhu5+fKbmu5J8ZXe/qaoOJPnW7r5+Sxp8DquqL0ryRd39nqr640keS/Itc+fidUm+t7u/aYuauVCq6mSSpe5ecS6U6Rv1307yjUmuTPLW7r7y7LVw8Uxfsz6Y5Mru/u2Z7a+La2NTVNXXJ3k+yTu6+8un2344yanuvn36YfTzu/sH5o77giRHkyxlMs/qY0le293PntUf4Byyyrm4OpMnaZ+umswaNX8upnUn8xKvZ6zNKufiLUme7+4feYnjzvi5i7VZ6VzM7f/RJM91920r7DsZ18WGWe2zbJK/lk14z9AzN3FFkuPdfaK7X0xyf5L9czX7k9w7XX4wyVVVtdKk6KxDd3+ou98zXf5/Sd6f5KKtbRVnsD+TN4/u7l9J8vLpCxmb56okH5gNcmyu7v4fmTytedbs+8K9mbxZz/uGJO/q7lPTN+N3Jdm3aQ1dACudi+7+xe4+PV39lUzmtWWTrXJdDDHkcxdr8FLnYvp59a8m+amz2qgF9RKfZTflPUOYm7goyZMz68v59ADxRzXTN4znklx4Vlq3oGoylPXPJnlkhd1/rqp+rap+vqpedVYbtng6yS9W1WNVdXCF/UOuHzbWgaz+puzaOHu+sLs/lEzevJP8yRVqXB9n399I8vOr7DvT6xkb4+bpkNd7VhlK5ro4u74uyYe7+7dW2e+62CRzn2U35T1DmJtYqYdtfvzpkBo2SFV9XpKfSfI93f2xud3vSfLK7v4zSf5Vkp892+1bMF/T3a9Jck2S754O5Zjl2jiLquq8JN+c5D+tsNu1sf24Ps6iqvoHmcxr+5OrlJzp9Yz1+zdJ/nSSVyf5UJIfXaHGdXF23ZCX7pVzXWyCM3yWXfWwFba95LUhzE0sJ7lkZv3iJE+tVlNVO5NckM9saAFnUFWfncl//p/s7v88v7+7P9bdz0+XjyT57Kp6xVlu5sLo7qemf38kyTszGR4za8j1w8a5Jsl7uvvD8ztcG2fdhz85pHj690dWqHF9nCXTBwV8U5Jv61UeCDDg9Yx16u4Pd/cnuvsPk/y7rPxv7Lo4S6afWf9ykgdWq3FdbLxVPstuynuGMDfxaJLLq+rS6W+9DyQ5PFdzOMknnyhzbSY3Wvst0gabjuu+O8n7u/vHVqn5U5+8X7Gqrsjk//FHz14rF0dVfe705t1U1ecmuTrJ++bKDif5jpr46kxusP7QWW7qIln1N6yujbNu9n3hxiT/ZYWah5JcXVWfPx1udvV0GxuoqvYl+YEk39zdH1+lZsjrGes0d8/0t2blf+Mhn7vYGH8xyW929/JKO10XG+8lPstuynvGzvU3efymT7+6OZN/rB1J7unuY1V1W5Kj3X04k5NyX1Udz6RH7sDWtfic9jVJvj3Jb8w8QvfvJ9mdJN3945mE6b9VVaeT/F6SA4L1pvnCJO+c5oOdSf5jd/9CVb0p+aPzcSSTJ1keT/LxJH99i9p6zquqz8nk6W/fObNt9ly4NjZJVf1UktcleUVVLSf5R0luT/LTVXVTkt9Jct20dinJm7r7b3b3qar6J5l8eE2S27rbqI51WOVcvDnJ+UneNX29+pXp06e/OMnbu/sbs8rr2Rb8COeMVc7F66rq1ZkMDTuZ6evV7LlY7XPXFvwI54yVzkV3350V7rF2XWy61T7Lbsp7hqkJAAAARsgwSwAAgBES5gAAAEZImAMAABghYQ4AAGCEhDkAAIAREuYAAABGSJgDAAAYIWEOAABghP4/AbRPxX7C3y8AAAAASUVORK5CYII=\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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}