{
"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"
]
}
],
"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"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def scalar_product_4(v1, v2):\n",
" \"\"\"Calculate mass scalar for Lorentz 4-momentum.\n",
" Arguments:\n",
" vector: Input Lorentz momentum vector.\n",
" \"\"\"\n",
" v0 = v1[:,3]*v2[:,3]\n",
" \n",
" vs = scalar_product_3(v1, v2)\n",
" output = (v0-vs)\n",
" \n",
" return output\n",
"\n",
"def scalar_product_3(v1, v2):\n",
" \"\"\"Calculate scalar product of two 3-vectors.\n",
" Arguments:\n",
" vec1: First vector.\n",
" vec2: Second vector.\n",
" \"\"\"\n",
" output = (v1[:,0]*v2[:,0] + v1[:,1]*v2[:,1] + v1[:,2]*v2[:,2])\n",
" return output\n",
"\n",
"def get_costheta_l(p1, p2): \n",
" \n",
" num = scalar_product_3(p1, p2)\n",
"\n",
" den1= np.sqrt(scalar_product_3(p1, p1))\n",
" den2= np.sqrt(scalar_product_3(p2, p2))\n",
" \n",
" costheta_l = num/(den1*den2)\n",
" \n",
" return costheta_l\n",
"\n",
"def beta_l(q2, m):\n",
" \n",
" return 1-4*(np.square(m)/q2)\n",
"\n",
"\n",
"def lambda_function(m1, m2, q2):\n",
" \n",
" return np.square(np.square(m1)-np.square(m2)) -2*q2*(np.square(m1)+np.square(m2)) + np.square(q2)\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"WARNING:tensorflow:From /Users/davide/miniconda3/envs/zfit_env2/lib/python3.7/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": [
"n_events=1000000\n",
"\n",
"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",
"z=np.array([[0.,0.,1.,0.] for i in range(n_events)])\n",
"weights_np, particles_np = B.generate(n_events=n_events)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(1000000, 4)"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"z.shape"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"105.0000404384508 5279.999786972526\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA2wAAAE6CAYAAABqJNGSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd3jV53n/8c9zlvY+Otp7s8TeGw/sGhMnsWMnduMMO27iJmmSJk3SjGa0PzepM5ukrmMnthM7XvHADhgwe5lhtkBsJECbLbTO+f7+OELIqgABko7G+3VdukDSF7iNAemj+3nu21iWJQAAAABA32MLdAEAAAAAgM4R2AAAAACgjyKwAQAAAEAfRWADAAAAgD6KwAYAAAAAfRSBDQAAAAD6qIAGNmPMU8aYKmPMzi48+zNjzNbWl1JjzKneqBEAAAAAAsUEcg+bMWa6pHOSnrEsa9g1/Lh/lDTKsqxP91hxAAAAABBgAe2wWZa1UlJd+7cZY3KMMQuNMZuNMauMMYWd/ND7JD3fK0UCAAAAQIA4Al1AJ56Q9IhlWfuMMRMk/UbS7IvvNMZkSMqS9G6A6gMAAACAXtGnApsxJlzSZEkvGWMuvjmow2P3SnrZsixvb9YGAAAAAL2tTwU2+Y9onrIsa+QVnrlX0hd6qR4AAAAACJg+Ndbfsqwzkg4ZY+6WJONXfPH9xpgCSTGS1gWoRAAAAADoNYEe6/+8/OGrwBhTboz5jKRPSPqMMWabpF2S5rf7IfdJesEK5GhLAAAAAOglAR3rDwAAAAC4vD51JBIAAAAAcAmBDQAAAAD6qIBNiXS73VZmZmagfnkAAAAACKjNmzfXWJYVf6VnAhbYMjMztWnTpkD98gAAAAAQUMaYI1d7hiORAAAAANBHEdgAAAAAoI8isAEAAABAH0VgAwAAAIA+isAGAAAAAH3UVQObMSbNGLPMGFNijNlljPlSJ8/MNMacNsZsbX35bs+UCwAAAACDR1fG+rdI+qplWVuMMRGSNhtjFluWtbvDc6ssy7qj+0sEAAAAgMHpqh02y7JOWJa1pfX7ZyWVSErp6cIAAAAAYLC7pjtsxphMSaMkbejk3ZOMMduMMX8zxgzthtoAAAAAYFDrypFISZIxJlzSK5K+bFnWmQ7v3iIpw7Ksc8aY2yW9Jimvk5/jYUkPS1J6evp1Fz3QNbZ4VXe+STVnm1R7vlENzV41tvjU1OKT12fJYbfJ5bDJZTdy2m0KD3IoOtSlqBCnokOdCnbaA/2fAAAAAKAbGMuyrv6QMU5JCyQtsizr8S48f1jSWMuyai73zNixY61NmzZdQ6kDj2VZOnbqgjYertO2stPaX3VOB6rP6cTphhv6eYMcNrnDg5QYFazEqGAlRQa3fT81JlSZcaGKDnV1038FAAAAgOthjNlsWdbYKz1z1Q6bMcZI+r2kksuFNWNMoqRKy7IsY8x4+Y9a1l5HzQOeZVnaWnZKb+84oYW7KlRWd0GSFOayK9cTrknZccp0hyk+IkhxYS7FhbsU4nTI5bApyGGT3Wbk9VlqbPGp2evvup1rbNGp+madvtCsUxeadKq+WdVnG1VxukG7j5/R0pJKNTT7PlBHVIhTmXGhyogLa/s2Kz5MeZ5wRQQ7A/FbAwAAAKCDrhyJnCLpAUk7jDFbW9/2LUnpkmRZ1u8kfVTSPxhjWiRdkHSv1ZXW3SDS1OLTa1uP6anVh7Sn4qycdqOpuW59dmq2xmXGqiAxQnab6ZFf27IsnbnQohNnLqis7oKO1J7Xkdp6Ha49r61lp7Rg+3H52v3fSo4KVl5ChPITwlu/jVCeJ1xhQV0+QQsAAACgG3TpSGRPGCxHIi3L0sKdFfqPv+3R0bp6FSZG6MHJmbpteJKiQvpGJ6upxafyk/U6UH1epZVnta/yrEorz2l/9Tk1tVzqzKXGhKgoKVJDkyM1NDlKQ5MjlRQVLH8TFgAAAMC16JYjkbh+teca9Y1XdmhJSaXyE8L19IPjNLMgvs8FHJfDpuz4cGXHh+vmIQltb/f6LB2tq28LcXsqzmr3iTNaUlKpizk/JtSpIe0C3JCkSGXHh/dYtxAAAAAYTAhsPWRr2Sk99Mwmna5v1rdvL9KnpmTKYb+mLQoBZ7cZZbnDlOUO061DE9vefr6xxR/ejp/WruNntPvEGf1h7eG2blyw06ahyVEqTo1WcZr/24y40D4XVAEAAIC+jiORPWDJ7ko9+vwWxUcE6YkHxqooKTLQJfW4Zq9PB6rPaffxM9p57Ix2HDulHcdOtw07iQpxakRqlEamRWtEarSKU6PkiQwOcNUAAABA4HTlSCSBrZutLK3WZ/64UUVJkXrqwXFyhwcFuqSAafH6VFp5TtvLT2lb+SltLTut0sqz8rZOOEmKClZxarRGpPmDXHFqNINNAAAAMGgQ2HrZzmOndffv1inTHaYXHpqoqNC+MVSkL7nQ5NWu46e1rfy0tpWd0vbyUzpcWy9JshmpKClSo9NjNCYjRqPTY5QWG8JRSgAAAAxIBLZedPpCs+b9arWavT698ehUxUcM3s7atTp5vklby0/p/SMntfnoSW09ekrnm7ySJHd4kMZkRLeFuGEpUQp22gNcMQAAAHDjmBLZi/71tZ06fuqC/vK5SYS1axQT5tKsAo9mFXgk+adT7q04qy1HT2pLa4hbtKtSkuS0Gw1NjrrUhcuIVlJUSCDLBwAAAHoMga0bLNtTpTe3HddXbs7XmIyYQJfT79ltRkOSIzUkOVL3T8yQJNWca9SWIye15egpbTlyUn/acERPrTkkyb/oe2xmrMZlxWpcZozyPRGysVYAAAAAAwBHIm9QQ7NXc/5rhUJcdr39xWlyOfrX6P7+qqnFp5ITZ7Tl6EltOnJSmw7XqfJMoyT/RMqxGTEamxmr8VkxGp4Szf8XAAAA9DkciewFz793VMdOXdBzn5lAKOhFLodNxWnRKk6L1qemZMmyLJXVXdDGw3XaeLhO7x2u09I9VZKkoNZnx7d24UanRysimIEwAAAA6PvosN2AC01eTfvPZcrzhOv5hycGuhx0UHuuURsPn2wLcbuOn5HXZ7VNoxyXGet/yYqRJ4KdcAAAAOhddNh62Euby1RzrlG/vX90oEtBJ+LCgzR3WKLmDkuUJJ1vbNH7R0/pvcN12nS4Tn/ZWKY/rD0sScqMC20LcBOyY5UeG8o6AQAAAAQcge06WZalZ9YdUXFqlMZlxga6HHRBWJBDU/PcmprnliQ1e33aeey0Nh0+qfcO12lJSaVe2lwuSUqMDNaE7FhNyIrTxOxYZbnDCHAAAADodQS267TuQK32V53TT+8uDnQpuE5Ou02j0mM0Kj1GD03Pls9n6UD1Oa0/VKcNB2u19kCtXt96XJIUHxGkCVmxmpAdp4lZscr1hBPgAAAA0OMIbNfpxU1ligpx6o4RSYEuBd3EZjPKS4hQXkKEHpiYIcuydKjmvNYfrNOGQ7XacLBOC7afkCS5w10an+XvwE3IjmWVAAAAAHoEge06NDR7tXh3peYVJyvYaQ90Oeghxhhlx4crOz5cH5+QLsuydLSuXusP+sPb+oO1entHhSQpJtT5gQBXlBhJgAMAAMANI7Bdh+V7q3S+yas7RiQHuhT0ImOMMuLClBEXpo+N8we48pMX/AHukD/ALdpVKUmKDHZofOv9twlZcRqSHCk7AQ4AAADXiMB2HRZsP6G4MJcmZjNsZDAzxigtNlRpsaG6e2yaJOnYqQvacLEDd6hWS0r8AS4iyKFxWbGakBWrSTlxGpocRYADAADAVRHYrpHXZ2llabVuHZooh51F2figlOgQfXh0qj48OlWSVHG6QRsO1frvwR2s1buty7wjgh2akBWnSTlxmpQdp8JE7sABAADg/yKwXaPt5ad0pqFF0/PjA10K+oHEqGDNH5mi+SNTJElVZxq07mCt1h2o1bqDlzpwMaFOTcy+FOCYQgkAAACJwHbNVpbWyBhpSq470KWgH/JEfjDAHTt1wR/eDtRq/cFa/W2nf4iJOzyoLbxNzolTRhyLvAEAAAYjAts1Wr2/WsNTohQb5gp0KRgAUqJD9NExqfromNS2KZQXu2/rDtTqzW3+PXBJUcGadLEDlxOn1JjQAFcOAACA3kBguwaNLV5tKz+tT07KCHQpGIDaT6G8d7x/CuWB6vNad7BW6w/UanlptV59/5gkKS02pLX75taknDglRAYHuHoAAAD0BALbNdh9/IyaWnwanR4T6FIwCBhjlOsJV64nXA9MzJDPZ6m06qzWHajV2gO1WrizQi9uKpckZbvD2rpvE7Pj5A4PCnD1AAAA6A4Etmvw/tFTkqRRBDYEgM1mVJgYqcLESH1qSpa8PkslJ860Brgavb71uP604agkqSAhoi28TcyOVXQoR3gBAAD6IwLbNXi/7JSSo4KVGMXxMwSe3WY0LCVKw1Ki9ND0bLV4fdpx7HTb/bcXNh7VH9YeljHSkKRITcl1a3JOnMZnxSrUxV99AACA/oDP2q7B1rKTGpkeHegygE457DaNSo/RqPQYfX5mrppafNpWfkpr9/s7cE+vOaQnVh6U0240Ki1Gk3PjNCXXrZFp0XKyUxAAAKBPIrB10dmGZpXVXdC949IDXQrQJS6HTeMyYzUuM1ZfuilPF5q82ni4Tmv212jNgRr9Yuk+/XzJPoW67BqfFaspOW5Nzo1TUWIkS7wBAAD6CAJbF5VWnpXkvxsE9EchLrum58e3LX0/Vd+kdQdqteZAjdbur9XyvSWSpNgwl38CZW6cpuS42QEHAAAQQAS2LtpbcU6SVJBIYMPAEB3q0m3Dk3Tb8CRJ0onTF7Rmf63Wtnbg3tpxQpJ/V9zkHP/xycm5cfJEcIcTAACgtxDYumhvxRmFuexKiQ4JdClAj0iK+uAS7wPV57X2QI3W7K/Rol0Vemmzf4VAfkK4Jue4NSXXrQnZsYoMdga4cgAAgIGLwNZFeyrOKj8xgrs9GBTa74D7+0mZ8vos7Tp+2t+BO1DTNoHSZqQRqdGa0np8cnRGjIKd9kCXDwAAMGAQ2Lpof9U53VSUEOgygICw24xGpEZrRGq0/mFmjhpbvNpy5FRbB+53Kw7qv5cdUJDDprGZMW0duOEpUbLzRQ4AAIDrRmDrgnONLao936QMd2igSwH6hCCHXZNy4jQpJ05fvaVAZxua9d6hurYO3E8W7dVPFu1VRLBDE7PjNCUnTpNz3crzhDPABAAA4BoQ2LrgSO15SVJGbFiAKwH6pohgp+YUJWhOaxe65lyj1h64NMBk8e5KSVJ8RJCmtA4wmZrnVlIUd0IBAACuhMDWBUdr6yVJGXF02ICucIcH6c7iZN1ZnCxJKqurb93/VqtV+2r02tbjkqSc+DBNy4vXlFy3JmbHKoIBJgAAAB9AYOuCI3X+wJZOYAOuS1psqO4dn657x6fL57O0p+Ks1uyv0ar9lwaY2G1GI9OiNSXXrWl5bo1Mi5bTbgt06QAAAAFFYOuCI7X1igl1Mr4c6AY2m9GQ5EgNSY7UQ9Oz1dji1eYjJ7Vmf41W76vRr97dp18u3acwl91//601wOVy/w0AAAxCBLYuOFp3Xulx3F8DekKQw67JOW5NznHrn2+VTtU3ad2BWq3eX6PV+2u0dE+VJCkhMsh/9631xRPJAm8AADDwEdi6oKzugkakRgW6DGBQiA516bbhSbpteJIk//23i+Ft2Z4qvbrlmCT/Au+pufGamhenCVlxCgvinzMAADDw8BnOVViWpYozDbo1ih1sQCCkxYbqvvHpuq/1/tvuE2f8AW5fjZ7bcERPrTkkh81odHqMpub5978Vp0bJwf03AAAwAFw1sBlj0iQ9IylRkk/SE5Zl/aLDM0bSLyTdLqle0oOWZW3p/nJ736n6ZjW1+JTA8Ssg4Gw2o2EpURqWEqVHZuSoodmrTYdPtnbgqvWzJaV6fHGpIoIcmpgTp2mtAS7bHcb9NwAA0C91pcPWIumrlmVtMcZESNpsjFlsWdbuds/cJimv9WWCpN+2ftvvVZxpkCQlRhHYgL4m2GnX1Dz/TjepUHXnm7T2QI1/AuW+S/vfkqOC23a/Tcl1yx0eFNjCAQAAuuiqgc2yrBOSTrR+/6wxpkRSiqT2gW2+pGcsy7IkrTfGRBtjklp/bL/WFtjosAF9XmyYS3eMSNYdI5JlWZaO1tVr1T5/gHtnd6Ve2lwuSSpMjGjrvk3IilOIyx7gygEAADp3TXfYjDGZkkZJ2tDhXSmSytq9Xt76tn4f2CpP+wMbRyKB/sUYo4y4MGXEhen+iRny+iztPHa67f7bH9ce0f+uOiSX3abRGdFtC7yHp0TJbuP4JAAA6Bu6HNiMMeGSXpH0ZcuyznR8dyc/xOrk53hY0sOSlJ6efg1lBs7FDhuBDejf7Daj4rRoFadF6wuzcnWhyav3Dte1HZ/8yaK9+smivYoMdmhyjlvT8t2anhevtNjQQJcOAAAGsS4FNmOMU/6w9ifLsl7t5JFySWntXk+VdLzjQ5ZlPSHpCUkaO3bs/wl0fVHlmQbFhbnkcjBxDhhIQlx2zciP14z8eElSzbnGtuXdq/fXaOGuCklSZlyopuXFa1qeW5Ny4hQR7Axk2QAAYJDpypRII+n3kkosy3r8Mo+9IelRY8wL8g8bOT0Q7q9JUsXpBrprwCDgDg/S/JEpmj8yRZZl6UD1ea3aV61V+2r0ypZyPbv+iOw2o1Fp/uOT0/LdGpHC+gAAANCzutJhmyLpAUk7jDFbW9/2LUnpkmRZ1u8kvS3/SP/98o/1/1T3lxoYNeea5IlkohwwmBhjlOsJV64nXJ+akqWmFp+2HD3ZFuB+vrRUP1tSqshgR9v0SY5PAgCAntCVKZGr1fkdtfbPWJK+0F1F9SV155uUlxAe6DIABJDLYdPE7DhNzI7TP9/q/3fh4vHJlfuq9bedHJ8EAAA945qmRA5GdeebFBvqCnQZAPqQ2DCX5hUna15x8hWPT45Oj9bUXI5PAgCA60dgu4ILTV5daPYqNpzABqBz13p88mIHjuOTAACgKwhsV1BX3yRJdNgAdBnHJwEAQHcisF1B3bnWwBZGYANwfa7l+OS0vHhNzeP4JAAAuITAdgVtHTYCG4BucLXjkz9bUqrHF3N8EgAAXEJgu4K6842SCGwAekZXj09mucM0Lc+tqbkcnwQAYLAhsF1B3flmSQQ2AL3jSscnX95crmfWffD45LQ8t0akRstuu+LmFQAA0I8R2K6g7nyj7DajSL6aDaCXXevxyRn58ZqeH6/k6JBAlw4AALoRge0KTtY3KzrEKRtfvQYQYJc7PrlqX7VWlta0HZ/M9YRrel68pue7NTE7TsFOe4ArBwAAN4LAdgVnLjQrMoTuGoC+p+PxyX1V57SytForSqv13IYjemrNIQU5bBqfFdvWfcvzhMsYvgAFAEB/QmC7grMNLYoI5rcIQN9mjFF+QoTyEyL02WnZamj2asOhOq3YW62V+6r1o7dKpLdKlBQVrGl5bs3I92hqrltRoXxBCgCAvo40cgVnG5q5vwag3wl22jUjP14z8uMlScdPXdDKUn94W7izQi9uKpfNSMVp0a3HJ+M1Mo3hJQAA9EUEtis429CihMjgQJcBADckOTpE945P173j09Xi9Wlb+WmtKK3WytJq/erdffrF0n2KDHZoWuvdt+n58UqKYngJAAB9AYHtCs40NHMkEsCA4rDbNCYjRmMyYvSVm/N1qr5Jq/fX+DtwpTV6a8cJSVKeJ1zTW7t047NiGV4CAECAkEauwH+HjSORAAau6FCX7hiRrDtG+IeXlFaeazs++ez6I/r9av/wkgnZcZqe518fkMvwEgAAeg2B7TJavD7VN3m5wwZg0DDGqCAxQgWJEXpoerYuNHm14VCtVpbWaEVplX70Vol+9FaJkqOCNb118uSUHIaXAADQkwhsl3G2oUWSOBIJYNAKcdk1s8CjmQUeSUN07OLwktJqvbXjhF7YWCabkUamRbcFuOJUhpcAANCdSCOXQWADgA9KiQ7RfePTdV/b8JJTWlFaoxWl1frF0n36+ZJ9igpxamqeWzNap08mRjG4CQCAG0EauYwzDc2SxB02AOiEf3hJrMZkxOorN+fr5Pl2w0v2Veut7f7hJfkJ4ZqeF68ZBfEal8nwEgAArhWB7TIuBrbIEH6LAOBqYsJcmlecrHnFl4aXrCit0srSGj2z/oieXH1IwU6bJmTF+XfEFcQr2x3G8BIAAK6CNHIZF49EMnQEAK5N++ElD0/P0YUmr9YfqtXK0mqtKK3WDxbslhZIabEhmpnv0Yz8eE3OjVOoiw9JAAB0xEfHy+AOGwB0jxCXXbMKPJpV4JEkldXVa0VptZbvrdYrW8r17PojctltGpcVo5n5Hs0sYHUAAAAXkUYuo77JH9j4ii8AdK+02FDdPzFD90/MUGOLV5sPn9Ty0mqt2FutH79doh+/XaKU6JC2xd1TcuO4TwwAGLRII5dR3+SVJIW6uCAPAD0lyGHX5Fy3Jue69a3bi3S8dXXA8r3VWrDtuJ5/76gcNqMxGTGaWeA/PlmUFEH3DQAwaBDYLuNiYAthohkA9Jrk6BDdOz5d945PV7PXpy1HLnXfHlu4R48t3CNPRJBm5MdrZoFHU/Pcigqh+wYAGLgIbJdxoalFwU6bbCyABYCAcNptmpAdpwnZcfrG3EJVnWnw330rrdaiXRV6aXO57DajUWnRmlkQrxn5Hg1NjuTfbQDAgEJgu4z6Ji/31wCgD/FEBuvusWm6e2xa2+Lu5Xv9kyd/+k6pfvpOqdzhrra9b9Py4hUb5gp02QAA3BASyWVcaPJyfw0A+qj2i7u/ekuBas41atU+/923ZXur9Or7x2SMVJwa3Xp8Ml4jUqNlp/sGAOhnCGyXUU9gA4B+wx0epLtGpequUany+iztOHZay/dWaUVptX757j79Yuk+xYQ6NS3PH96m5cUrPiIo0GUDAHBVBLbLqG/2KoQjkQDQ79htRiPTojUyLVpfvilfJ883adX+Gi3fW6WVpdV6Y9txSdLwlKjWu2/xGpkWLYfdFuDKAQD4v0gkl1Hf2KJQJkQCQL8XE+bSncXJurM4WT6fpd0nzrR1336z/IB+9e5+RQY7NK317tuM/HglRAYHumwAACQR2C6rvsmrpChGRQPAQGKzGQ1LidKwlCg9OjtPpy80a01r921FabXe2nFCklSUFKlZBfGaVejRKLpvAIAAIrBdxoVmr0K4wwYAA1pUiFO3D0/S7cOTZFmW9lSc1YrSai3bU6UnVh7Ub5YfUGSwQ9Pz4zWrwKMZBfFyh3P3DQDQewhsl1Hf1MLQEQAYRIwxKkqKVFFSpB6ZkaMzDc1as69Gy/ZWadneai3YfkLGSCNSojSzwKNZhR6NSIli7xsAoEcR2C6DPWwAMLhFBjt12/Ak3TY86QN335btrdavWidPxoW5/GsDCj2anudWdCh73wAA3YtEchkXmjgSCQDw63j37eT5Jq3c5z86eXHvm81IYzJi/N23Ao+KkiJkDN03AMCNIbB1oqnFpxafpTACGwCgEzFhLs0fmaL5I1Pk9VnaVn5Ky/f4u28/WbRXP1m0VwmRQZpV4NHMAo+m5rkVHsSHXADAteOjRycuNHkliT1sAICrstuMRqfHaHR6jL5yS4GqzjZoxd5qLd/rnzr5wsYyOe1G4zJjNavAo1mF8cqJD6f7BgDoEhJJJ+qbWySJoSMAgGvmiQjW3WPTdPfYNDV7fdpy5KSW7a3W8r1V+vHbJfrx2yVKjQlpC2+Tst0cwQcAXBaBrRMNzT5JUrCTvTsAgOvntNs0ITtOE7Lj9C+3Fer4qQtavrday/ZW6ZUt5Xp2/RG5HDZNyo5r2/uWERcW6LIBAH0Iga0TjS3+I5HBDr7iCQDoPsnRIfr4hHR9fEK6Glu82njoZOvagCp9/83d+v6bu5XtDmtdGxCv8VmxCuJjEQAMalcNbMaYpyTdIanKsqxhnbx/pqTXJR1qfdOrlmX9oDuL7G2NrR02l4MOGwCgZwQ57Jqa59bUPLe+c8cQHa4537Y24LkNR/TUmkMKddk1JdfdOrwkXsnRIYEuGwDQy7rSYfuDpF9LeuYKz6yyLOuObqmoD2hs8Qc2vqoJAOgtme4wPejO0oNTsnShyat1B2u0bE+13t1TpcW7KyVJhYkRmlng0exCj0anR8th5wuLADDQXTWwWZa10hiT2fOl9B0Xj0QGcYcNABAAIS67ZhcmaHZhgn5gWTpQfU7L9vjvvj256qB+t+KAokKcmpEfrzlFHs3Ij2dpNwAMUN11h22SMWabpOOSvmZZ1q7OHjLGPCzpYUlKT0/vpl+6+zW1ddgIbACAwDLGKNcToVxPhB6anq2zDc1ava9GS/dUadmeKr2x7Xjb0u7ZhQmaU+RRnoe1AQAwUHRHYNsiKcOyrHPGmNslvSYpr7MHLct6QtITkjR27FirG37tHsGRSABAXxUR7NRtw5N02/Ak+XyWth87rXdLKrV0T5UeW7hHjy3co5ToEM0p8h+dnJgdp2AnH88AoL+64cBmWdaZdt9/2xjzG2OM27Ksmhv9uQPl4pFIho4AAPoym81oZFq0RqZF6yu3FKjidIOW7a3S0pIqvbSpXM+sO6IQp39wyZwij2YVeJQYFRzosgEA1+CGA5sxJlFSpWVZljFmvCSbpNobriyALk6J5EgkAKA/SYwK1n3j03Xf+HQ1NHu1/mCt3t3jD3BLSvyDS4YmR2pOoUezixI0IiVKNhtHJwGgL+vKWP/nJc2U5DbGlEv6niSnJFmW9TtJH5X0D8aYFkkXJN1rWVafPe7YFY3cYQMA9HPBTrtmFng0s8Cjf7vT0r6qc1paUqV391Tq18v265fv7pc73KWZBR7NKfRoap5bEcHOQJcNAOigK1Mi77vK+38t/9j/AaNt6Ahn/gEAA4AxRvkJEcpPiNA/zMzRyfNNWrmvWktL/CsDXt5cLqfdaHxWbOt0So+y3GGBLhsAoO6bEjmgtI31p8MGABiAYsJcmj8yRfNHpqjF69OWo6e0dE+l3i2p0g8X7NYPF+xWtjtMswv9g0vGZsZyrxsAAoTA1onGFp+MkRyc6wcADHAOu03js2I1PitW37ytSDWlSVAAACAASURBVGV19f57b3uq9My6I3py9SFFBDk0PT9eswo9mlkQL3d4UKDLBoBBg8DWicYWn4IcNnbYAAAGnbTYUH1ycqY+OTlT5xtbtGZ/jd7dU6V391TprR0nZIxUnBrdOrjEoyFJkXy8BIAeRGDrRGOzlx1sAIBBLyzIoVuGJuqWoYmyLEu7jp/xDy7ZW6X/Wlyq/1pcqqSoYM0p8mhOUYImsfMNALodga0TTV4f99cAAGjHGKNhKVEalhKlL92Up6qzDVq+t1pLSyr16pZjem79UYW67JqW59ZNRf7BJXEcnQSAG0Zg60Rjs09BTgIbAACX44kI1j1j03TP2DQ1NHu17mCtlpZUasnuKi3aVSljpNHpMZpT5NHNRQnK9YRzdBIArgOBrRONLT657AQ2AAC6Ithp16wCj2YVePTD+f6jk0tKKrWkpFL/uXCv/nPhXqXHhuqmogTdVOTRuKxYOfk4CwBdQmDrRGMLd9gAALge7Y9OfvmmfJ04fUFLS6q0tKRSz204oqfWHFJEsEMzCzy6qcijmfkeRYWysBsALofA1onGFo5EAgDQHZKiQnT/xAzdPzFD9U0tWrWvRkt2V+rdPVV6c9tx2W1G4zNj/UcnhyQoI46F3QDQHoGtExfH+gMAgO4T6nLo1qGJunVoorw+S1vLTvnvvZVU6kdvlehHb5UozxOuOUUJunmIRyPTYmRnJyqAQY7A1onGFp+iQzieAQBAT7HbjMZkxGhMRoy+PrdQR2vrtaSkUkv3VOrJVQf1uxUHFBvm0uxC/9HJaXnxCgvi0xYAgw//8nWiqcXHZWgAAHpRelyoPj01S5+emqXTF5q1srRaS0oq9c6uCr28uVwuu02TcuJ00xD/4JKkqJBAlwwAvYLA1okWr09OO0cwAAAIhKgQp+YVJ2tecbKavT5tOnyy7ejkd17bqe+8Jg1NjvQfnSxK0LCUSFYGABiwCGydaPbSYQMAoC9wtnbWJuXE6dt/V6QD1ef9RydLKvXrd/fpl0v3KSkqWDcPSdAtQxI1IZuVAQAGFgJbJ5q9lhx02AAA6FOMMcr1hCvXE65HZuSo7nyTlpZUavHuSr24qUzPrDuiiGCHZhf6J07OyI9XRDB30gH0bwS2TrT4fHLa+OocAAB9WWyYS3ePTdPdY9N0ocmr1ftrtHh3hZaUVOn1rcfb7r3dMtR/dNITGRzokgHgmhHYOtHsteR00GEDAKC/CHHZdfOQBN08JEFen6UtR0/qnV0Vemd3pb7915369l93amRatG4ekqBbhyYoJz6ce28A+gUCWyeavT456LABANAv2W1G4zJjNS4zVt+6vUj7qs5p8W7/xMmfLNqrnyzaqyx3WOu9twSNSmffG4C+i8DWiRavxZRIAAAGAGOM8hMilJ8QoS/MylXF6QYtbr339vSaQ3pi5UHFhbl0U5G/Ozc1z61gpz3QZQNAGwJbJ5gSCQDAwJQYFawHJmbogYkZOtPQrBV7q/XO7kq9veOE/rKpTCFOu6bnu3XLkETNLvQoJswV6JIBDHIEtg4sy1KLz5KDwAYAwIAWGXxp31tTi0/rD9Zq8W5/923RrsrWo5UxunlIom4ZkqC02NBAlwxgECKwddDisyRJTs6yAwAwaLgcNk3Pj9f0/Hj9YP5Q7Th2Wu/s8oe3Hy7YrR8u2K3CxAjdMtQf3oYms6wbQO8gsHXQ7PVJkpwOOmwAAAxGxhiNSI3WiNRofe3WAh2pPd86tOTSsu6U6BDNHZaoW4cmakwGQ0sA9BwCWwfNXn+HzcE/vAAAQFJGXJg+Oy1bn52WrdpzjVpaUqVFuyr07Loj+v3qQ3KHu3TzkETdOjRBk3PccvFFXwDdiMDWQcvFDht32AAAQAdx4UG6Z1ya7hmXprMNzVq+t1oLd1Xoja3H9Px7RxUR7NCcQo/mDkvU9Px4hbr4VAvAjeFfkQ7a7rAR2AAAwBVEtBta0tDs1Zr9NVq4s0JLSir12tbjCnbaND0vXnOHJWpOYYKiQp2BLhlAP0Rg66Cpxd9hc7CHDQAAdFGw0645RQmaU5SgFq9P7x2u06KdFVq0q1Lv7K6Uw2Y0KSdOtw5N1C1DE+SJCA50yQD6CQJbB5c6bAQ2AABw7Rx2mybnuDU5x63vzRuq7cdOa+HOCi3aVaF/fW2nvvP6To1Jj9GtQ/1DS9LjWBcA4PIIbB1whw0AAHQXm81oZFq0RqZF6xtzC7Sv6pwW7qzQwp0V+vHbJfrx2yUakhTZNnEyPyGcdQEAPoDA1kFTa2Bz2AhsAACg+xhjlJ8QofyECH1xTp7K6uq1aJc/vP1sSakeX1yqLHeYbh2aqLnDEjUiJUo2plYDgx6BrYMWL0ciAQBAz0uLDW1bF1B1tkGLd1dq4c4KPbnqoH634oASI4M1d1iibh+exK43YBAjsHXQ4rs4dIQOGwAA6B2eiGB9YkKGPjEhQ6frm/Xu3kr9bUeFnn/vqP6w9rDiI4I0d2iibhueqPGZsXyeAgwiBLYOmlrosAEAgMCJCnXqrlGpumtUqs41tmjZnir9becJvbS5TM+uP6K4MJduHZao24claWI24Q0Y6AhsHVzssDF0BAAABFp4kKNt11t9U4uW763W2ztO6LX3j+nPG44qJtSpW4b4O29Tct18/gIMQAS2Di7eYXNwThwAAPQhoS6Hbh+epNuHJ6mh2asVpf7w9taOE/rLpjJFhTh185AE3d4a3oIc9kCXDKAbENg6aGKsPwAA6OOCnfa2PW4NzV6t3lejt3ec0KJdFXp5c7kighy6aUiCbh+epGl5bgU7CW9Af0Vg6+DSlEgCGwAA6PuCnXbdNCRBNw1JUGOLV2v31+rtHSf0zu5K/fX9Ywpz2TWnyN95m5HvUYiL8Ab0JwS2Di5NieRIJAAA6F+CHHbNKvRoVqFH/+71ae2BWv2ttfP2xrbjCnX533/7sCTNKoxXqItPBYG+jr+lHTS1+AObiw4bAADox5x2m2bkx2tGfrx+9KFh2nCoTm/tOKFFOyv01vYTCnbaNLvQoztGJGt2oYdjk0AfRWDroMXXOnSEDhsAABggHHabpuS6NSXXrR/OH6b3DtXprR3H9bcdFXp7R4VCXXbdVJSgecXJmp7PwBKgL7lqYDPGPCXpDklVlmUN6+T9RtIvJN0uqV7Sg5ZlbenuQntLS+vQEYeNDhsAABh47DajSTlxmpQTp+/PG6r1B+u0YPtxLWw9NhkR5NDNQxM0b0SypuS65XLwOREQSF3psP1B0q8lPXOZ998mKa/1ZYKk37Z+2y81e1mcDQAABgeH3aapeW5NzXPrhx8aptX7a7Rg2wm9s7tCr245pqgQp+YOTdQdxUmalB3Hkm4gAK4a2CzLWmmMybzCI/MlPWNZliVpvTEm2hiTZFnWiW6qsVd5W49E2tnDBgAABhGn3aZZBR7NKvCosWWYVpXWaMH241qw/bj+sqlMcWEuzR2WqDtGJGt8ViyfKwG9pDvusKVIKmv3ennr2/pnYLMuLs7mK0gAAGBwCnJcWhXQ0OzV8r1VenP7Cb2ypVx/2nBUnogg3T48SXeMSNLo9BjZCG9Aj+mOwNbZ31Cr0weNeVjSw5KUnp7eDb9097vYYSOvAQAA+Pe8zR2WpLnDklTf1KKlJVVasP24/vzeUf1h7WElRwX7w1txsopTo+QfbwCgu3RHYCuXlNbu9VRJxzt70LKsJyQ9IUljx47tNNQF2sXARocNAADgg0JdDs0rTta84mSdbWjWkpJKLdh2Qn9cd1hPrj6ktNgQzRuRrPkjU1SQGBHocoEBoTsC2xuSHjXGvCD/sJHT/fX+mnRprD+dfQAAgMuLCHbqrlGpumtUqk7XN2vR7got2H5C/7PyoH6z/IAKEyM0f2SK5hUnKTUmNNDlAv1WV8b6Py9ppiS3MaZc0vckOSXJsqzfSXpb/pH+++Uf6/+pniq2N/h8luw2QzsfAACgi6JCnbpnbJruGZummnONenvHCb2+9bgeW7hHjy3co3GZMbpzZIr+bniSYsNcgS4X6Fe6MiXyvqu835L0hW6rKMBafJbshDUAAIDr4g4P0t9PytTfT8pUWV293th2XK+9f0zfeW2n/u2NXZqW59b8kSm6eUiCwoK647AXMLDxt6QDn2UxphYAAKAbpMWG6guzcvX5mTnaU3FWr289rje3HdeX/7JVwU6bbh6SqPnFyZqeH8+CbuAyCGwdtHgJbAAAAN3JGKOipEgVJUXq67cWaPPRk3p96zG9tf2E3tx2XFEhTt0+PEnzRyZrfGYsawKAdghsHdBhAwAA6Dk2m9G4zFiNy4zV9+YN1ep9NXp96zG9vvWYnn/vqBIjg3XnyGTdWZysocmRzBXAoEdg66DF5yOwAQAA9AKn3aZZhR7NKvSovqlFS0qq9MbWY3pq9SE9sfKgcj3humtUij40KkUp0SGBLhcICAJbB16fCGwAAAC9LNTl0J3F/s7ayfNNenvnCb32/jH9ZNFe/fSdvZqYFacPj07RbcOTFM6wEgwi/GnvwOvzMSUSAAAggGLCXPrEhAx9YkKGjtbW66/vH9Or75frn1/eru+8vlO3Dk3UXaNSNDXXLYedYSUY2AhsHdBhAwAA6DvS40L1pZvy9MU5udpy9JRe3VKuBdv9e97iI4I0vzhZHx6dqiHJkYEuFegRBLYOvNxhAwAA6HOMMRqTEaMxGTH67rwhWranSq9uOaY/rjusJ1cfUmFiRNt9t4TI4ECXC3QbAlsHXktyENgAAAD6rCCHXXOHJWnusCSdPN+kBduP65Utx/Qff9ujxxbu0ZRctz48OkW3Dk1UqItPd9G/8Se4A6/Px+4PAACAfiImzKUHJmXqgUmZOlh9Tq+9f0yvvn9M//SXbQp17dTcYYm6e0yaJmSx3w39E4GtA6/PosMGAADQD2XHh+srtxToyzfla9ORk3p1S7ne2n5Cr245prTYEN09Jk0fGZPKigD0KwS2Drw+SzamRAIAAPRbNpvR+KxYjc/yL+detKtCL24q0+OLS/WzJaWamuvWR8ek6tahiQp22gNdLnBFBLYOvD5LDjuBDQAAYCAIcdn1odZhJGV19Xp5c7le3lyuL72wVRHB/t1v94xN04jUKBm+aI8+iMDWQQsdNgAAgAEpLTZU/3Rzvr40J0/rD9bqpc3lemVLuf604ajyE8J195g0fWhUiuIjggJdKtCGwNaBz+IOGwAAwEBmsxlNznVrcq5b/zZ/qBZsO6GXNpfpx2+X6LGFezSr0KN7xqZpZkG8nCzmRoAR2Dpo8VpMEAIAABgkIoOd+viEdH18Qrr2V53VS5vK9er7x7R4d6Xc4S59ZHSqPjYuTdnx4YEuFYMUga0Dn2XxlRQAAIBBKNcToW/eXqR/vrVAK0qr9eKmMv1+9SH9z8qDmpgdq/vGpzOoBL2OwNZBi89SsJMOGwAAwGDlsNs0pyhBc4oSVHW2QS9vLtcL75XpSy9sVXSoUx8Znar7xqcp1xMR6FIxCBDYOvD5LNk5EgkAAABJnohgfX5mrh6ZnqO1B2r1/MajembdYf1+9SGNy4zRvePS9Xcjkui6occQ2Dpo8VmyMyUSAAAA7dhsRlPz3Jqa51bNuUa9uqVcz79Xpq++tE3/9uYu3TUqRfdNSFdhYmSgS8UAQ2DrwOtj6AgAAAAuzx0epIen5+ihadlaf7BOL2w8quffK9Mf1x3RyLRofXx8uu4oTlKoi0+1ceP4U9SB18dYfwAAAFydMUaTcuI0KSdO35vXpFe3lOuFjWX6+ivb9cO3duujY1J1/8QM5TBhEjeAwNaB16LDBgAAgGsTG+bSZ6dl6zNTs7Tx8Ek9t/6Inlt/RE+vOaypuW7dPzFDNxV55GAaOa4Rga0DOmwAAAC4XsYYjc+K1fisWFWdLdKLG8v05w1H9chzm5UYGayPT0jXvePT5IkIDnSp6CcIbB14GToCAACAbuCJCNajs/P0yIwcvbunSs+uP6LHF5fql0v3ae6wRD0wMUPjs2Jl+NwTV0Bg64Cx/gAAAOhODrtNtwxN1C1DE3Ww+pz+tOGoXtpUpgXbT6ggIUL3T8rQXaNSFB7Ep+b4vzhE20ELgQ0AAAA9JDs+XN+5Y4g2fOsmPfaR4XLYjb7z2k5N/Pel+sGbu3W0tj7QJaKPIcZ34LMIbAAAAOhZIS67PjYuXfeMTdP7Zaf0x7WH9cy6w3p67SHdVJSgT0/J0sRsjkuCwPZ/0GEDAABAbzHGaHR6jEanx+hbtxfp2XVH9Of3jmrx7koVJkbo01OzdGdxsoKd9kCXigDhSGQHXgIbAAAAAiAhMlhfu7VAa/9lth77yHBZlvT1l7dryv97V48vLlXV2YZAl4gAoMPWAVMiAQAAEEjBzkvHJdcdqNVTaw7pV+/u02+X79e8Ecn69NQsDUuJCnSZ6CUEtg68Pkt2O4ENAAAAgWWM0eRctybnunWo5rz+uPawXtpUplffP6aJ2bH63IwczcyP557bAMeRyA7osAEAAKCvyXKH6ft3DtW6b83Rt28v0pHaen3q6Y267Rer9MrmcjW1+AJdInoIga0DL1MiAQAA0EdFBjv10PRsrfjnWXr8nmJZlvTVl7Zpxk+W6clVB3WusSXQJaKbEdjasSxLliXZ6LABAACgD3M5bPrw6FQt/PI0Pf2pccqIC9WP3irRpP9YqscW7lHVGQaUDBTcYWvHsvzfEtgAAADQHxhjNKvAo1kFHm0rO6UnVh7U/6w4oN+vOqQPj07Rw9OzlR0fHugycQMIbO34WhMbJyIBAADQ3xSnReu/PzFah2vO68nVB/XSpnK9uKlMd4xI1hdm5aogMSLQJeI6cCSyHe/FwEZiAwAAQD+V6Q7Tjz40XKu/MVsPTc/WkpJK3frzlXrk2c3aeex0oMvDNaLD1g5HIgEAADBQxEcE6Zu3FemR6Tl6es0hPb32sBbuqtDsQo8enZ2r0ekxgS4RXUCHrR2ORAIAAGCgiQlz6Su3FGj1N2bra7fka8vRk/rwb9bq/ic3aP3B2kCXh6sgsLXjo8MGAACAASoqxKlHZ+dpzTdm65u3FWpPxRnd+8R63fvEOm06XBfo8nAZXQpsxpi5xpi9xpj9xph/6eT9Dxpjqo0xW1tfPtv9pfa8ix028hoAAAAGqrAghz43I0ervzFb371jiPZXnddHf7dOn3zqPe0o545bX3PVwGaMsUv6b0m3SRoi6T5jzJBOHv2LZVkjW1+e7OY6e4XVuiCeDhsAAAAGumCnXZ+emqWVX5+pf7mtUNvKT2ner1frc89u0t6Ks4EuD6260mEbL2m/ZVkHLctqkvSCpPk9W1ZgeLnDBgAAgEEm1OXQIzNytOrrs/Tlm/K0dn+t5v5ipb70wvs6VHM+0OUNel0JbCmSytq9Xt76to4+YozZbox52RiT1i3V9bKLRyLtJDYAAAAMMhHBTn35pnyt/PosPTIjR+/sqtRNj6/QN1/drqozDYEub9DqSmDrLL1YHV5/U1KmZVkjJC2R9MdOfyJjHjbGbDLGbKqurr62SnvBpTtsBDYAAAAMTjFhLn1jbqFWfn2WHpiYoZc3l2vmT5frZ4tLdb6xJdDlDTpdCWzlktp3zFIlHW//gGVZtZZlNba++r+SxnT2E1mW9YRlWWMtyxobHx9/PfX2KPawAQAAAH7xEUH6/p1DteQrMzSrwKNfLN2nmT9drj9vOKoWry/Q5Q0aXQlsGyXlGWOyjDEuSfdKeqP9A8aYpHav3imppPtK7D3sYQMAAAA+KCMuTP/9idF69fOTlREbqm/9dYfm/mKVlpZUyrI6HrxDd7tqYLMsq0XSo5IWyR/EXrQsa5cx5gfGmDtbH/uiMWaXMWabpC9KerCnCu5J7GEDAAAAOjc6PUYvPTJJv7t/jLw+S5/54yZ9/H83MFGyhzm68pBlWW9LervD277b7vvflPTN7i2t9/l87GEDAAAALscYo7nDEjWnyKM/bziqxxeX6vZfrtIDEzP0TzfnKyrEGegSB5wuLc4eLJgSCQAAAFyd027TJydnatnXZupj49L0x3WHNfuny/XixrK2Jgi6B4GtHY5EAgAAAF0XG+bSv981XG8+OlWZ7jB9/ZXtuus3a7S17FSgSxswCGztXBrrH+BCAAAAgH5kWEqUXn5kkh6/p1jHTzfort+s0Xde26kzDc2BLq3fI7C1Y7VNiSSxAQAAANfCGKMPj07Vsq/N1IOTM/WnDUd08+MrtGhXRaBL69cIbO1wJBIAAAC4MeFBDn1v3lD99fNTFBPq0uee3azPPbtJFacbAl1av0Rga8frYw8bAAAA0B2K06L15j9O1TfmFmr53mrd/PgKPbf+CLvbrhGBrZ22xdkkNgAAAOCGOe02/cPMHL3zT9M1Ii1K//raTn3y6Y10264Bga0diyORAAAAQLfLiAvTc5+ZoB9+aJg2HqrTLT9bode3HqPb1gUEtnbaOmzkNQAAAKBbGWP0wMQM/e1L05TrCdeXXtiqR//8vurONwW6tD6NwNYOQ0cAAACAnpXpDtNLj0zW1+cW6J3dFZr785Vaf7A20GX1WQS2dtjDBgAAAPQ8u83o8zNz9doXpig8yKGP/+96/XLpvrYhgLiEwNaOz8ceNgAAAKC3DE2O0hv/OFXzipP1+OJS/f1TG1R9tjHQZfUpBLZ2LgZ6O5fYAAAAgF4RHuTQzz82Uv/vw8O16fBJ3f7LVdp0uC7QZfUZBLZ2OBIJAAAA9D5jjO4dn67XH52iMJdd9/3ver24sSzQZfUJBLZ2Lk2JJLEBAAAAva0wMVKvfWGKJmTF6euvbNf339ilFq8v0GUFFIGtHfawAQAAAIEVHerSHz41Tp+ekqU/rD2sB5/eqDMNzYEuK2AIbO2whw0AAAAIPIfdpu/OG6L//OgIrT9Yq3v/Z72qzjYEuqyAILC1c3GMqKHDBgAAAATcPWPT9PsHx+lw7Xl95LdrdajmfKBL6nUEtnYspkQCAAAAfcqM/Hg9/9BEnW/06qO/Xaudx04HuqReRWBrhyORAAAAQN9TnBatlx+ZpGCnXff/foN2HR88oY3A1o6PoSMAAABAn5QdH67nH5qoUKddn3hyg3YfPxPoknoFga0d9rABAAAAfVd6XKheeHiSQpx2feLJ9dpfdTbQJfU4Als7FnvYAAAAgD7NH9omym6z6ZNPbVTVmYE9PZLA1s7FnXwENgAAAKDvyogL09MPjtPJ+iZ96g8bda6xJdAl9RgCWzsXj0Ta+V0BAAAA+rThqVH670+M1p6Ks/rS8+/Ld3EgxQBDNGnn0h02OmwAAABAXzerwKPv3jFES/dU6TfL9we6nB5BYGvHYkokAAAA0K/8/aQMzR+ZrP9aXKpV+6oDXU63I7C1wx42AAAAoH8xxug/PjxceZ5wfeXFbTpV3xTokroVga0d9rABAAAA/U+oy6GffWykTp5v0g/e3B3ocroVga2dixcVyWsAAABA/zI0OUqfn5WrV98/pmV7qgJdTrchsLVzaUokiQ0AAADobx6dlauc+DD9cMFuNV/c2dXPEdja4UgkAAAA0H+5HDZ9+++KdLDmvJ5bfyTQ5XQLAls7l8b6B7gQAAAAANdlVoFHU3Lj9Kt39+tCkzfQ5dwwAls7VtuUSBIbAAAA0B8ZY/TF2XmqO9+kl7eUB7qcG0Zga4cjkQAAAED/Nz4rViPTovXkqoNtgwX7KwJbO14fe9gAAACA/s4Yo09OztCR2nptPnoy0OXcEAJbO22Ls0lsAAAAQL92y5BEhTjtenXLsUCXckMIbO1YHIkEAAAABoSwIIdmF3m0pKSybVZFf0Rga6etw0ZeAwAAAPq9abluVZ9t1P6qc4Eu5boR2Nph6AgAAAD+f3v3FitXVcdx/PtLK0W5X6zcL0FjBKOIDUWIpglagRCqBrTES70FUTHyYAJqgqS+gIpGjdGgkCBBwBvaB1CbaOITSCEoIirVFKnUcikBEQ2p/H2YXRzHmdNpO5fdM99P0pyZvdf0/M/8z9pr/c9ee4/mj1OPOxiAX23YMuVIdp4FWxc/h02SJEmaP4488IXsvWghf/zb36cdyk4bqmBLckaSPyRZn+TSPvsXJbm52X9HkmNGHegkPPecn8MmSZIkzRdJOG7x3qx/dB4viUyyAPgacCZwPHB+kuN7mn0AeKKqXgp8Cbhy1IFOwrYlkQss2CRJkqR5YfE+i3j86WenHcZOG+YM28nA+qr6c1U9C9wErOhpswK4rnn8feD0ZPerelwSKUmSJM0vB+21B4//Y/ct2BYO0eZw4KGu5xuBpYPaVNXWJE8CBwGPdTdKcgFwAcBRRx21kyGPT1WRdE6dSpIkSdr9rTr1GM458bBph7HThinY+lUvvR9kMEwbqupq4GqAJUuWtO7DEN51ytEsP+GQaYchSZIkaURecei+0w5hlwxTsG0Ejux6fgTw8IA2G5MsBPYDdrt7Zy7ed08W77vntMOQJEmSJGC4a9juBF6W5NgkewArgTU9bdYAq5rH5wI/r93548QlSZIkqQW2e4atuSbtIuCnwALg2qq6L8lqYF1VrQGuAa5Psp7OmbWV4wxakiRJkmbBMEsiqapbgVt7tl3W9fhfwHmjDU2SJEmSZttQH5wtSZIkSZo8CzZJkiRJaikLNkmSJElqKQs2SZIkSWopCzZJkiRJaikLNkmSJElqKQs2SZIkSWqpVNV0vnHyKPDgVL55fwcDj007CD3PfLSHuWgPc9Ee5qI9zEV7mIv2MBftsb1cHF1VL57rP5hawdY2SdZV1ZJpx6EO89Ee5qI9zEV7mIv2MBftYS7aw1y0xyhy4ZJISZIkSWopCzZJkiRJaikLtv+6etoB6H+Yj/YwF+1hLtrDXLSHuWgPc9Ee5qI9djkXXsMmSZIkSS3lGTZJkiRJaqmZK9iSnJHkD0nWJ7m0z/5FSW5u9t+R5JjJRzn/JTkyyS+S3J/kviQf79NmWZInk9zTrpwDHQAABVpJREFU/LtsGrHOgiQbktzbvM/r+uxPkq80/eI3SU6aRpzzXZKXd/2+35PkqSQX97SxX4xRkmuTPJLkt13bDkyyNskDzdcDBrx2VdPmgSSrJhf1/DQgF59P8vvmOHRLkv0HvHbOY5p2zIBcXJ7kr13HorMGvHbOeZd2zIBc3NyVhw1J7hnwWvvFiAyax45rvJipJZFJFgB/BN4EbATuBM6vqt91tfkI8KqqujDJSuCtVfWOqQQ8jyU5FDi0qu5Osg9wF/CWnlwsAz5RVWdPKcyZkWQDsKSq+n5OSDMQfww4C1gKfLmqlk4uwtnTHK/+Ciytqge7ti/DfjE2Sd4APA18u6pe2Wz7HLClqq5oJpwHVNUlPa87EFgHLAGKzjHttVX1xER/gHlkQC6WAz+vqq1JrgTozUXTbgNzHNO0Ywbk4nLg6ar6whyv2+68SzumXy569l8FPFlVq/vs24D9YiQGzWOB9zKG8WLWzrCdDKyvqj9X1bPATcCKnjYrgOuax98HTk+SCcY4E6pqU1Xd3Tz+O3A/cPh0o9IcVtAZHKqqbgf2bw5WGp/TgT91F2sav6r6JbClZ3P3uHAdnUG515uBtVW1pRl01wJnjC3QGdAvF1X1s6ra2jy9HThi4oHNoAH9YhjDzLu0A+bKRTNffTtw40SDmkFzzGPHMl7MWsF2OPBQ1/ON/H+R8HybZlB4EjhoItHNqHSWnb4GuKPP7tcl+XWS25KcMNHAZksBP0tyV5IL+uwfpu9otFYyeNC1X0zWS6pqE3QGaWBxnzb2kcl7P3DbgH3bO6ZpNC5qlqdeO2Dpl/1isl4PbK6qBwbst1+MQc88dizjxawVbP3OlPWuCR2mjUYkyd7AD4CLq+qpnt13A0dX1auBrwI/mnR8M+S0qjoJOBP4aLPkopv9YoKS7AGcA3yvz277RTvZRyYoyaeBrcANA5ps75imXfd14DjgRGATcFWfNvaLyTqfuc+u2S9GbDvz2IEv67Ntzn4xawXbRuDIrudHAA8PapNkIbAfO7cMQNuR5AV0fslvqKof9u6vqqeq6unm8a3AC5IcPOEwZ0JVPdx8fQS4hc4ylm7D9B2NzpnA3VW1uXeH/WIqNm9bAtx8faRPG/vIhDQX6J8NvLMGXIg/xDFNu6iqNlfVv6vqOeCb9H+P7RcT0sxZ3wbcPKiN/WK0BsxjxzJezFrBdifwsiTHNn/BXgms6WmzBth2t5Zz6Vzc7F+DRqxZZ30NcH9VfXFAm0O2XT+Y5GQ6v6+PTy7K2ZBkr+aCWZLsBSwHftvTbA3wnnScQueC5k0TDnWWDPwrqf1iKrrHhVXAj/u0+SmwPMkBzdKw5c02jVCSM4BLgHOq6pkBbYY5pmkX9VzH/Fb6v8fDzLs0Gm8Efl9VG/vttF+M1hzz2LGMFwt3PeTdR3NXqYvovCkLgGur6r4kq4F1VbWGzpt/fZL1dM6srZxexPPaacC7gXu7bj/7KeAogKr6Bp2C+cNJtgL/BFZaPI/FS4BbmhpgIfCdqvpJkgvh+VzcSucOkeuBZ4D3TSnWeS/Ji+jcUe1DXdu6c2G/GKMkNwLLgIOTbAQ+A1wBfDfJB4C/AOc1bZcAF1bVB6tqS5LP0pmgAqyuKldn7IIBufgksAhY2xyzbm/u6nwY8K2qOosBx7Qp/AjzxoBcLEtyIp2lXBtojlnduRg075rCjzBv9MtFVV1Dn+ue7RdjNWgeO5bxYqZu6y9JkiRJu5NZWxIpSZIkSbsNCzZJkiRJaikLNkmSJElqKQs2SZIkSWopCzZJkiRJaikLNkmSJElqKQs2SZIkSWopCzZJkiRJaqn/AOMAKWgZDwyqAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 1080x360 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"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(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)\n",
"\n",
"val=np.sqrt(lambda_function(mB_mass, mKst_mass, q2_range))*beta_l(q2_range, mmu_mass)\n",
"plt.plot(q2_range/1e6, val);\n",
"fig=plt.gcf()\n",
"fig.set_size_inches(15,5)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA3MAAAEvCAYAAADvmpjfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAcB0lEQVR4nO3df7Cl9V0f8PfHXSE1bUkkq1Vgs7SgMxtN0+QKdvzRjFRc1LraQrPoKG3prKnSqe3YkbQjjVRnwF9MWmntNuAQtILFxt4ZVzEVa1tHKUuMJptIvdmusiGTEJYhxYi48dM/zsG5nty7+9zsuXvvs+f1mtnZ58f3ufd79uE557z5/qruDgAAAOPyGVtdAQAAADZOmAMAABghYQ4AAGCEhDkAAIAREuYAAABGSJgDAAAYoZ1bXYFZr3rVq3rPnj1bXQ0AAIAt8fjjj3+su3edqdy2C3N79uzJkSNHtroaAAAAW6Kqfm9IOd0sAQAARkiYAwAAGCFhDgAAYISEOQAAgBES5gAAAEZImAMAABghYQ4AAGCEhDkAAIAREuYAAABGSJgDAAAYIWEOAABghHZudQUYbs+tPz+47PE7vm4TawIAAGw1YW6LbSSgAQAAvESY2wTbIaBtVh20+AEAwPZgzBwAAMAICXMAAAAjpJslGzK0+6bumAAAsLm0zAEAAIzQoDBXVfuq6omqWqmqW9c4f2FVPTg9/2hV7Zke/8yquq+q3ltVH6iqt8y3+gAAAIvpjGGuqnYkuTvJdUn2JrmxqvbOFLs5ybPdfUWSu5LcOT1+Q5ILu/uLk7whybe/FPQAAAD49A1pmbsqyUp3H+vuF5M8kGT/TJn9Se6bbj+U5JqqqiSd5OVVtTPJn0vyYpKPz6XmAAAAC2zIBCiXJHly1f6JJFevV6a7T1XVc0kuziTY7U/y4SSfleSfdvfJ2V9QVQeTHEyS3bt3b/AlsB1tZJ07k6UAAMDGDQlztcaxHljmqiSfTPL5SV6Z5H9W1X/r7mN/pmD3oSSHkmRpaWn2Z3OeE/wAAGDjhnSzPJHkslX7lyZ5ar0y0y6VFyU5meSbk/xid/9xd380ya8lWTrbSgMAACy6IWHusSRXVtXlVXVBkgNJlmfKLCe5abp9fZJHuruT/H6Sr6qJlyf50iS/M5+qAwAALK4zdrOcjoG7JcnDSXYkube7j1bV7UmOdPdyknuS3F9VK5m0yB2YXn53kp9I8r5MumL+RHf/9ia8DhaELpkAADAxZMxcuvtwksMzx25btf1CJssQzF73/FrHAQAAODuDFg0HAABgexHmAAAARmhQN0sYI+PrAAA4n2mZAwAAGCEtcxCteAAAjI+WOQAAgBES5gAAAEZImAMAABghY+ZggzYyvm4jjMUDAGAjtMwBAACMkDAHAAAwQsIcAADACBkzB9vE0LF4xtYBAJBomQMAABglLXMwMhuZTVMrHgDA+UvLHAAAwAhpmYPzmFY8AIDzl5Y5AACAERLmAAAARkg3y4E20l0NAABgs2mZAwAAGCFhDgAAYIQGdbOsqn1J3pZkR5K3d/cdM+cvTPKOJG9I8kySN3X38ar6liT/fFXR1yZ5fXe/Zx6VB+bHzJcAAONyxpa5qtqR5O4k1yXZm+TGqto7U+zmJM929xVJ7kpyZ5J090919+u6+3VJvjXJcUEOAADg7A3pZnlVkpXuPtbdLyZ5IMn+mTL7k9w33X4oyTVVVTNlbkzy02dTWQAAACaGhLlLkjy5av/E9NiaZbr7VJLnklw8U+ZNWSfMVdXBqjpSVUeefvrpIfUGAABYaEPGzM22sCVJb6RMVV2d5BPd/b61fkF3H0pyKEmWlpZmfzawzWzWUh3G4gEADDekZe5EkstW7V+a5Kn1ylTVziQXJTm56vyB6GIJAAAwN0PC3GNJrqyqy6vqgkyC2fJMmeUkN023r0/ySHd3klTVZyS5IZOxdgAAAMzBGbtZdvepqrolycOZLE1wb3cfrarbkxzp7uUk9yS5v6pWMmmRO7DqR3xlkhPdfWz+1QfOJ5ZHAAAYbtA6c919OMnhmWO3rdp+IZPWt7Wu/e9JvvTTryIAAACzhnSzBAAAYJsZ1DIHsN3okgkALDotcwAAACMkzAEAAIyQbpbAeW8zFjnXdRMA2Gpa5gAAAEZImAMAABghYQ4AAGCEhDkAAIAREuYAAABGSJgDAAAYIWEOAABghKwzB/Bp2MjaddakAwA2g5Y5AACAERLmAAAARkg3S4BNpksmALAZtMwBAACMkDAHAAAwQsIcAADACAlzAAAAI2QCFIBtxGQpAMBQwhzASAl+ALDYBnWzrKp9VfVEVa1U1a1rnL+wqh6cnn+0qvasOvfaqvr1qjpaVe+tqpfNr/oAAACL6Yxhrqp2JLk7yXVJ9ia5sar2zhS7Ocmz3X1FkruS3Dm9dmeSn0zy5u5+TZI3JvnjudUeAABgQQ1pmbsqyUp3H+vuF5M8kGT/TJn9Se6bbj+U5JqqqiTXJvnt7v6tJOnuZ7r7k/OpOgAAwOIaMmbukiRPrto/keTq9cp096mqei7JxUm+IElX1cNJdiV5oLt/8KxrDcCGGF8HAOefIWGu1jjWA8vsTPLlSb4kySeS/HJVPd7dv/xnLq46mORgkuzevXtAlQDYLEODn9AHAFtrSDfLE0kuW7V/aZKn1iszHSd3UZKT0+O/2t0f6+5PJDmc5PWzv6C7D3X3Uncv7dq1a+OvAgAAYMEMCXOPJbmyqi6vqguSHEiyPFNmOclN0+3rkzzS3Z3k4SSvrarPmoa8v5Hk/fOpOgAAwOI6YzfL6Ri4WzIJZjuS3NvdR6vq9iRHuns5yT1J7q+qlUxa5A5Mr322qn40k0DYSQ539/CBGwBsW8bhAcDWGrRoeHcfzqSL5Opjt63afiHJDetc+5OZLE8AAADAnAxaNBwAAIDtZVDLHACcDV0yAWD+tMwBAACMkJY5ALYVrXgAMIyWOQAAgBES5gAAAEZImAMAABghYQ4AAGCEhDkAAIARMpslAKNl5ksAFpmWOQAAgBES5gAAAEZImAMAABghY+YAWAjG1wFwvtEyBwAAMELCHAAAwAjpZgkAM3TJBGAMhDkAOAuCHwBbRTdLAACAERLmAAAARkiYAwAAGCFj5gDgHBk6vs7YOgCGGNQyV1X7quqJqlqpqlvXOH9hVT04Pf9oVe2ZHt9TVX9YVe+Z/vnx+VYfAABgMZ2xZa6qdiS5O8lXJzmR5LGqWu7u968qdnOSZ7v7iqo6kOTOJG+anvtgd79uzvUGAABYaEO6WV6VZKW7jyVJVT2QZH+S1WFuf5K3TrcfSvJjVVVzrCcALAzLHQAwxJBulpckeXLV/onpsTXLdPepJM8luXh67vKq+s2q+tWq+oqzrC8AAAAZ1jK3VgtbDyzz4SS7u/uZqnpDkp+rqtd098f/zMVVB5McTJLdu3cPqBIAkGjFA1hkQ1rmTiS5bNX+pUmeWq9MVe1MclGSk939R939TJJ09+NJPpjkC2Z/QXcf6u6l7l7atWvXxl8FAADAghkS5h5LcmVVXV5VFyQ5kGR5psxykpum29cneaS7u6p2TSdQSVX95SRXJjk2n6oDAAAsrjN2s+zuU1V1S5KHk+xIcm93H62q25Mc6e7lJPckub+qVpKczCTwJclXJrm9qk4l+WSSN3f3yc14IQAAAItk0KLh3X04yeGZY7et2n4hyQ1rXPezSX72LOsIAADAjEGLhgMAALC9CHMAAAAjNKibJQAwfhtZxmAjLHkAsDW0zAEAAIyQMAcAADBCulkCAGdlI903dckEmB8tcwAAACMkzAEAAIyQMAcAADBCwhwAAMAICXMAAAAjZDZLAOCcMfMlwPwIcwDAtiT4AZyebpYAAAAjJMwBAACMkG6WAMDoDe2SqTsmcD7RMgcAADBCwhwAAMAICXMAAAAjJMwBAACMkDAHAAAwQmazBAAWhoXIgfPJoDBXVfuSvC3JjiRv7+47Zs5fmOQdSd6Q5Jkkb+ru46vO707y/iRv7e4fnk/VAQA2j+AHbHdn7GZZVTuS3J3kuiR7k9xYVXtnit2c5NnuviLJXUnunDl/V5JfOPvqAgAAkAwbM3dVkpXuPtbdLyZ5IMn+mTL7k9w33X4oyTVVVUlSVd+Y5FiSo/OpMgAAAEPC3CVJnly1f2J6bM0y3X0qyXNJLq6qlyf5niTfd/ZVBQAA4CVDwlytcawHlvm+JHd19/On/QVVB6vqSFUdefrppwdUCQAAYLENmQDlRJLLVu1fmuSpdcqcqKqdSS5KcjLJ1Umur6ofTPKKJH9SVS9094+tvri7DyU5lCRLS0uzQREAAIAZQ8LcY0murKrLk3woyYEk3zxTZjnJTUl+Pcn1SR7p7k7yFS8VqKq3Jnl+NsgBAACwcWcMc919qqpuSfJwJksT3NvdR6vq9iRHuns5yT1J7q+qlUxa5A5sZqUBALYTyxgAW6EmDWjbx9LSUh85cmSrq/EpNvImDQAwD4IfLKaqery7l85UbsgEKAAAAGwzwhwAAMAIDZkABQCALWAsHnA6WuYAAABGSJgDAAAYIWEOAABghIyZAwA4DxhfB4tHyxwAAMAICXMAAAAjJMwBAACMkDAHAAAwQiZAAQBYMEMnSzFRCmxvwhwAAGsyQyZsb7pZAgAAjJAwBwAAMELCHAAAwAgJcwAAACMkzAEAAIyQ2SwBADhrZr6Ec0/LHAAAwAgJcwAAACMkzAEAAIzQoDFzVbUvyduS7Ejy9u6+Y+b8hUnekeQNSZ5J8qbuPl5VVyU59FKxJG/t7nfOq/IAAIyP8XUwH2dsmauqHUnuTnJdkr1JbqyqvTPFbk7ybHdfkeSuJHdOj78vyVJ3vy7JviT/oapMugIAAHCWhnSzvCrJSncf6+4XkzyQZP9Mmf1J7ptuP5Tkmqqq7v5Ed5+aHn9Zkp5HpQEAABbdkFayS5I8uWr/RJKr1yvT3aeq6rkkFyf5WFVdneTeJK9O8q2rwh0AAJyWLpmwviEtc7XGsdkWtnXLdPej3f2aJF+S5C1V9bJP+QVVB6vqSFUdefrppwdUCQAAYLENCXMnkly2av/SJE+tV2Y6Ju6iJCdXF+juDyT5gyRfNPsLuvtQdy9199KuXbuG1x4AAGBBDQlzjyW5sqour6oLkhxIsjxTZjnJTdPt65M80t09vWZnklTVq5N8YZLjc6k5AADAAjvjmLnpGLhbkjycydIE93b30aq6PcmR7l5Ock+S+6tqJZMWuQPTy788ya1V9cdJ/iTJd3T3xzbjhQAAsNiMr2PRDFomoLsPJzk8c+y2VdsvJLlhjevuT3L/WdYRAACAGUO6WQIAALDNCHMAAAAjJMwBAACM0KAxcwAAcD4xWQrnAy1zAAAAIyTMAQAAjJAwBwAAMELCHAAAwAiZAAUAAE7DZClsV1rmAAAARkiYAwAAGCFhDgAAYISMmQMAgDkZOr7O2DrmQcscAADACAlzAAAAIyTMAQAAjJAxcwAAcI5Zu4550DIHAAAwQsIcAADACAlzAAAAIyTMAQAAjJAJUAAAYBszWQrrGdQyV1X7quqJqlqpqlvXOH9hVT04Pf9oVe2ZHv/qqnq8qt47/fur5lt9AACAxXTGMFdVO5LcneS6JHuT3FhVe2eK3Zzk2e6+IsldSe6cHv9Ykr/V3V+c5KYk98+r4gAAAItsSMvcVUlWuvtYd7+Y5IEk+2fK7E9y33T7oSTXVFV1929291PT40eTvKyqLpxHxQEAABbZkDB3SZInV+2fmB5bs0x3n0ryXJKLZ8r8nSS/2d1/9OlVFQAAgJcMmQCl1jjWGylTVa/JpOvltWv+gqqDSQ4mye7duwdUCQAAmGWylMUypGXuRJLLVu1fmuSp9cpU1c4kFyU5Od2/NMk7k3xbd39wrV/Q3Ye6e6m7l3bt2rWxVwAAALCAhoS5x5JcWVWXV9UFSQ4kWZ4ps5zJBCdJcn2SR7q7q+oVSX4+yVu6+9fmVWkAAIBFd8YwNx0Dd0uSh5N8IMnPdPfRqrq9qr5hWuyeJBdX1UqSf5bkpeULbklyRZLvrar3TP98ztxfBQAAwIIZtGh4dx9Ocnjm2G2rtl9IcsMa131/ku8/yzoCAAAwY9Ci4QAAAGwvg1rmAACA84uZL8dPyxwAAMAICXMAAAAjJMwBAACMkDAHAAAwQsIcAADACJnNEgAAOC0zX25PWuYAAABGSJgDAAAYIWEOAABghIQ5AACAETIBCgAAMDcmSzl3tMwBAACMkDAHAAAwQsIcAADACAlzAAAAIyTMAQAAjJAwBwAAMELCHAAAwAhZZw4AANgS1qQ7O1rmAAAARmhQmKuqfVX1RFWtVNWta5y/sKoenJ5/tKr2TI9fXFW/UlXPV9WPzbfqAAAAi+uMYa6qdiS5O8l1SfYmubGq9s4UuznJs919RZK7ktw5Pf5Cku9N8t1zqzEAAACDWuauSrLS3ce6+8UkDyTZP1Nmf5L7ptsPJbmmqqq7/6C7/1cmoQ4AAIA5GRLmLkny5Kr9E9Nja5bp7lNJnkty8TwqCAAAwKcaEuZqjWP9aZRZ/xdUHayqI1V15Omnnx56GQAAwMIaEuZOJLls1f6lSZ5ar0xV7UxyUZKTQyvR3Ye6e6m7l3bt2jX0MgAAgIU1ZJ25x5JcWVWXJ/lQkgNJvnmmzHKSm5L8epLrkzzS3YNb5gAAAE5n6Jp0i7Qe3RnDXHefqqpbkjycZEeSe7v7aFXdnuRIdy8nuSfJ/VW1kkmL3IGXrq+q40n+YpILquobk1zb3e+f/0sBAABYHENa5tLdh5Mcnjl226rtF5LcsM61e86ifgAAAKxh0KLhAAAAbC/CHAAAwAgJcwAAACMkzAEAAIyQMAcAADBCwhwAAMAICXMAAAAjJMwBAACMkDAHAAAwQsIcAADACO3c6goAAADMy55bf35w2eN3fN0m1mTzaZkDAAAYIWEOAABghIQ5AACAERLmAAAARkiYAwAAGCFhDgAAYISEOQAAgBES5gAAAEZImAMAABghYQ4AAGCEhDkAAIARGhTmqmpfVT1RVStVdesa5y+sqgen5x+tqj2rzr1levyJqvqa+VUdAABgcZ0xzFXVjiR3J7kuyd4kN1bV3pliNyd5truvSHJXkjun1+5NciDJa5LsS/Lvpj8PAACAszCkZe6qJCvdfay7X0zyQJL9M2X2J7lvuv1QkmuqqqbHH+juP+ru/5tkZfrzAAAAOAtDwtwlSZ5ctX9iemzNMt19KslzSS4eeC0AAAAbtHNAmVrjWA8sM+TaVNXBJAenu89X1RMD6nUuvSrJx7a6EiRxL7YT92L7cC+2D/di+3Avtg/3YvtwL2bUnVv2q890L1495IcMCXMnkly2av/SJE+tU+ZEVe1MclGSkwOvTXcfSnJoSIW3QlUd6e6lra4H7sV24l5sH+7F9uFebB/uxfbhXmwf7sX2Ma97MaSb5WNJrqyqy6vqgkwmNFmeKbOc5Kbp9vVJHununh4/MJ3t8vIkVyb532dbaQAAgEV3xpa57j5VVbckeTjJjiT3dvfRqro9yZHuXk5yT5L7q2olkxa5A9Nrj1bVzyR5f5JTSb6zuz+5Sa8FAABgYQzpZpnuPpzk8Myx21Ztv5DkhnWu/YEkP3AWddwOtm0X0AXkXmwf7sX24V5sH+7F9uFebB/uxfbhXmwfc7kXNekNCQAAwJgMGTMHAADANiPMTVXVvqp6oqpWqurWNc5fWFUPTs8/WlV7zn0tz39VdVlV/UpVfaCqjlbVP1mjzBur6rmqes/0z21r/Szmo6qOV9V7p//WR9Y4X1X1b6bPxm9X1eu3op7nu6r6wlX/zb+nqj5eVd81U8azsUmq6t6q+mhVvW/Vsc+uqndV1e9O/37lOtfeNC3zu1V101plGG6de/FDVfU70/egd1bVK9a59rTvZ2zMOvfirVX1oVXvQ1+7zrWn/d7FxqxzLx5cdR+OV9V71rnWczEn632P3czPC90sk1TVjiT/J8lXZ7KcwmNJbuzu968q8x1JXtvdb66qA0m+qbvftCUVPo9V1ecl+bzufndV/YUkjyf5xpl78cYk393dX79F1VwoVXU8yVJ3r7kWyvSD+h8n+dokVyd5W3dffe5quHim71kfSnJ1d//equNvjGdjU1TVVyZ5Psk7uvuLpsd+MMnJ7r5j+mX0ld39PTPXfXaSI0mWMlln9fEkb+juZ8/pCziPrHMvrs1kJu1TVZNVo2bvxbTc8Zzm/YyNWedevDXJ8939w6e57ozfu9iYte7FzPkfSfJcd9++xrnj8VzMxXrfY5P8vWzS54WWuYmrkqx097HufjHJA0n2z5TZn+S+6fZDSa6pqrUWRecsdPeHu/vd0+3/l+QDSS7Z2lpxBvsz+fDo7v6NJK+Yvpmxea5J8sHVQY7N1d3/I5PZmldb/blwXyYf2LO+Jsm7uvvk9AP5XUn2bVpFF8Ba96K7f6m7T013fyOTdW3ZZOs8F0MM+d7FBpzuXky/r/7dJD99Tiu1gE7zPXbTPi+EuYlLkjy5av9EPjVA/GmZ6QfGc0kuPie1W1A16cr615I8usbpv15Vv1VVv1BVrzmnFVs8neSXqurxqjq4xvkhzw/zdSDrfyh7Ns6dz+3uDyeTD/Akn7NGGc/HufcPkvzCOufO9H7GfNwy7fJ67zrdyTwX59ZXJPlId//uOuc9F5tg5nvspn1eCHMTa7WwzfY/HVKGOamqP5/kZ5N8V3d/fOb0u5O8urv/apJ/m+TnznX9FsyXdffrk1yX5DunXTlW82ycQ1V1QZJvSPKf1zjt2dh+PB/nUFX9y0zWtf2pdYqc6f2Ms/fvk/yVJK9L8uEkP7JGGc/FuXVjTt8q57mYszN8j133sjWOnfG5EOYmTiS5bNX+pUmeWq9MVe1MclE+va4FnEFVfWYmD8BPdfd/mT3f3R/v7uen24eTfGZVveocV3NhdPdT078/muSdmXSPWW3I88P8XJfk3d39kdkTno1z7iMvdSme/v3RNcp4Ps6R6WQBX5/kW3qdCQEGvJ9xlrr7I939ye7+kyT/MWv/G3suzpHpd9a/neTB9cp4LuZrne+xm/Z5IcxNPJbkyqq6fPp/vQ8kWZ4ps5zkpVllrs9koLX/izRn037d9yT5QHf/6Dpl/tJL4xWr6qpM/jt+5tzVcnFU1cunA3hTVS9Pcm2S980UW07ybTXxpZkMsP7wOa7qIln3/7B6Ns651Z8LNyX5r2uUeTjJtVX1yml3s2unx5ijqtqX5HuSfEN3f2KdMkPezzhLM2Omvylr/xsP+d7FfPzNJL/T3SfWOum5mK/TfI/dtM+LnWdX5fPDdParWzL5B9uR5N7uPlpVtyc50t3LmdyY+6tqJZMWuQNbV+Pz2pcl+dYk7101he6/SLI7Sbr7xzMJ0/+oqk4l+cMkBwTrTfO5Sd45zQc7k/yn7v7Fqnpz8qf343AmM1muJPlEkr+/RXU971XVZ2Uy+9u3rzq2+l54NjZJVf10kjcmeVVVnUjyr5LckeRnqurmJL+f5IZp2aUkb+7uf9jdJ6vqX2fy5TVJbu9uvTrOwjr34i1JLkzyrun71W9MZ5/+/CRv7+6vzTrvZ1vwEs4b69yLN1bV6zLpHnY80/er1fdive9dW/ASzhtr3YvuvidrjLH2XGyq9b7HbtrnhaUJAAAARkg3SwAAgBES5gAAAEZImAMAABghYQ4AAGCEhDkAAIAREuYAAABGSJgDAAAYIWEOAABghP4/u2Bt7+bKFA0AAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 1080x360 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.hist(q2_array/1e6,bins=70, weights=weights_np,density=True);\n",
"fig=plt.gcf()\n",
"fig.set_size_inches(15,5);"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"def lorentz_boost(fvector_to_boost, boost_vector):\n",
" \n",
" boost_vector_spatial=-boost_vector[:,0:3]\n",
" boost_vector_time_component=boost_vector[:,3].reshape(-1,1)\n",
" \n",
" beta = boost_vector_spatial/boost_vector_time_component\n",
" beta2 = scalar_product_3(beta,beta).reshape(-1,1)\n",
" \n",
" gamma = 1./np.sqrt(1.-beta2)\n",
" gamma_ = (gamma-1.)/beta2\n",
" \n",
" ve = fvector_to_boost[:,3].reshape(-1,1)\n",
" vp = fvector_to_boost[:,0:3]\n",
" \n",
" ve2 = gamma*(ve + scalar_product_3(vp,beta).reshape(-1,1))\n",
" vp2 = vp + (gamma*ve + gamma_*scalar_product_3(vp, beta).reshape(-1,1))*beta\n",
" \n",
" boosted_fvector=np.concatenate([vp2,ve2],axis=1)\n",
" return boosted_fvector"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"p1_q = lorentz_boost(l1, q_np)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"p2_q = lorentz_boost(l2, q_np)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(True, True, True)"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"(p1_q[:,0:3]+p2_q[:,0:3]).mean()<1e-10, np.sqrt(scalar_product_4(p1_q,p1_q)).mean()-105<1e-10, np.sqrt(scalar_product_4(p2_q,p2_q)).mean()-105<1e-10"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"cos_theta_l_np=get_costheta_l(p1_q,z)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"def dGamma_weights():\n",
" \n",
" Kst=particles_np['Kst']\n",
" l1 =particles_np['l1']\n",
" l2 =particles_np['l2']\n",
" q=particles_np['q']\n",
" \n",
" \n",
" q2 = scalar_product_4(q,q)\n",
" modq = np.sqrt(q2)\n",
" \n",
" l1_q = lorentz_boost(l1, q)\n",
" l2_q = lorentz_boost(l2, q)\n",
" \n",
"\n",
" cos_theta_l=get_costheta_l(l1_q,z)\n",
"\n",
" #G_weights = 1/(1-cos_theta_l*(mmu_mass)**2)\n",
" #print(cos_theta_l)\n",
" G_weights = (q2/(1-beta_l(q2,mmu_mass))) #*(1+cos_theta_l+cos_theta_l**2)\n",
" \n",
" return G_weights\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"dGamma=dGamma_weights()"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA3MAAAEvCAYAAADvmpjfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAWq0lEQVR4nO3df7Dl5X0X8PdHNqEmdUgKtDb86KJQZ5Iaa1xJnVrNlEJJqNlWwWzitKg4GC2OndoxREeM2M4stW1GLaODgoNpDUQ0ujPdlMbG0ZlOgrukSQgQzBa3ZQOTkMIQMaWU9OMf55C5OTl391z2/jjPPa/XzM495/t9vvc+5z77Pef7vs/zfZ7q7gAAADCWP7DTFQAAAGDjhDkAAIABCXMAAAADEuYAAAAGJMwBAAAMSJgDAAAY0J6drsCsc845p/fu3bvT1QAAANgR999//xe7+9xTlVu6MLd3794cPXp0p6sBAACwI6rqNxcpZ5glAADAgIQ5AACAAQlzAAAAAxLmAAAABiTMAQAADEiYAwAAGJAwBwAAMCBhDgAAYEDCHAAAwICEOQAAgAEJcwAAAAPas9MVAABgOe298ZcWLnv84FVbWBNgHj1zAAAAAxLmAAAABiTMAQAADEiYAwAAGJAJUAAAltRGJiDZiK2YrGSrJksxCQusT5gDAFgxWxUSge0lzAEAsK2ESdgcwhwAALuCIZmsGhOgAAAADEjPHADANjLEENgseuYAAAAGJMwBAAAMyDBLAABWjslS2A2EOQCAOVzsA8vOMEsAAIABCXMAAAADMswSAOA0WW4A2Al65gAAAAakZw4AAE7CZDgsKz1zAAAAA9IzBwAMb9GeE70mwG4izAEAK8NEJWw1f1hgOxlmCQAAMCBhDgAAYEDCHAAAwIDcMwcALCX3twGcnDAHAADbzNp1bAbDLAEAAAYkzAEAAAxImAMAABiQMAcAADAgYQ4AAGBAwhwAAMCAhDkAAIABWWcOANg2FgIH2DzCHABwWgQ0gJ2x0DDLqrqyqh6pqmNVdeOc/WdW1d3T/fdV1d7p9pdV1Z1V9UBVPVxV797c6gMAAKymU/bMVdUZSW5NcnmSE0mOVNWh7n5oTbHrkjzd3RdX1YEktyR5W5JrkpzZ3X+8ql6R5KGqen93H9/sFwIAALvRRnq/jx+8agtrwrJZpGfu0iTHuvvR7n4+yV1J9s+U2Z/kzunje5JcVlWVpJO8sqr2JPmDSZ5P8qVNqTkAAMAKWyTMnZfksTXPT0y3zS3T3S8keSbJ2ZkEu/+X5Ikkv5XkZ7r7qdkfUFXXV9XRqjr65JNPbvhFAAAArJpFwlzN2dYLlrk0yVeSvCbJRUn+XlX9ka8r2H1bd+/r7n3nnnvuAlUCAABYbYvMZnkiyQVrnp+f5PF1ypyYDqk8K8lTSd6R5Je7+/eSfKGqfi3JviSPnm7FAYCtY4ZKgOW3SM/ckSSXVNVFVfXyJAeSHJopcyjJtdPHVyf5SHd3JkMrv7cmXpnku5J8ZnOqDgAAsLpOGeam98DdkOTeJA8n+UB3P1hVN1fVW6fFbk9ydlUdS/LjSV5cvuDWJN+Y5NOZhMJ/192f2uTXAAAAsHIWWjS8uw8nOTyz7aY1j5/LZBmC2eOenbcdAACA07PQouEAAAAsF2EOAABgQAsNswQAAJbfRmaiPX7wqi2sCdtBzxwAAMCAhDkAAIABCXMAAAADcs8cAKyIjdxLA8DyE+YAYGACGsDqMswSAABgQMIcAADAgIQ5AACAAQlzAAAAAxLmAAAABiTMAQAADMjSBAAAsII2srTJ8YNXbWFNeKmEOQBYMtaOA2ARhlkCAAAMSJgDAAAYkDAHAAAwIGEOAABgQMIcAADAgMxmCQDbxCyVAGwmPXMAAAADEuYAAAAGJMwBAAAMSJgDAAAYkDAHAAAwILNZAsBpMEMlADtFzxwAAMCAhDkAAIABGWYJAACc1EaGlB8/eNUW1oS19MwBAAAMSJgDAAAYkDAHAAAwIGEOAABgQMIcAADAgMxmCQAzLAQOwAj0zAEAAAxImAMAABiQMAcAADAgYQ4AAGBAwhwAAMCAzGYJwEowQyUAu42eOQAAgAEJcwAAAANaKMxV1ZVV9UhVHauqG+fsP7Oq7p7uv6+q9q7Z9/qq+mhVPVhVD1TVN2xe9QEAAFbTKe+Zq6ozktya5PIkJ5IcqapD3f3QmmLXJXm6uy+uqgNJbknytqrak+QXkvxwd3+yqs5O8nub/ioAAIClsJF7lI8fvGoLa7L7LdIzd2mSY939aHc/n+SuJPtnyuxPcuf08T1JLquqSnJFkk919yeTpLt/u7u/sjlVBwAAWF2LhLnzkjy25vmJ6ba5Zbr7hSTPJDk7ybcn6aq6t6o+XlV///SrDAAAwCJLE9Scbb1gmT1J/mySP53ky0l+taru7+5f/ZqDq65Pcn2SXHjhhQtUCQAsNwDAalukZ+5EkgvWPD8/yePrlZneJ3dWkqem2/9Hd3+xu7+c5HCSN8z+gO6+rbv3dfe+c889d+OvAgAAYMUsEuaOJLmkqi6qqpcnOZDk0EyZQ0munT6+OslHuruT3Jvk9VX1imnI+/NJHgoAAACn5ZTDLLv7haq6IZNgdkaSO7r7waq6OcnR7j6U5PYk76uqY5n0yB2YHvt0Vf1cJoGwkxzubmNiAAAATtMi98yluw9nMkRy7bab1jx+Lsk16xz7C5ksTwAAp+Q+OABYzEKLhgMAALBchDkAAIABCXMAAAADEuYAAAAGJMwBAAAMSJgDAAAY0EJLEwDA6bDcAABsPj1zAAAAAxLmAAAABiTMAQAADEiYAwAAGJAwBwAAMCBhDgAAYECWJgDgJbHcAADsLD1zAAAAAxLmAAAABmSYJQAAsCM2MmT/+MGrtrAmY9IzBwAAMCBhDgAAYEDCHAAAwICEOQAAgAGZAAWAr2H9OAAYg545AACAAQlzAAAAAxLmAAAABiTMAQAADEiYAwAAGJAwBwAAMCBhDgAAYEDWmQNYAdaOA4DdR88cAADAgIQ5AACAAQlzAAAAAxLmAAAABiTMAQAADEiYAwAAGJClCQAGZbkBAFhteuYAAAAGJMwBAAAMSJgDAAAYkDAHAAAwIGEOAABgQMIcAADAgCxNALBELDcAACxKzxwAAMCA9MwBAABLb9HRK8cPXrXFNVkeC/XMVdWVVfVIVR2rqhvn7D+zqu6e7r+vqvbO7L+wqp6tqp/YnGoDAACstlOGuao6I8mtSd6c5LVJ3l5Vr50pdl2Sp7v74iTvTXLLzP73JvnQ6VcXAACAZLGeuUuTHOvuR7v7+SR3Jdk/U2Z/kjunj+9JcllVVZJU1Q8meTTJg5tTZQAAABa5Z+68JI+teX4iyRvXK9PdL1TVM0nOrqrfSfKuJJcnMcQSWElmqAQAtsIiPXM1Z1svWOafJHlvdz970h9QdX1VHa2qo08++eQCVQIAAFhti/TMnUhywZrn5yd5fJ0yJ6pqT5KzkjyVSQ/e1VX100leleT3q+q57v75tQd3921JbkuSffv2zQZFAAAAZiwS5o4kuaSqLkryuSQHkrxjpsyhJNcm+WiSq5N8pLs7yfe8WKCq3pPk2dkgBwAAwMadMsxN74G7Icm9Sc5Ickd3P1hVNyc52t2Hktye5H1VdSyTHrkDW1lpAACAVbfQouHdfTjJ4ZltN615/FySa07xPd7zEuoHAADAHAstGg4AAMByWahnDoCvZbkBAGCn6ZkDAAAYkDAHAAAwIGEOAABgQMIcAADAgEyAAjBlUhMAYCR65gAAAAYkzAEAAAxImAMAABiQMAcAADAgYQ4AAGBAwhwAAMCAhDkAAIABWWcO2PWsHwcA7EZ65gAAAAYkzAEAAAxImAMAABiQMAcAADAgYQ4AAGBAZrMEhmSGSgBg1emZAwAAGJAwBwAAMCBhDgAAYEDCHAAAwICEOQAAgAGZzRJYGmaoBABYnJ45AACAAQlzAAAAAxLmAAAABiTMAQAADEiYAwAAGJAwBwAAMCBLEwBbynIDAABbQ88cAADAgIQ5AACAAQlzAAAAAxLmAAAABiTMAQAADMhslgAAwK6xkZm0jx+8agtrsvWEOWDDLDcAALDzDLMEAAAYkDAHAAAwIGEOAABgQMIcAADAgBYKc1V1ZVU9UlXHqurGOfvPrKq7p/vvq6q90+2XV9X9VfXA9Ov3bm71AQAAVtMpZ7OsqjOS3Jrk8iQnkhypqkPd/dCaYtclebq7L66qA0luSfK2JF9M8he6+/Gq+o4k9yY5b7NfBHD6zFAJADCWRXrmLk1yrLsf7e7nk9yVZP9Mmf1J7pw+vifJZVVV3f3r3f34dPuDSb6hqs7cjIoDAACsskXC3HlJHlvz/ES+vnftq2W6+4UkzyQ5e6bMX0ry6939uy+tqgAAALxokUXDa8623kiZqnpdJkMvr5j7A6quT3J9klx44YULVAkAAGC1LdIzdyLJBWuen5/k8fXKVNWeJGcleWr6/PwkH0zyI939G/N+QHff1t37unvfueeeu7FXAAAAsIIW6Zk7kuSSqrooyeeSHEjyjpkyh5Jcm+SjSa5O8pHu7qp6VZJfSvLu7v61zas2sAiTmgAA7F6n7Jmb3gN3QyYzUT6c5APd/WBV3VxVb50Wuz3J2VV1LMmPJ3lx+YIbklyc5B9V1Sem/755018FAADAilmkZy7dfTjJ4ZltN615/FySa+Yc95NJfvI06wgAAMCMhRYNBwAAYLkIcwAAAANaaJglsDxMagIAQKJnDgAAYEjCHAAAwICEOQAAgAEJcwAAAAMS5gAAAAZkNktYEmapBABgI/TMAQAADEiYAwAAGJAwBwAAMCBhDgAAYEAmQIEtZFITAAC2ip45AACAAQlzAAAAAxLmAAAABiTMAQAADMgEKLBBJjUBAGAZ6JkDAAAYkDAHAAAwIGEOAABgQO6Zg7gPDgCA8eiZAwAAGJAwBwAAMCBhDgAAYEDCHAAAwIBMgMKuZVITAAB2Mz1zAAAAAxLmAAAABmSYJUMxdBIAACb0zAEAAAxImAMAABiQMAcAADAg98yx49wHBwAAG6dnDgAAYEDCHAAAwIAMs2RLGDoJAABbS88cAADAgIQ5AACAARlmyYYYPgkAAMtBzxwAAMCAhDkAAIABGWaJoZMAADAgPXMAAAAD0jO3S+ltAwCA3W2hnrmqurKqHqmqY1V145z9Z1bV3dP991XV3jX73j3d/khVff/mVR0AAGB1nbJnrqrOSHJrksuTnEhypKoOdfdDa4pdl+Tp7r64qg4kuSXJ26rqtUkOJHldktck+W9V9e3d/ZXNfiGrQG8bAADwokWGWV6a5Fh3P5okVXVXkv1J1oa5/UneM318T5Kfr6qabr+ru383yf+pqmPT7/fRzan++AQ0AADgpVgkzJ2X5LE1z08keeN6Zbr7hap6JsnZ0+0fmzn2vJdc20EIaAAAwFZbJMzVnG29YJlFjk1VXZ/k+unTZ6vqkQXqtZ3OSfLFna4ESbTFMtEWy0NbLA9tsTy0xfLQFstDW8yoW3bsR5+qLb5tkW+ySJg7keSCNc/PT/L4OmVOVNWeJGcleWrBY9PdtyW5bZEK74SqOtrd+3a6HmiLZaItloe2WB7aYnloi+WhLZaHtlgem9UWi8xmeSTJJVV1UVW9PJMJTQ7NlDmU5Nrp46uTfKS7e7r9wHS2y4uSXJLkf51upQEAAFbdKXvmpvfA3ZDk3iRnJLmjux+sqpuTHO3uQ0luT/K+6QQnT2US+DIt94FMJkt5IcmPmskSAADg9C20aHh3H05yeGbbTWseP5fkmnWO/akkP3UadVwGSzsEdAVpi+WhLZaHtlge2mJ5aIvloS2Wh7ZYHpvSFjUZDQkAAMBIFrlnDgAAgCUjzE1V1ZVV9UhVHauqG+fsP7Oq7p7uv6+q9m5/LXe/qrqgqv57VT1cVQ9W1d+dU+ZNVfVMVX1i+u+med+LzVFVx6vqgenv+uic/VVV/2J6bnyqqt6wE/Xc7arqj635P/+JqvpSVf3YTBnnxhapqjuq6gtV9ek1276pqj5cVZ+dfn31OsdeOy3z2aq6dl4ZFrdOW/yzqvrM9D3og1X1qnWOPen7GRuzTlu8p6o+t+Z96C3rHHvS6y42Zp22uHtNOxyvqk+sc6zzYpOsdx27lZ8XhlkmqaozkvzvJJdnspzCkSRv7+6H1pT520le393vrKoDSX6ou9+2IxXexarqW5N8a3d/vKr+UJL7k/zgTFu8KclPdPcP7FA1V0pVHU+yr7vnroUy/aD+O0nekuSNSf55d79x+2q4eqbvWZ9L8sbu/s01298U58aWqKo/l+TZJP++u79juu2nkzzV3QenF6Ov7u53zRz3TUmOJtmXyTqr9yf5U9399La+gF1knba4IpOZtF+omqwaNdsW03LHc5L3MzZmnbZ4T5Jnu/tnTnLcKa+72Jh5bTGz/2eTPNPdN8/ZdzzOi02x3nVskr+aLfq80DM3cWmSY939aHc/n+SuJPtnyuxPcuf08T1JLquqeYuicxq6+4nu/vj08f9N8nCS83a2VpzC/kw+PLq7P5bkVdM3M7bOZUl+Y22QY2t19//MZLbmtdZ+LtyZyQf2rO9P8uHufmr6gfzhJFduWUVXwLy26O5f6e4Xpk8/lsm6tmyxdc6LRSxy3cUGnKwtpterfznJ+7e1UivoJNexW/Z5IcxNnJfksTXPT+TrA8RXy0w/MJ5Jcva21G5F1WQo659Mct+c3X+mqj5ZVR+qqtdta8VWTyf5laq6v6qun7N/kfOHzXUg638oOze2z7d09xPJ5AM8yTfPKeP82H5/PcmH1tl3qvczNscN0yGvd6wznMx5sb2+J8nnu/uz6+x3XmyBmevYLfu8EOYm5vWwzY4/XaQMm6SqvjHJf0ryY939pZndH0/ybd39J5L8yyT/Zbvrt2K+u7vfkOTNSX50OpRjLefGNqqqlyd5a5L/OGe3c2P5OD+2UVX9w0zWtf3FdYqc6v2M0/evkvzRJN+Z5IkkPzunjPNie709J++Vc15sslNcx6572JxtpzwvhLmJE0kuWPP8/CSPr1emqvYkOSsvbWgBp1BVL8vkBPjF7v7Ps/u7+0vd/ez08eEkL6uqc7a5miujux+ffv1Ckg9mMjxmrUXOHzbPm5N8vLs/P7vDubHtPv/ikOLp1y/MKeP82CbTyQJ+IMlf6XUmBFjg/YzT1N2f7+6vdPfvJ/k3mf87dl5sk+k1619Mcvd6ZZwXm2ud69gt+7wQ5iaOJLmkqi6a/tX7QJJDM2UOJXlxVpmrM7nR2l+RNtl0XPftSR7u7p9bp8wffvF+xaq6NJP/x7+9fbVcHVX1yukNvKmqVya5IsmnZ4odSvIjNfFdmdxg/cQ2V3WVrPsXVufGtlv7uXBtkv86p8y9Sa6oqldPh5tdMd3GJqqqK5O8K8lbu/vL65RZ5P2M0zRzz/QPZf7veJHrLjbH9yX5THefmLfTebG5TnIdu2WfF3tOr8q7w3T2qxsy+YWdkeSO7n6wqm5OcrS7D2XSMO+rqmOZ9Mgd2Lka72rfneSHkzywZgrdf5DkwiTp7n+dSZj+W1X1QpLfSXJAsN4y35Lkg9N8sCfJf+juX66qdyZfbY/DmcxkeSzJl5P8tR2q665XVa/IZPa3v7lm29q2cG5skap6f5I3JTmnqk4k+cdJDib5QFVdl+S3klwzLbsvyTu7+29091NV9U8zuXhNkpu726iO07BOW7w7yZlJPjx9v/rYdPbp1yT5t939lqzzfrYDL2HXWKct3lRV35nJ8LDjmb5frW2L9a67duAl7Brz2qK7b8+ce6ydF1tqvevYLfu8sDQBAADAgAyzBAAAGJAwBwAAMCBhDgAAYEDCHAAAwICEOQAAgAEJcwAAAAMS5gAAAAYkzAEAAAzo/wN/Di51jrFbBAAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 1080x360 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.hist(q2_array/1e6,bins=70, weights=weights_np*dGamma,density=True);\n",
"fig=plt.gcf()\n",
"fig.set_size_inches(15,5);"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA2MAAAE6CAYAAACbL4ofAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd3jV5f3/8dfnnOydkL0JEMKSFYYgQ0Xce9fRai122Wprx69qtbW1ta22OFpX1Vrr3jgQRZEhG9kEErJDJtk755z79wfo11pHkCSf5Jzn47q8TgIn8YUX5pzX574/79syxggAAAAAMLAcdgcAAAAAAF9EGQMAAAAAG1DGAAAAAMAGlDEAAAAAsAFlDAAAAABsQBkDAAAAABv0WxmzLOtRy7JqLMva2YvnzrUsa4tlWS7Lsi74zO+lW5a1zLKsPZZl7bYsK7O/MgMAAADAQOnPlbHHJZ3Sy+eWSvqWpKc+5/eekPRnY8wYSdMl1fRFOAAAAACwU7+VMWPMSkn1n/41y7JGWJa11LKszZZlrbIsK+fwc4uNMdsleT7z/LGS/Iwx7xx+Xqsxpr2/MgMAAADAQBnoe8YeknSdMWaqpBsl/f0rnp8tqdGyrJcsy/rIsqw/W5bl7PeUAAAAANDP/AbqX2RZVpikWZKetyzr418O/Iov85M0R9JkHdrK+KwObWf8Z/+kBAAAAICBMWBlTIdW4RqNMZOO4GvKJX1kjCmUJMuyXpE0U5QxAAAAAEPcgG1TNMY0SyqyLOtCSbIOmfgVX7ZRUrRlWXGHPz9B0u5+jAkAAAAAA8IyxvTPN7aspyXNlxQrqVrSrZLek/QPSUmS/CU9Y4z5rWVZ0yS9LClaUqekKmPMuMPf5yRJd0myJG2WtMgY090voQEAAABggPRbGQMAAAAAfLGBnqYIAAAAABBlDAAAAABs0S/TFGNjY01mZmZ/fGsAAAAAGPQ2b95cZ4yJ+7Ln9EsZy8zM1KZNm/rjWwMAAADAoGdZVslXPYdtigAAAABgA8oYAAAAANiAMgYAAAAANqCMAQAAAIANKGMAAAAAYAPKGAAAAADYgDIGAAAAADagjAEAAACADShjAAAAAGADv948ybKsYkktktySXMaY3P4MBQAAgK/PGKMul0dtXS61d7vV1u1SW5db7d2HPne5jdzGyOMxcnsOfez2GDktS/5+lvydDvk7HQo4/Bjk71BEsL/Cg/wUEeSvkACnLMuy+48JDHm9KmOHHW+Mqeu3JAAAAPhSHd1uHWjqUGVjpw40dai6qVMH27pV/9l/2rvV7fL0Ww6nw1J4kJ+igv0VFx6o+PCgQ48Rhz5OigxSekyIkiKD5OdkIxbwRY6kjAEAAKAfGWNU3dylwrpWFdW1qai2TcUH23WgsUOVTR1qaO/5n68JD/RTTFiAYkIDlBQZpHHJEYoJC1BEkL/CAv0UEuA89Bjop9AAp4L8nfJ3OuR0SA7LktNhffLoMUY9bqMet0fdLo963B71uI06etxq6exRc4fr0GNnj1o6XWpo71FNc6f2VDbrg31dau1y/Vc2P4el1OhgpcWEKGNYiLJiw5STGK7RieEaFhY4UP9ZgUGrt2XMSFpmWZaR9KAx5qF+zAQAAODVPB6jsoZ25VW1KK+yRfuqW1RY16biujZ19Lg/eV6gn0MZw0KUGh2iyelRSo4KVlJkkJKjgpUcGayEyEAF+jlt/JP8t/Zul2pbulTR2KGy+naVHGxXSX27Sg+2a1vZATV3/l9Ziw0L1OjEMI1OiNAxqZGanB6l9JgQtj/Cp1jGmK9+kmUlG2MOWJYVL+kdSdcZY1Z+5jmLJC2SpPT09KklJSX9kRcAAGBI6exxa9eBJu2saFZeVbPyqlq0t6pF7d2HSpdlSekxIcqKDdXw2DANjwtVVmyoMmNDlRQRJIfDO8qJMUa1rV3ae/jPv7eqRXurDxXRzp5DWypjQgM0MTVSk9KiNSUjSrkZMQoOGDxlEzgSlmVt/qpZG70qY5/5prdJajXG/OWLnpObm2s2bdp0RN8XAABgqHN7jPJrWrS9rElbyxu1raxRe6ta5PIcer8VFeKvnMRw5SRGHHpMilB2QphCAnz3zhGX26N91a3aWtaorWUN+qi0UQW1rTJG8ndampwWrWNHDNOsEcM0KT1qUK0EAl+mT8qYZVmhkhzGmJbDH78j6bfGmKVf9DWUMQAA4AvaulzaUtqgDUX12lBUrx0VTZ+seIUH+WliapQmpkVqYmqUjkmNUkJEINvweqG5s0dbShq0tvCg1u4/qB0VTTJGCvJ3aNaIWJ04Jl4n5iQoMTLI7qjAF+qrMpYl6eXDn/pJesoY8/sv+xrKGAAA8EYNbd3aWHyoeG0srtfOA81ye4wcljQuOVJTM6I/KV+Zw0K9Zouh3Zrae7S+6KDWFNTpvb01KqvvkCSNS47QiWMStHBsgsYlR1B0Maj0yzbF3qCMAQAAb9DR7daG4nqtzq/Vqvw65VW1SJIC/ByalBal6Zkxmj48RlMyohUW6LtbDQeSMUb5Na16d0+13ttToy2lDfIYKSs2VGdMTNZZE5M0Mj7c7pgAZQwAAOBIeDxGuyubtSq/TqsLarWxuEHdLo8CnA7lZkZr9shYTR8eo2NSI7l3aZA42NqlZburtWTbAa0tPChjpDFJETpzYpLOn5KqhAi2MsIelDEAAICv0NTRo5X7arV8T7VW5tepvq1bkpSTGK7jRsZqTnacpmcy1W8oqGnu1Bs7KrVk2wFtKW2U02Hp+NFxuig3TcfnxMufA6gxgChjAAAAn6Oork3L91Rr+Z4abSyul8tjFBMaoPnZcZqTHavZI2MVH86KylBWVNem5zaV6YXN5apt6VJceKAumJqqy2akKzU6xO548AGUMQAAAB0aOb+puF7vHi5ghXVtkqTRCeGHJvONidektGg5GbjhdVxuj97fW6tnN5bp/b01MsbolPGJunr2cE3NiGboB/oNZQwAAPgsl9ujDUX1enNnpZburFZda5cCnA7NyIrRgjEJOiEnXmkxrJD4korGDj2xtlhPry9Vc6dLx6RG6urZw3X6MUlsYUSfo4wBAACf0uP2aO3+g3prZ6Xe3lWt+rZuBfs7dXxOnE4dn6Tjc+KZegi1d7v00pYKPbqmSIW1bUqNDtb35o/QBVNTGcyCPkMZAwAAXs/tMfpwf51e23pAy3ZXq6mjR6EBTp04JkGnTUjUvOx4hm/gc3k8Ru/vrdG97xVoa1mjEiOCdO28LF0yLZ2/MzhqlDEAAOCVjDHaUdGkVz46oCXbD6i2pUvhgX46aWyCTp2QpDmjYhXkz5tp9I4xRmsKDure9/K1vqhesWEB+sHxI/WNGemslOFro4wBAACvUlzXple3HtCrWytUWNemAKdDx+fE6ZxJKTo+J54ChqO2oahef31nn9YWHlRqdLB+ujBbZ09MkYPhLjhClDEAADDkNbR167VtB/TyRxXaWtYoy5JmDI/ROZNSdOr4JEWG+NsdEV7GGKNV+XW6c2medh1oVk5iuH5xao7mZ8cxfRG9RhkDAABDkttjtCq/Vs9vLtc7u6rV7fZoTFKEzpmUrLMmJSspMtjuiPABHo/R6zsq9Ze396q0vl3zsuP06zPHakRcmN3RMARQxgAAwJBScrBNz28q14tbylXZ1KnoEH+dMzlFF05N09jkCLvjwUd1uzx6Ym2xFr+br06XW1fPHq7rThzFZE58KcoYAAAY9Nq7XXprR5We21Sm9UX1cljS3Ow4XZSbphPHxDNAAYNGbUuX/rQ0T89vLldceKB+dVqOzpmUwtZFfC7KGAAAGLT2VrXoyXUlevmjCrV2uZQxLEQX5abpvCkpbEPEoLa1rFG3vrZL28oaNWdUrO44dwIHiON/UMYAAMCg0uVy660dVXpyXYk2lTQowM+h0yck6ZJpaZo+PIYVBgwZHo/Rf9aX6I9v5cljpJ8uzNZVs4fLydRFHEYZAwAAg0LJwTY9tb5Uz28uV31btzKGheiyGem6YGqaYkID7I4HfG0HGjt08ys79V5ejSamRupPF0zU6MRwu2NhEKCMAQAA27jcHi3Pq9GT60q0Kr9OToelk8Yk6LKZ6Zo9IpZzm+A1jDF6bdsB/WbJbrV2uvSzk0fr28cN5++4j+tNGWMEDAAA6FON7d16ZmOZ/r22RBWNHUqMCNINC7J18bQ0JUYG2R0P6HOWZensSSk6bmSs/t9LO/T7N/fovbwa3XXRRCVHcf8jvhgrYwAAoE/kV7fosQ+L9dKWcnX2eDQzK0bfmjVcC8bEy8/psDseMCCMMXpuU5l+s2S3nA5LvztnvM6elGJ3LNiAlTEAANCvPB6jFftq9NiaYq3Kr1OAn0PnTErWVbOHa0wS54LB91iWpYunpWtm1jDd8OxW/fiZrVqVX6fbzx6v4ACOacB/o4wBAIAj1trl0gubyvSvtSUqqmtTQkSgblyYrUunp2tYWKDd8QDbZQwL1XPXHqt7lufr3vcLtKO8SfdfNkUj48PsjoZBhG2KAACg1yqbOvTYmmI9vb5ULV0uTUqL0lWzM3XahCT5sxUR+Fwf7KvVDc9uVVePW384/xidNTHZ7kgYAGxTBAAAfWJvVYseWlmo17ZVyO0xOm1Ckr593HBNTo+2Oxow6M3LjtMbPzpO1z31kX709EfaWFSvW84YqwA/LmD4OsoYAAD4XMYYrSus14Mr92vF3loF+zt12YwMffu44UqLCbE7HjCkJEUG6+lFM/WnpXl6eFWR9lW36O+XTWFbr4+jjAEAgP/i9hgt3VmlB1fu1/byJg0LDdBPT8rW5TMzFM0BzcDX5u906KbTx2pccqR+/uJ2nX3/Gj18ZS7DbnwYZQwAAEiSOrrdemFzmR5eVaTS+nYNjw3V788dr/OnpCrInylwQF85Z3KKhseGatG/N+n8f3youy+apFPGJ9odCzZggAcAAD6upbNHT6wt0aOri3SwrVuT0qL03XlZOmlsopwOy+54gNeqbu7Uon9v1rayRv3ilBx9d16WLIv/57wFAzwAAMAXamjr1mNrivT4h8Vq7nRpXnacvj9/hKYPj+ENITAAEiKC9Oyimbrx+W26c2meqpo69Oszx3ERxIdQxgAA8DE1LZ16ZFWRnlxXovZut04el6AfHj9KE1Ij7Y4G+Jwgf6fuuWSykiKD9PCqIlU1d2rxJZPZGuwjKGMAAPiIisYOPfjBfj2zsUwut0dnTkzW9+eP1OjEcLujAT7N4bB00+ljlRQZrNvf2K3LHlmvR67MZWCOD6CMAQDg5Yrq2vSPFQV6aUuFLEs6b3Kqvjd/hDJjQ+2OBuBTrj5uuBIjg3T9s1t14YNr9Z9rZighIsjuWOhHlDEAALzU/tpW3bM8X0u2HZC/06HLZ2Zo0dwsJUcF2x0NwBc4bUKSYkID9O3HN+rCBw4VMs71815MUwQAwMsU1bXpnuX5enVrhQL9nLri2Ax9Z06W4sI5XBYYKj4qbdA3H92gsEA/PXnNDGXFhdkdCUeoN9MUKWMAAHiJkoNtumd5gV7+qFwBfg5dMTND184bodgwShgwFO0+0Kwr/rlelmXpP9fM4P7OIYYyBgCADyg92K5738vXSx9VyM9h6fKZGbp2Xpbiw7nXBBjqCmpadNkj69Xt8uip78zUmKQIuyOhlyhjAAB4sbL6dt3/foFe2Fwuh8PSN6an6/vzRyieG/4Br1JysE2XPLRO3S6Pnlk0U6MSWCEbCihjAAB4oYrGDt3/foGe31QmS5YunZ6m780fqcRIShjgrYrq2nTRg2slSc8umsk9ZEMAZQwAAC9S29Kl+98v0FPrS2VkdPG0NH1//kimIwI+Ir+6RZc8tE7+ToeevXamMoZxPMVgRhkDAMALNHX06OGVhXp0TZG6XB5dODVV1504SimUMMDn7Kls1qUPr1NogJ+e++6x/BwYxChjAAAMYR3dbv1rbbH+sWK/mjp6dPoxSfrJSdkawfYkwKftrGjSpQ+vU3x4oF747ixFhwbYHQmfozdlzHEE38xpWdZHlmW9fvTRAADAF+lxe/TvdSWa9+f39ce38jQ5PUqvX3ec7v/GFIoYAI1PidTDV+aqrKFDVz2+Ue3dLrsj4WvqdRmT9GNJe/orCAAAvs7jMXrlowqdeNcHuuWVnUqPCdFz1x6rx6+arvEpkXbHAzCIzMwapnsumaTt5Y36wX+2qMftsTsSvoZelTHLslIlnS7pkf6NAwCA7zHG6N3d1TrtnlW6/tmtCg3006PfytXz3z1W04fH2B0PwCB1yvgk3X7OeL2/t1a/fHGH+uP2I/Qvv14+72+Sfi6JQw0AAOhDG4rq9ce39mhLaaMyh4Xonksn64wJSXI4LLujARgCLpuRodqWLv3t3XwlRQbpxpNH2x0JR+Ary5hlWWdIqjHGbLYsa/6XPG+RpEWSlJ6e3mcBAQDwRgU1rbpzaZ7e2V2thIhA3XHuBF2Ymyp/55HcQQAA0o9PHKXKxk7d936BsuJCdd6UVLsjoZd6szI2W9JZlmWdJilIUoRlWU8aYy7/9JOMMQ9Jekg6NE2xz5MCAOAFalo69bd38/XsxjIF+zv1s5NH6+rZwxUc4LQ7GoAhyrIs3X7OeJXUt+mXL+5QekyIcjPZ4jwUHNFo+8MrYzcaY874sucx2h4AgP/W1uXSQysL9fCqQnW7PLp8ZoauO2GkhoUF2h0NgJdobO/WuX//UE0dPXr1B7OVFhNidySf1qej7QEAwJHrcXv05LoSzfvzCi1enq/jR8fr3Z/M021njaOIAehTUSEB+uc3c+X2GF39+EY1d/bYHQlfgUOfAQDoB8YYLdtdrTuX5qmwtk3TM2P0/07L0eT0aLujAfByHxbU6cpHN2hedpwevjKXgUA2YWUMAAAbbC5p0IUPrNW1/94sS9LDV+bq2WtnUsQADIhZI2P16zPHanleje57v8DuOPgSvR1tDwAAvkJxXZvuXJqnt3ZWKS780ITEi3JT5ceERAAD7IqZGfqotFF/fXefJqRG6vjR8XZHwuegjAEAcJSaOnp033v5evzDYvk7HbphQbaumTNcoYG8zAKwh2VZuuPcCcqratH1z2zVkh8ep/RhDPQYbLhUBwDA1+Rye/TvtcU6/i8r9MjqIp07OUUrbpyvHy8YRREDYLvgAKceuHyKjDH67pOb1dHttjsSPoMyBgDA17Bib41OXbxKt7y6S9kJYVryw+P0pwsmKj4iyO5oAPCJjGGhWnzJZO2patatr+20Ow4+g8t2AAAcgfzqFv3ujT36YF+tMoeF6MErpmrh2ARZFtPKAAxOx+fE6/vzR+j+9/fruFFxOmtist2RcBhlDACAXqhv69Zf39mnpzaUKiTAqZtPH6Mrj81UgB+bTAAMftcvyNa6wnr96qUdmpQaxf1jgwSvIAAAfIkul1sPryzUvD+/r6c2lOqyGen64GfH65o5WRQxAEOGv9OhxZdMksOSrnt6i7pdHrsjQZQxAAA+lzFGS3dWaeFfV+r3b+7R1IxoLf3xHP327PGKCQ2wOx4AHLHU6BDdef4x2lbepLuW7bU7DsQ2RQAA/sfOiibd/vpurS+q16j4MP3r6umalx1ndywAOGqnTkjSZTPS9eDKQs0eGau5/GyzFWUMAIDDDrZ26S/L9umZjaWKDgnQ7eeM16XT0ji0GYBXueWMsVpfVK+fv7Bdb18/V5Eh/nZH8lm8ugAAfJ7L7dHja4p0/F9W6LlNZbpq1nC9f+N8XTEzgyIGwOsE+Tt190UTVdvapduW7LI7jk9jZQwA4NPWFNTpN0t2aV91q+aMitWvzxirUQnhdscCgH51TGqUfnj8SC1enq+TxyXolPFJdkfySZQxAIBPKqtv1+/f2KOlu6qUFhOsh66YqpM4LwyAD/nhCSO1PK9aN728U7mZMYoNC7Q7ks9h7wUAwKd0dLt197K9WnD3B/pgX61+dvJovXPDPC0cl0gRA+BT/J0O3X3RJLV0ufSrl3bIGGN3JJ/DyhgAwCcYY/T69kr94c09OtDUqbMnJeuXp+YoKTLY7mgAYJvshHDduDBbd7yZpyXbK3XWxGS7I/kUyhgAwOvtPtCs25bs0oaieo1LjtDiSydrWmaM3bEAYFD49nFZemN7pX67ZJfmjopVVAhnKQ4UtikCALxWQ1u3bn5lh864d5UKalp1x7kT9NoPj6OIAcCnOB2W/nDeMWpo79Edb+6xO45PYWUMAOB1XG6PntpQqruW7VNrl0tXHpupGxZkc5YOAHyBsckR+s6cLD3wwX6dMzlFs0bE2h3JJ1DGAABeZWNxvW55Zafyqlo0e+Qw3XrmOGUzqh4AvtL1C0bprZ2V+tVLO7T0+rkK8nfaHcnrsU0RAOAValu69JPnturCB9aqpdOlBy6foie/PYMiBgC9FOTv1B3nTlDxwXbd+16+3XF8AitjAIAhzeX26N/rSnT3sn3qcnn0g+NH6AfHj1RIAC9xAHCkZo+M1flTUvXQykKdNyVVI+LC7I7k1VgZAwAMWRuL63XGvav1myW7NSk9Skuvn6OfnZxDEQOAo/DLU3MU5O/Uba/t4uyxfkYZAwAMOZ/ektjc0aMHLp+iJ66eriyu4ALAUYsLD9QNC7K1Kr9Oy3ZX2x3Hq3HpEAAwZHx6S2Kny82WRADoJ1ccm6FnNpbq9td3a152HMM8+gkrYwCAIeGzWxLfvn4uWxIBoJ/4Ox267axxKm/o0AMf7Lc7jtfiFQwAMKjVtnTpD2/t0UtbKpQcGaQHLp+ik8clyrIsu6MBgFebNSJWpx+TpH+s2K/zp6QqLSbE7kheh5UxAMCg5HJ79NiaIp3wlxVasu2Avj9/hN796TydMj6JIgYAA+Sm08bIsqQ7l+bZHcUrsTIGABh0Pn1w85xRsbrtrHGMVwYAGyRHBes7c7J073sFumZOoyalRdkdyauwMgYAGDQ+OyXxH5cdmpJIEQMA+1w7b4RiwwJ0xxt7GHXfx1gZAwDYzu0xempDqf60NE+dPW59f/4I/fAEpiQCwGAQFuinHy/I1i2v7NQ7u6u1cFyi3ZG8Bq9yAABb7axo0k2v7NS2skbNHjlMvz17PCthADDIXDItTY+tKdIfl+bphJx4+TnZYNcX+K8IALBFa5dLt7++W2fdt1oVDe3628WT9OS3Z1DEAGAQ8nc69MtTclRY26ZnNpbZHcdrsDIGABhQxhi9vatav1myS5VNnfrGjHT94uQcRYb42x0NAPAlThqboOnDY/S3d/N13pQUtpL3AVbGAAADpryhXdf8a5O+++RmRQb768XvzdId506giAHAEGBZln5xymjVtXbpibUldsfxCtRZAEC/63F79M/VRVr8br6kQ+fWXDU7k3sOAGCImZoRo3nZcXrwg/26fGaGwgKpE0eDV0EAQL/aVFyvM+5ZrT++lafjRsXq3Z/O03fmZlHEAGCIuuGkbDW09+jxNUV2RxnyqLIAgH7R2N6tP76Vp2c2lik5MkgPXTGVccgA4AUmpUVpwZh4PbSyUFccm6nIYLaaf11clgQA9CljjF7cXK4T7vpAz28u16K5WXrnJ/MoYgDgRW44KVvNnS79czWrY0eDlTEAQJ8pqGnVza/s0LrCek1Oj9Lvz5mgsckRdscCAPSxccmROnV8oh5dXaSrZmUqOjTA7khD0leujFmWFWRZ1gbLsrZZlrXLsqzfDEQwAMDQ0dnj1l3L9urUxSu1+0Cz7jh3gl787iyKGAB4sesXZKu1y6XHPyy2O8qQ1ZuVsS5JJxhjWi3L8pe02rKst4wx6/o5GwBgCFi5r1a3vLpTJQfbde7kFP3qtDGKCw+0OxYAoJ+NTgzXSWMT9PiHxfrO3CwmK34NX7kyZg5pPfyp/+F/TL+mAgAMejXNnbru6Y905aMb5LAs/eeaGfrrxZMoYgDgQ74/f4SaOnr09PpSu6MMSb2qr5ZlOSVtljRS0v3GmPX9mgoAMGi5PUZPrS/Rn5buVZfboxsWZOvaeVkK8nfaHQ0AMMAmp0dr1ohhenhVoa6claFAP14LjkSvpikaY9zGmEmSUiVNtyxr/GefY1nWIsuyNlmWtam2travcwIABoGdFU067+9rdMuru3RMWqTevn6ufrxgFEUMAHzYD44fqZqWLr24ucLuKEPOEY22N8Y0Sloh6ZTP+b2HjDG5xpjcuLi4PooHABgMWrtc+u2S3TrrvtWqaOzQ4ksm6clvz9Dw2FC7owEAbDZrxDBNTI3UAx/sl8vtsTvOkNKbaYpxlmVFHf44WNICSXn9HQwAYD9jjJburNSCuz7QYx8W6dLp6Vr+k/k6e1KKLMuyOx4AYBCwLEvfP36kSuvb9caOSrvjDCm9uWcsSdK/Dt835pD0nDHm9f6NBQCwW1l9u259bZfey6vRmKQI/f3yKZqSHm13LADAIHTSmARlxYXqn6uLdNbEZC7Y9dJXljFjzHZJkwcgCwBgEOhxe/TIqiItXr5PDsvSzaeP0bdmZcrPeUQ72wEAPsThsHTV7OG65ZWd2lzSoNzMGLsjDQm8sgIAPrGxuF6n37NKdy7N09xRcXrnJ/N0zZwsihgA4CudPyVFkcH+enRNkd1RhgxOZgMAqKGtW398K0/PbipTcmSQHr4yVyeNTbA7FgBgCAkJ8NOl09P10Mr9Km9oV2p0iN2RBj0udQKADzPG6IXN5Trx7g/0wpZyXTs3S+/8ZB5FDADwtVx5bIYsy9K/Piy2O8qQwMoYAPiogpoW3fTyTq0vqteU9Cj9/twJGpMUYXcsAMAQlhwVrNMmJOmZjWX68YJshQVSN74MK2MA4GM6e9z6y9t7deriVdpT2aw7zp2gF747iyIGAOgTV8/OVEunSy9sKrM7yqBHVQUAH/LBvlrd8spOlda367zJKfrV6WMUGxZodywAgBeZnB6tiWlRenJ9qb45K5Mx91+ClTEA8AE1zZ364VNb9M1HN8jPYempa2bo7osnUcQAAP3ishnpKqhp1cbiBrujDGqUMQDwYm6P0RNri3XiXR9o2e5q3bAgW29dP0ezRsbaHQ0A4MXOPCZZ4UF++s/6ErujDGpsUwQAL7Wzokk3vbxD28qbdNzIWN1+zngNjw21OxYAwAcEBzh1/pRUPbW+VL8+o0vD2N4qZUUAACAASURBVInxuVgZAwAv09LZo98s2aWz7lutisZOLb5kkv797ekUMQDAgLpsRrq63R69sLnc7iiDFitjAOAljDFaurNKty3ZpZqWLl02I10/OzlHkcH+dkcDAPigUQnhmj48Rk9tKNV35mTJ4WCQx2exMgYAXqCsvl1XP75R3/vPFg0LDdRL35ul350zgSIGALDVZTPSVXKwXasL6uyOMiixMgYAQ1i3y6NHVhfqnuX5cliWbj59jL41K1N+Tq61AQDsd8r4REWH+Ou5TWWamx1nd5xBhzIGAEPUxuJ63fTyDu2rbtXJ4xJ065njlBwVbHcsAAA+Eejn1FkTk/X0xjI1tfcoMoQdG5/GpVMAGGIa2rr1ixe268IH1qqty61HrszVg1fkUsQAAIPSBVPT1O3y6PUdB+yOMuiwMgYAQ4QxRi9sLtcdb+5Rc6dL187L0o9PHKWQAH6UAwAGr/EpERqdEK4XNpfrshkZdscZVHgFB4AhoKCmRTe9vFPri+o1NSNavz93vHISI+yOBQDAV7IsS+dPTdEdb+Zpf22rRsSF2R1p0GCbIgAMYp09bv3l7b06dfEq5VW16A/nTdDz1x5LEQMADCnnTEqR02HpRc4c+y+sjAHAILVib41+/eoulda367zJKfrV6WMUGxZodywAAI5YfESQ5mXH6aUtFfrpwtFycuaYJFbGAGDQqW7u1A+e2qJvPbZRfk5LT31nhu6+eBJFDAAwpJ0/JVVVzZ36cD9njn2MlTEAGCTcHqMn15XoL2/vVZfbo5+clK1r52Up0M9pdzQAAI7aiWPiFRbopyXbDmjOKM4ckyhjADAo7Chv0k2v7ND28ibNGRWr288er8zYULtjAQDQZ4L8nVo4NkFLd1bpd+dMUIAfm/T4LwAANmru7NFtr+3S2fev1oHGTi2+ZJKeuHo6RQwA4JXOnJis5k6XVuXX2h1lUGBlDABsYIzRGzsq9dslu1Xb2qUrZmbopwtHKzLY3+5oAAD0m9kjYxUZ7K8l2w7oxDEJdsexHWUMAAZYycE23fLqLq3cV6txyRF6+MpcTUyLsjsWAAD9LsDPoVPHJ2rJtgPq7HEryN+374tmmyIADJAul1v3LM/XSX9dqS0lDbr1zLF69QezKWIAAJ9y5sRktXW79X5ejd1RbMfKGAAMgA8L6nTzqztVWNum0yck6ZYzxioxMsjuWAAADLiZWcMUGxaoJdsP6NQJSXbHsRVlDAD6UW1Ll+54c49e/qhC6TEhevyqaZo/Ot7uWAAA2MbpsHTahEQ9u7FM7d0uhQT4biVhmyIA9AOPx+g/60t04l0r9Pr2A7ruhJFadsNcihgAAJJOGZeoLpdHq/J9+wBo362hANBPdh9o1k2v7NBHpY2amRWj350zQSPjw+yOBQDAoDFteIwig/21bFe1Th6XaHcc21DGAKCPtHa59Ld39umxD4sVFeyvuy+aqHMnp8iyLLujAQAwqPg7HToxJ17L86rlcnvk5/TNDXuUMQA4SsYYvb2rWr9ZskuVTZ26dHq6fnHKaEWFBNgdDQCAQWvhuAS99FGFNhY36NgRw+yOYwvKGAAchbL6dt322i4tz6tRTmK47vvGFE3NiLY7FgAAg97c7DgF+jm0bHcVZQwA0HvdLo8eWV2oe5bny2FZuvn0MfrWrEyf3WYBAMCRCgnw05xRsVq2q1q/PmOsT27rp4wBwBHaUFSvm17eofyaVp08LkG3njlOyVHBdscCAGDIWTg2Ue/uqdHuymaNS460O86Ao4wBQC/Vt3XrD2/u0fOby5USFaxHrszVgrEJdscCAGDIOnFMvCxLend3DWUMAPC/PB6j5zaV6c6leWrpdOm780boRyeO9OlDKgEA6AvDwgJ1TGqUVuyr0Y8XjLI7zoDjnQQAfImdFU265dWd+qi0UdMyo/W7cyZodGK43bEAAPAa87PjdM97+Wpo61Z0qG9NIuZOcwD4HE0dPbr11Z06677VKqtv110XTtRz1x5LEQMAoI/NGx0nY6RVBXV2RxlwrIwBwKcYY/TSlgr94a09qm/r1uUzM/TThaMVGexvdzQAALzSxNQoRYf4a8XeGp01MdnuOAPqK8uYZVlpkp6QlCjJI+khY8zi/g4GAANtb1WLbnl1pzYU1WtSWpQev2q6xqf43s3EAAAMJKfD0pxRcVq5r1Yej5HD4Tsj7nuzMuaS9FNjzBbLssIlbbYs6x1jzO5+zgYAA6K1y6XF7+7To2uKFR7kpz+eN0EX5ab51IsBAAB2mj86Tq9tO6Ddlc0+dSH0K8uYMaZSUuXhj1ssy9ojKUUSZQzAkGaM0Rs7KnX767tV3dylS6en6ecn5/jczcMAANhtbnacJGnF3hrK2BexLCtT0mRJ6z/n9xZJWiRJ6enpfRANAPrP/tpW3frqLq0uqNO45Ag9cPlUTU6PtjsWAAA+KTYsUBNSIrVyX51+eILvjLjvdRmzLCtM0ouSrjfGNH/2940xD0l6SJJyc3NNnyUEgD7U0e3Wve/l6+FVhQryd+q3Z4/TZTMy5GRLIgAAtpo1cpgeXV2kjm63ggOcdscZEL0abW9Zlr8OFbH/GGNe6t9IAND3jDF6e1eVFtz9gf6+Yr/OnJis9346X1cem0kRAwBgEJg1IlY9bqNNJfV2RxkwvZmmaEn6p6Q9xpi7+z8SAPSt0oPtum3JLr2XV6PRCeF67tpjNX14jN2xAADAp+RmRMvPYenD/Qc1Z1Sc3XEGRG+2Kc6WdIWkHZZlbT38a78yxrzZf7EA4Oh19rj14AeF+vuKAvk5LN18+hh9c1am/J2cdw8AwGATGuinSWlR+nD/QbujDJjeTFNcLYk9PACGDGOM3t1To9++vktl9R0645gk3Xz6WCVGBtkdDQAAfIlZI4bpvvcL1NzZo4ggf7vj9DsuDwPwKoW1rfrWYxv1nSc2KcjPqaeumaH7vjGFIgYAwBAwc8QweYy0scg37hs7otH2ADBYtXW5dO97Bfrn6kIF+Tl1yxljdeWxGWxJBABgCJmSHq0AP4c+3H9QJ45JsDtOv6OMARjSjDF6bdsB/eHNPFU1d+qCqan6xSk5igsPtDsaAAA4QkH+TuVmRGutj9w3RhkDMGTtqWzWra/t0oaiek1IidT9l03R1AwObgYAYCiblhmje9/LV0tnj8K9/L4xyhiAIaepvUd3v7NX/15Xoshgf91x7gRdPC2N88IAAPACuZnR8hjpo9JGzc327hH3lDEAQ4bHY/TcpjL96e29amzv1uUzM/STk7IVFRJgdzQAANBHJqdHy2FJm0oaKGMAMBh8VNqgW1/bpe3lTZqWGa3fnDVDY5Mj7I4FAAD6WFign3ISI7SlpMHuKP2OMgZgUKtt6dKflubp+c3lig8P1OJLJumsicmyLLYkAgDgrXIzo/Xi5nK53B75efFkZMoYgEGpx+3RE2tL9Ld39qnT5da1c7N03YmjFBbIjy0AALzd1IxoPbG2RHlVLRqfEml3nH7DuxoAg877e2v0u9d3a39tm+aMitVtZ43TiLgwu2MBAIAB8vF05M0lDZQxABgIBTWt+t0bu7Vib60yh4XokStzdeKYeLYkAgDgY1KigpUYEaTNJQ365qxMu+P0G8oYANs1dfRo8bv5emJtsYL9nbrptDH65qxMBfh57x5xAADwxSzL0tSMaG0p9e4hHpQxALZxe4ye3lCqu9/Zp4b2bl0yLU0/XThasWGBdkcDAAA2OyY1Um/sqFRDW7eiQ73zGBvKGABbfFhQp9++vlt5VS2aPjxGvz5jrFfvCQcAAEdmQuqh9wU7Kpq89rwxyhiAAVV6sF2/f3O33t5VrZSoYP39sik6dXwi94UBAID/8vFF2u3ljZQxADgarV0u3f9+gf65qkh+Tks3LszWNXOyFOTvtDsaAAAYhCKC/JUVG6rt5U12R+k3lDEA/crjMXpxS7n+9PZe1bZ06bzJKfr5KTlKjAyyOxoAABjkJqRGakNRvd0x+g1lDEC/WV94UL9/c4+2lzdpcnqUHrpiqianR9sdCwAADBETUiL16tYDqmnpVHy4913IpYwB6HNFdW36w5t7tGx3tZIig/TXiyfq7Ikpcji4LwwAAPTeMalRkqSdFU06IYcyBgBfqKGtW4uX5+vJdSUK9HPoxoXZ+vZxWQoO4L4wAABw5MYlR8iypG1lTTohJ8HuOH2OMgbgqHW53HriwxLd+16+Wrtcunhaum44aZRXbicAAAADJzTQT8NjQ7W7stnuKP2CMgbgazPG6M0dVbpzaZ5K69s1LztOvzptjEYnhtsdDQAAeIkxiRHaUeGdExUpYwC+li2lDfr9G3u0uaRBOYnheuLq6V57BggAALDP6MRwvbGjUm1dLoUGeld98a4/DYB+V1bfrjuX5un17ZWKCw/UnedP0AVT0+RkOAcAAOgHOYd33OytbtEUL5vKTBkD0CtNHT36+/sFemxNsRwO6UcnjtK1c7O87goVAAAYXHISIyRJe6soYwB8TJfLrSfXleq+9/LV2NGj86ek6saFozm0GQAADIjU6GCFBji1t6rF7ih9jjIG4HN5PEavbqvQX97ep4rGDs0ZFatfnJKj8SmRdkcDAAA+xOGwlJ0Yrj1eOFGRMgbgvxhj9MG+Wt25dK/2VDZrfEqE7jz/GB03KtbuaAAAwEflJEborZ2VMsbIsrznPnXKGIBPbCtr1B/fytPawoNKjwnRPZdO1hkTkuRgOAcAALBRTmK4nt5QqurmLq+6VYIyBkDFdW3687K9emN7pWJCA3TbmWP1jRkZCvBz2B0NAABAo+LDJEkFNa2UMQDeobalS/csz9fTG0rl73ToRyeM1HfmZik8yN/uaAAAAJ/IijtUxgrrWr3q1gnKGOCDWrtcenhloR5eVagul0eXTk/Tj04cpfhw77nSBAAAvEdCRKBCA5wqrG2zO0qfoowBPqSzx60n15XoHyv262Bbt06bkKgbF47+5GoTAADAYGRZlobHhWp/bavdUfoUZQzwAT1uj57fVK57luerqrlTs0cO040LR2uylx2cCAAAvFdWbJg2lzTYHaNPUcYAL+b2GC3ZdkB/fXefSg62a3J6lO6+aKJmjfSevdYAAMA3ZMWFasn2A+rscSvI32l3nD5BGQO8kDFGb++q1t3v7NW+6laNSYrQP7+ZqxNy4r3qbA4AAOA7suLCZIxUVNemMUkRdsfpE5QxwIsYY7Qyv053Ldur7eVNyooL1X3fmKzTxnNWGAAAGNqyYkMlSYW1lDEAg8zG4nr9+e292lBUr5SoYP3pgmN03uQU+Tk5KwwAAAx9WXEflzHvGeJBGQOGuO3ljbr7nX1asbdWceGB+u3Z43TxtDQF+nnHXmoAAABJCgnwU0JEoErq2+2O0mcoY8AQtb28UYvfzdfyvBpFhfjrl6fm6JvHZio4gBIGAAC8U2p0iMp8qYxZlvWopDMk1Rhjxvd/JABf5rMl7Gcnj9aVx2YoPMjf7mgAAAD9Ki06WBuLvWe8fW9Wxh6XdJ+kJ/o3CoAvQwkDAAC+Li0mRK9tO6Aet0f+XnBf/FeWMWPMSsuyMvs/CoDP8+kSFhlMCQMAAL4rNTpYHiNVNXUqLSbE7jhHrc/uGbMsa5GkRZKUnp7eV98W8Fk7ypu0ePk+vbvnUAm7cWG2vjkrkxIGAAB8Vlr0oQJWVt9OGfs0Y8xDkh6SpNzcXNNX3xfwNVvLGnXfe/mUMAAAgM/4uICVN3TYnKRvME0RGASMMVpbeFD3v1+gNQUHKWEAAACfIzEySA5LKmvwjomKlDHARsYYvZdXo/vfL9CW0kbFhQfqV6fl6BszMhQWyP+eAAAAn+bvdCgpMthrxtv3ZrT905LmS4q1LKtc0q3GmH/2dzDAm7k9Rm/trNT97+/XnspmpUQF6/ZzxuvCqakK8uecMAAAgC+SEh2sA02ddsfoE72ZpnjpQAQBfEGP26OXP6rQAyv2q7CuTVlxofrLhRN19qRkrxjPCgAA0N8SIoK0o7zR7hh9gn1QwADo7HHruU1levCDQlU0dmhsUoT+ftkUnTwuUU6HZXc8AACAISMxIlDvNHfKGCPLGtrvoyhjQD9qaOvWE2tL9MTaYh1s61ZuRrR+d+54zc+OG/I/PAAAAOyQEBGkzh6Pmjtdigwe2oPOKGNAPyirb9cjqwr17KYydfZ4dEJOvBbNzdKM4TGUMAAAgKOQEBEkSapu7qSMAfg/28sb9eDKQr21o1JOh6VzJqVo0dwsjUoItzsaAACAV/h0Gcse4u+xKGPAUTLGaMXeWj24cr/WFdYrPMhPi+aO0FWzMz/5YQEAAIC+kRARKEmqbu6yOcnRo4wBX1O3y6PXth3QQyv3a191q5Iig3Tz6WN08bQ0DmoGAADoJ59eGRvqKGPAEapr7dJ/1pXqyfUlqm3pUk5iuO6+aKLOnMh4egAAgP4W5O9UZLA/ZQzwJbsONOmxNcV6besBdbs9mj86TlfPHq45o2IZygEAADCA4sMDVcM2RcC7uT1G7+yu1qNrirShqF7B/k5dPC1N35qdqRFxYXbHAwAA8EkxoQGqb++2O8ZRo4wBn6Opo0fPbyrT4x8Wq7yhQylRwfrVaTm6ODddkSHcDwYAAGCnmNAAFdS02h3jqFHGgE/Jr27Rv9eV6IXN5Wrvdmt6ZoxuPn2MFoxJkB/3gwEAAAwKMaEBqm9jZQwY8rpdHr29q0pPrivR+qJ6BTgdOmNikq6ePVzjUyLtjgcAAIDPiAkNUEN7tzweI4dj6N67TxmDz6po7NDT60v1zMYy1bV2KS0mWL88NUcXTk3VsLBAu+MBAADgC0SHBMhjpObOHkWFBNgd52ujjMGneDxGqwrq9OS6Ei3fUy0j6YTR8bp8ZobmZsfJOYSvrAAAAPiKmNBDBay+rZsyBgx2B1u79OKWcv1nfalKDrZrWGiAvjtvhC6dnq60mBC74wEAAOAIRH+qjGXF2RzmKFDG4LXcHqNV+bV6blOZ3tldrR630bTMaP3kpGydMj5RgX5OuyMCAADgaxj2qTI2lFHG4HXKG9r1/KZyPb+pTAeaOhUd4q8rj83UxdPSlJ0Qbnc8AAAAHKWPV8YahvhZY5QxeIUul1vv7K7WsxvLtLqgTpJ03MhY3XT6WC0YG88qGAAAgBeJCDpUY1o6XTYnOTqUMQxZxhjtrmzWi5sr9PJH5Wpo71FKVLB+dMIoXZibqtRo7gUDAADwRqEBh2pMM2UMGFhVTZ16ZWuFXt5Sob3VLfJ3Wlo4NlEXTUvTcSNjmYgIAADg5RwOS2GBfmqljAH9r63Lpbd3VemlLRVas79OxkhT0qN0+znjdcaEpE/2DQMAAMA3hAX6qbWrx+4YR4UyhkHL7TH6cH+dXtpSoaU7q9TR41ZaTLCuO2GUzp2couGxoXZHBAAAgE3CgvzU2sXKGNBnPB6jLaUNen17pd7YUanali6FB/npnMkpOm9KinIzomVZbEMEAADwdWGBfgzwAI6WMUbby5v0+vYDemN7pQ40dSrAz6ETRsfrrEnJOiEnXkH+TEMEAADA/wlnZQz4ej6ehPj69kq9sb1SpfXt8ndamjsqTj87ZbQWjElQeJC/3TEBAAAwSIUF+qmqqdPuGEeFMoYB4/EYba9o0rJdVVq6q0qFtW1yOizNHhmrH54wUiePTVRkCAUMAAAAX+3QAA9WxoAv1OP2aH1hvZbtrtKyXdWqau6U02FpZlaMrjkuS6eMT1QMkxABAABwhEIpY8D/au92aeW+Oi3bVaXleTVq6uhRkL9D87Lj9PNxo3VCTryiQihgAAAA+PoC/R3qdnnsjnFUKGPoE+UN7Xp/b61W5NVozf46dfZ4FBXirwVjEnTyuATNGRWn4ACGcAAAAKBvBPo51eXyyBgzZKdtU8bwtfS4PdpU3KAVe2v0Xl6N8mtaJUnpMSG6ZFq6Fo5N0PThMfJzOmxOCgAAAG8U6HfofWaXyzNkJ29TxtBrVU2dWplfqxV7a7RqX51aulzyd1qaMXyYLp6WpuNz4pUVGzpkr0wAAABg6Pi4gFHG4JWaOnq0rvCgPiyo0+qCOu2vbZMkJUYE6YyJSZo/Ol6zR8YqLJC/RgAAABhYn6yM9bil4KE5kZt30fhEl8utzSUNWlNQpzUFB7W9vFEeIwX7OzV9eIwumZau2SNjNSYpnNUvAAAA2OrT2xSHKsqYD2vtcmlLSYM2FtdrQ1G9tpY1qsvlkdNhaWJqpH54/EjNHhmryenRCvDj3i8AAAAMHv+3TdFtc5KvjzLmQw62dmlj8f+Vr92VzXJ7jByWNDY5Qt+Yka7ZI2I1IytG4UFDc6kXAAAAvuHjlbHOHlbGMMh0udzKq2zR1rJGbStr1NayRhXWHbrnK8DPoUlpUfrevBGaNjxGU9KjKF8AAAAYUvwPT+12eYzNSb4+ypgX8HiMig+2aVt5o7aVNemjskbtOdCsbvehqwSxYYGalBalC3JTNT0zRhNSIxXoNzQnzgAAAACS5HQcmmHg9rAyhgHS0e1WXlWz9lS2aHdlk/ZUtiivsllt3Yf2ygb7OzUhNVJXzc7UxLQoTUyLUnJkEAM3AAAA4FX8Dpcxl5uVMfSxbpdHxQfbVFDTqv01rdpb3aLdlc0qrmvTxyux4YF+GpMUoQumpmpscoSOSY3SqPgwDloGAACA1/t4ZYxtivhajDGqbe1SWX27CmvbVFDbqv01bdpf26rS+na5P/UXKy0mWGMSI3TmMckamxyhsUkRSo0OZsULAAAAPsnPV+4ZsyzrFEmLJTklPWKM+WO/pvISxhjVt3WrqrlTFQ0dKq1vV/nhx7L6dpU1tP/X9Bd/p6XhsaEakxSuM45J0sj4MI2IC1NWXKhCAujNAAAAwMfSooN1yxljlRUbaneUr+0r3+FbluWUdL+kkySVS9poWdZrxpjd/R1uMPJ4jFo6Xapv71ZDe7ca27vV0Naj2tYuVTV1qqalU1VNnapu7lJtS9cnQzQ+Fhbop7SYEA2PDdW87DilxYQoPSZEGcMOPbLFEAAAAPhq8RFB+vZxw+2OcVR6s9wyXVKBMaZQkizLekbS2ZKGTBlr73ZpdX6d3B6jHo+R2+NRj9vI7TFyuT1yeYy6XR61dbvV0e1SW7db7V2HH7tdautyq7mzR43tPWps79YXrYSGB/opPiJQCRFBmjE8RvERQUo8/HlKdLDSokMUFeLP1kIAAAAAvSpjKZLKPvV5uaQZn32SZVmLJC2SpPT09D4J11fqWrq16N+be/XckACnQgL8FBp4+DHAqYhgf6VEBSsqxF8xoQGKCglQdIi/okMCFHX4MS48UKGBbCUEAAAA0Du9aQ+ft4zzP2tDxpiHJD0kSbm5uYPqLrqEyEC9ft1x8nNa8nNY8nM45HRY8nd+/Hjo42B/pxwOVq0AAAAA9L/elLFySWmf+jxV0oH+idM/Av2cGp8SaXcMAP+/vbuLkaus4zj+/YUCiWgERBEBkRjihUYRG6ohmiZoLYSAGF9KjOJbFBUjFyb4kiipN/iCiXqhUWmCBhFF0V4UoYkmXkEoTRVqUaopWmhataZIMDHVvxdz2kzGOdstnZlzmPl+ks3unuc5m//Mf5+3Pc85K0mSpMOW87SI+4Hzkpyb5ARgHbBxumFJkiRJ0nw74pWxqjqY5FrgbgaPtt9QVdunHpkkSZIkzbFlPXGiqjYBm6YciyRJkiQtDP+plSRJkiR1wMWYJEmSJHXAxZgkSZIkdcDFmCRJkiR1wMWYJEmSJHXAxZgkSZIkdcDFmCRJkiR1IFU1+R+a/BV4dOI/+Ok7Dfhb10EIMBd9Yi76w1z0h7noD3PRH+aiP8xFfywnF+dU1fOXqjCVxVjfJNlSVSu7jkPmok/MRX+Yi/4wF/1hLvrDXPSHueiPSeXCbYqSJEmS1AEXY5IkSZLUgUVZjH276wB0mLnoD3PRH+aiP8xFf5iL/jAX/WEu+mMiuViIe8YkSZIkqW8W5cqYJEmSJPXKXC3GkqxN8vskO5N8akz5iUlub8rvS/KS2Uc5/5KcneRXSXYk2Z7kE2PqrE5yIMm25uNzXcS6CJLsSvJg8z5vGVOeJF9v2sVvk1zQRZzzLsnLhn7ftyV5Isl1I3VsF1OSZEOSfUkeGjp2apLNSR5pPp/Scu7VTZ1Hklw9u6jnU0suvpzk4aYPujPJyS3nLtmf6ei05OKGJI8N9UOXtpy75JxLR6clF7cP5WFXkm0t59ouJqhtHjutMWNutikmOQ74A/AmYDdwP3BVVf1uqM5HgVdW1TVJ1gFXVtU7Owl4jiU5AzijqrYmeQ7wAPCWkVysBj5ZVZd1FObCSLILWFlVY/8XRjPQfhy4FFgFfK2qVs0uwsXT9FePAauq6tGh46uxXUxFkjcATwLfq6pXNMe+BOyvqhubyeQpVXX9yHmnAluAlUAx6M9eU1X/mOkLmCMtuVgD/LKqDib5IsBoLpp6u1iiP9PRacnFDcCTVfWVJc474pxLR2dcLkbKbwIOVNX6MWW7sF1MTNs8FngvUxgz5unK2IXAzqr6U1X9G/ghcMVInSuAW5qv7wAuTpIZxrgQqmpPVW1tvv4nsAM4s9uotIQrGHT+VVX3Aic3HZGm52Lgj8MLMU1XVf0a2D9yeHhMuIXBYDvqzcDmqtrfDKabgbVTC3QBjMtFVd1TVQebb+8Fzpp5YAuopV0sx3LmXDoKS+Wimau+A7htpkEtqCXmsVMZM+ZpMXYm8Jeh73fz/wuAw3WaTv8A8LyZRLegMtgK+mrgvjHFr0vymyR3JXn5TANbLAXck+SBJB8aU76ctqPJWkf7oGq7mJ3Tq2oPDAZf4AVj6tg+Zu/9wF0tZUfqzzQZ1zZbRje0bMWyXczW64G9VfVIS7ntYkpG5rFTGTPmaTE27grX6B7M5dTRhCR5NvAT4LqqemKkeCtwTlW9CvgG8LNZx7dALqqqC4BLgI81WyGG2S5mKMkJdfh5WwAAAqlJREFUwOXAj8cU2y76x/YxQ0k+CxwEbm2pcqT+TMfum8BLgfOBPcBNY+rYLmbrKpa+Kma7mIIjzGNbTxtzbMm2MU+Lsd3A2UPfnwU83lYnyQrguTy9y/M6giTHM/gFvrWqfjpaXlVPVNWTzdebgOOTnDbjMBdCVT3efN4H3Mlge8mw5bQdTc4lwNaq2jtaYLuYub2HtuQ2n/eNqWP7mJHmRvfLgHdVyw3ty+jPdIyqam9V/aeq/gt8h/Hvse1iRpr56luB29vq2C4mr2UeO5UxY54WY/cD5yU5t/nL8zpg40idjcChp5q8jcHNwv4lZ8Kavc03Azuq6qstdV546H69JBcy+F38++yiXAxJTmpuPiXJScAa4KGRahuB92TgtQxuEN4z41AXSetfOG0XMzc8JlwN/HxMnbuBNUlOabZrrWmOaYKSrAWuBy6vqqda6iynP9MxGrln+ErGv8fLmXNpMt4IPFxVu8cV2i4mb4l57FTGjBXHHnI/NE9gupbBCz4O2FBV25OsB7ZU1UYGb+z3k+xkcEVsXXcRz7WLgHcDDw49hvUzwIsBqupbDBbDH0lyEPgXsM6F8VScDtzZzO9XAD+oql8kuQYO52ITgycp7gSeAt7XUaxzL8mzGDx97MNDx4ZzYbuYkiS3AauB05LsBj4P3Aj8KMkHgD8Db2/qrgSuqaoPVtX+JF9gMPkEWF9V7qg4Bi25+DRwIrC56a/ubZ58/CLgu1V1KS39WQcvYW605GJ1kvMZbK3aRdNfDeeibc7VwUuYG+NyUVU3M+YeY9vF1LXNY6cyZszNo+0lSZIk6ZlknrYpSpIkSdIzhosxSZIkSeqAizFJkiRJ6oCLMUmSJEnqgIsxSZIkSeqAizFJkiRJ6oCLMUmSJEnqgIsxSZIkSerA/wCjPbkTb9u8hQAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 1080x360 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"q2_array=scalar_product_4(q_np,q_np)\n",
"q2_range=np.arange(4*mmu_mass**2, (mB_mass-mKst_mass)**2,1000)\n",
"\n",
"val=(q2_range/(1-beta_l(q2_range,mmu_mass)))*np.sqrt(lambda_function(mB_mass, mKst_mass, q2_range))*beta_l(q2_range, mmu_mass)\n",
"plt.plot(q2_range/1e6, val);\n",
"fig=plt.gcf()\n",
"fig.set_size_inches(15,5)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#class dGamma(zfit.pdf.ZPDF):\n",
"# \n",
"# _PARAMS = ['ml','mB','mKst']\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",
"# pKstx, pKsty, pKstz, pKstE, qx, qy, qz, qE, p1x, p1y, p1z, p1E, p2x, p2y, p2z, p2E = x.unstack_x()\n",
"# \n",
"# pKst = np.array([pKstx, pKsty, pKstz, pKstE])\n",
"# q = np.array([qx, qy, qz, qE])\n",
"# p1 = np.array([p1x, p1y, p1z, p1E])\n",
"# p2 = np.array([p2x, p2y, p2z, p2E])\n",
"# \n",
"# #z = tf.expand_dims(tf.stack([0., 0., 1., 0.], axis=0), axis=0)\n",
"# #z = [0.,0.,1.]\n",
"# q=p1+p2\n",
"# q2=scalar_product_4(q,q)\n",
"# modq = ztf.sqrt(q2)\n",
"#\n",
"# #q_in_q_rf = lorentz_vector([0.,0.,0.], modq)\n",
"# #p1q = lorentz_boost(p1, q_in_q_rf)\n",
"# #costheta_l= get_costheta_l(p1, z)\n",
"# \n",
"# pdf= lambda_function(mB, mKst, q2)*beta_l(q2, ml)\n",
"# return pdf\n",
"#"
]
},
{
"cell_type": "code",
"execution_count": 138,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEFCAYAAAD36MwKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAKc0lEQVR4nO3dX4ymZ1nH8d9FV0FQB+MSE/uHqZYUN0ZEJ4iQmEZMbFMXIsSE+ifRNOyJIBqJLp5pYuwBMZqAmg1WYiQlWnvg2gZMENMeEEILoq21pimFbotpCTL+OakNlwfzAtt1uzPZGfZ55+rnc7Izz+47ufJk97v33u8z91Z3B4A5XrD0AAAcLGEHGEbYAYYRdoBhhB1gmCNLD5AkR48e7c3NzaXHADhU7rvvvi9298vOvb4WYd/c3My999679BgAh0pVfe58123FAAwj7ADDLBr2qjpeVae2t7eXHANglEXD3t2nu/vExsbGkmMAjGIrBmAYYQcYRtgBhhF2gGHW4huU9mPz5J0X/dpHb7nxACcBWA9W7ADDCDvAMMIOMIywAwwj7ADDOCsGYBhnxQAMYysGYBhhBxhG2AGGEXaAYYQdYBhhBxhG2AGGEXaAYYQdYBhhBxhG2AGGEXaAYYQdYBhhBxhG2AGG8R9tAAzjP9oAGMZWDMAwwg4wjLADDCPsAMMIO8Awwg4wjLADDCPsAMMIO8Awwg4wjLADDCPsAMMIO8Awwg4wjLADDCPsAMMIO8Awwg4wjLADDCPsAMMsGvaqOl5Vp7a3t5ccA2CURcPe3ae7+8TGxsaSYwCMYisGYBhhBxhG2AGGEXaAYYQdYBhhBxhG2AGGEXaAYYQdYBhhBxhG2AGGEXaAYYQdYBhhBxhG2AGGEXaAYYQdYBhhBxhG2AGGEXaAYYQdYBhhBxhG2AGGEXaAYYQdYBhhBxhG2AGGEXaAYYQdYBhhBxhG2AGGEXaAYYQdYBhhBxjmwMNeVddV1T1V9SdVdd1Bf30ALmxPYa+qW6vqyaq6/5zr11fVQ1X1cFWdXF3uJP+d5EVJzhzsuADsZq8r9g8kuf7sC1V1WZL3JbkhybEkN1XVsST3dPcNSX4zyW8f3KgA7MWewt7ddyf50jmXX5Pk4e5+pLufTvKhJG/q7q+sfv4/krzwwCYFYE+O7OO1lyd57KzPzyT5kap6c5KfTPLSJO99rhdX1YkkJ5Lkqquu2scYAJxtP2Gv81zr7r4jyR27vbi7TyU5lSRbW1u9jzkAOMt+noo5k+TKsz6/IskT+xsHgP3az4r9k0leUVVXJ3k8yVuT/OyBTHWJbJ6886Jf++gtNx7gJAAHZ6+PO96W5ONJrq2qM1V1c3c/k+TtST6S5MEkf9ndD3zjRgVgL/a0Yu/um57j+l1J7jrQiQDYl0WPFKiq41V1ant7e8kxAEZZNOzdfbq7T2xsbCw5BsAoDgEDGEbYAYYRdoBhhB1gGE/FAAzjqRiAYWzFAAwj7ADDCDvAMMIOMIywAwwj7ADDeI4dYBjPsQMMYysGYBhhBxhG2AGG2dP/ecr/t3nyzn29/tFbbjygSQCezYodYBhhBxjGc+wAw3iOHWAYWzEAwwg7wDDCDjCMsAMMI+wAwwg7wDDCDjCMsAMMI+wAwzhSAGCYRY/t7e7TSU5vbW29bck5lrCfY38d+QtciK0YgGGEHWAYYQcYRtgBhhF2gGGEHWAYYQcYRtgBhhF2gGGEHWCYRY8UqKrjSY5fc801S45x6DiOALiQRVfs3X26u09sbGwsOQbAKLZiAIYRdoBhhB1gGGEHGEbYAYYRdoBhFn2OnUvPM/AwnxU7wDDCDjCMsAMMI+wAwwg7wDDCDjCMxx3Zs/08Kpl4XBIulUVX7FV1vKpObW9vLzkGwCjOYwcYxh47wDDCDjCMsAMMI+wAw3jckUvGyZJwaVixAwwj7ADDCDvAMMIOMIw3TzkUvPEKe2fFDjCMsAMMI+wAw9hjZzz78zzfWLEDDCPsAMMIO8Awwg4wjDdP4QK88cphZMUOMIywAwyzaNir6nhVndre3l5yDIBRFt1j7+7TSU5vbW29bck54BthP/vziT16Lp6tGIBhhB1gGI87wpryqCUXy4odYBhhBxjGVgwMZBvn+c2KHWAYYQcYxlYM8Cy2cQ4/K3aAYYQdYBhhBxjGHjtwYBx8th6s2AGGEXaAYWzFAGvDo5YHQ9iBEfyl8HW2YgCGsWIHnvemrfat2AGGEXaAYWzFAOzDOm7jWLEDDCPsAMMIO8Awwg4wjLADDCPsAMMIO8Awwg4wjLADDFPdvfQMqaqnknzuIl9+NMkXD3Ccadyf3blHF+b+7G6pe/Ty7n7ZuRfXIuz7UVX3dvfW0nOsK/dnd+7Rhbk/u1u3e2QrBmAYYQcYZkLYTy09wJpzf3bnHl2Y+7O7tbpHh36PHYBnm7BiB+Aswg4wzKEOe1VdX1UPVdXDVXVy6XnWSVVdWVUfq6oHq+qBqnrn0jOto6q6rKo+XVV/u/Qs66iqXlpVt1fVv65+L/3o0jOtk6r6tdWfr/ur6raqetHSMyWHOOxVdVmS9yW5IcmxJDdV1bFlp1orzyT59e7+viSvTfLL7s95vTPJg0sPscb+MMmHu/uVSV4V9+prquryJL+SZKu7vz/JZUneuuxUOw5t2JO8JsnD3f1Idz+d5ENJ3rTwTGuju7/Q3Z9affxf2fkDefmyU62XqroiyY1J3r/0LOuoqr49yY8l+dMk6e6nu/vLy061do4k+ZaqOpLkxUmeWHieJIc77Jcneeysz89EuM6rqjaTvDrJJ5adZO38QZLfSPKVpQdZU9+T5Kkkf7barnp/Vb1k6aHWRXc/nuQ9ST6f5AtJtrv775adasdhDnud55pnN89RVd+a5K+T/Gp3/+fS86yLqvqpJE92931Lz7LGjiT5oSR/3N2vTvI/SbyXtVJV35GdXYKrk3x3kpdU1c8vO9WOwxz2M0muPOvzK7Im/wxaF1X1TdmJ+ge7+46l51kzr0/yxqp6NDvbeD9eVX+x7Ehr50ySM9391X/p3Z6d0LPjJ5J8truf6u7/TXJHktctPFOSwx32TyZ5RVVdXVXfnJ03Lf5m4ZnWRlVVdvZGH+zu3196nnXT3e/u7iu6ezM7v3f+vrvXYrW1Lrr735M8VlXXri69Icm/LDjSuvl8ktdW1YtXf97ekDV5c/nI0gNcrO5+pqrenuQj2Xk3+tbufmDhsdbJ65P8QpJ/rqp/XF37re6+a8GZOHzekeSDq8XTI0l+aeF51kZ3f6Kqbk/yqew8hfbprMnRAo4UABjmMG/FAHAewg4wjLADDCPsAMMIO8AlVlW3VtWTVXX/Hn7ty6vqo1X1T1X1D6ujMC5I2AEuvQ8kuX6Pv/Y9Sf68u38gye8k+b3dXiDsAJdYd9+d5EtnX6uq762qD1fVfVV1T1W9cvVTx5J8dPXxx7KHww6FHWA9nEryju7+4STvSvJHq+ufSfKW1cc/neTbquo7L/SFDu13ngJMsTqs73VJ/mrndIIkyQtXP74ryXur6heT3J3k8ex8p+tzEnaA5b0gyZe7+wfP/YnufiLJm5Ov/QXwlu7e3u2LAbCg1ZHan62qn0l2DvGrqletPj5aVV9t9buT3Lrb1xN2gEusqm5L8vEk11bVmaq6OcnPJbm5qj6T5IF8/U3S65I8VFX/luS7kvzurl/fIWAAs1ixAwwj7ADDCDvAMMIOMIywAwwj7ADDCDvAMP8HJcE65xtjqmgAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.hist(dGamma,bins = 20, log=True);"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.8"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Davide Lancierini