{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from math import pi\n",
"import matplotlib\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"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"
]
},
{
"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 tensorflow as tf\n",
"import zfit\n",
"\n",
"ztf = zfit.ztf\n",
"ztyping = zfit.util.ztyping\n",
"ztypes = zfit.settings.ztypes"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"zfit.run.check_numerics = False\n",
"zfit.settings.set_verbosity(10)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def beta_l(q2, m):\n",
" beta_l = ztf.sqrt(1-4*(ztf.square(m)/q2))\n",
" return beta_l\n",
"\n",
"def lambda_phase(m1, m2, q2):\n",
" \n",
" lambd = ztf.square(ztf.square(m1)-ztf.square(m2)) -2*q2*(ztf.square(m1)+ztf.square(m2)) + ztf.square(q2)\n",
" \n",
" return lambd\n",
"\n",
"def beta_l_np(q2, m):\n",
" beta_l = np.sqrt(1-4*(np.square(m)/q2))\n",
" return beta_l\n",
"\n",
"def lambda_phase_np(m1, m2, q2):\n",
" \n",
" lambd = np.square(np.square(m1)-np.square(m2)) -2*q2*(np.square(m1)+np.square(m2)) + np.square(q2)\n",
" \n",
" return lambd\n",
"def scalar_product_3_np(v1,v2):\n",
" \n",
" output = v1[:,1]*v2[:,1]+v1[:,2]*v2[:,2]+v1[:,3]*v2[:,3]\n",
" return output\n",
"\n",
"def scalar_product_4_np(v1,v2):\n",
" \n",
" scalar_3 = scalar_product_3_np(v1,v2)\n",
" output = v1[:,0]*v2[:,0]- scalar_3\n",
" return output\n",
"\n",
"\n",
"def a(q2, m, m1, m2):\n",
" \n",
" step1=q2*(1+ztf.square(beta_l(q2, m)))+ lambda_phase(m1, m2, q2)/2 + 2*m*(ztf.square(m1)-ztf.square(m2)+q2)+4*ztf.square(m)*ztf.square(m1)\n",
" #step2 = beta_l(q2, m)*ztf.sqrt(lambda_phase(m1, m2, q2))*step1\n",
" \n",
" return step1\n",
"\n",
"def b(q2, m, m1, m2):\n",
" \n",
" step1=2*(q2*(1+ztf.square(beta_l(q2, m)))+ m*(ztf.sqrt(lambda_phase(m1, m2, q2))*beta_l(q2,m)+(ztf.square(m1)-ztf.square(m2)+q2)))\n",
" #step2 = beta_l(q2, m)*ztf.sqrt(lambda_phase(m1, m2, q2))*step1\n",
" \n",
" return step1\n",
"\n",
"def c(q2, m, m1, m2):\n",
" \n",
" step1= q2*(1+ztf.square(beta_l(q2, m)))- lambda_phase(m1, m2, q2)/2 + 2*m*ztf.sqrt(lambda_phase(m1, m2, q2))*beta_l(q2, m)\n",
" #step2 = beta_l(q2, m)*ztf.sqrt(lambda_phase(m1, m2, q2))*step1\n",
" \n",
" return step1\n",
"\n",
"\n",
"def matrix_elt(q2, m, m1, m2, cos_theta_l):\n",
" \n",
" out = a(q2, m, m1, m2)+b(q2, m, m1, m2)*cos_theta_l+c(q2, m, m1, m2)*ztf.square(cos_theta_l)\n",
" \n",
" return out \n",
" \n",
"def phase_space(q2, m, m1, m2):\n",
" \n",
" phase_space_term = beta_l(q2, m)*ztf.sqrt(lambda_phase(m1, m2, q2))\n",
" \n",
" return phase_space_term\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"mmu_mass = 0.105\n",
"mB_mass=5.280\n",
"mKst_mass=0.892\n",
"\n",
"lower_q2= 4*mmu_mass*mmu_mass\n",
"upper_q2=(mB_mass-mKst_mass)*(mB_mass-mKst_mass)\n",
"cutoff=0.01"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"class dGamma(zfit.pdf.ZPDF):\n",
" \n",
" _PARAMS = ['mB','mKst','ml']\n",
" _N_OBS = 2\n",
" VALID_LIMITS={\n",
" 'cos_theta_l':(-1,1),\n",
" 'q2': (lower_q2, upper_q2),\n",
" #'Eg': (cutoff, 10)\n",
" }\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",
" q2, cos_theta_l = x.unstack_x()\n",
" pdf = matrix_elt(q2, ml, mB, mKst, cos_theta_l)*phase_space(q2, ml, mB, mKst)\n",
" #pdf= phase_space(q2, mB, mKst, ml)\n",
" \n",
" return pdf\n"
]
},
{
"cell_type": "code",
"execution_count": 7,
"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": [
"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": 8,
"metadata": {},
"outputs": [],
"source": [
"lower = ((lower_q2, -1.),)\n",
"upper = ((upper_q2, 1.),)\n",
"\n",
"obs = zfit.Space([\"q2\",\"cos_theta_l\"], limits=(lower,upper))\n"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"pdf = dGamma(obs=obs, ml = mlepton_par, mB = mB_par, mKst = mKst_par)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"n_events=100000"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"WARNING:tensorflow:From /Users/davide/miniconda3/envs/zfit_env/lib/python3.6/site-packages/zfit/core/sample.py:163: to_int32 (from tensorflow.python.ops.math_ops) is deprecated and will be removed in a future version.\n",
"Instructions for updating:\n",
"Use tf.cast instead.\n"
]
}
],
"source": [
"sampler = pdf.create_sampler(n=n_events)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"for i in range(1):\n",
" sampler.resample()\n"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"a=sampler.to_pandas()"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"q2vals_MC = np.array([a['q2'][i] for i in range(n_events)])\n",
"q2vals = q2vals_MC.reshape(n_events,1)\n",
"\n",
"cos_theta_vals_MC = np.array([a['cos_theta_l'][i] for i in range(n_events)])\n",
"cos_theta_vals = cos_theta_vals_MC.reshape(n_events,1)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA3YAAAEvCAYAAAAJs1ObAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAYEklEQVR4nO3df6yldZ0f8PdnQU2qm4oFKQu4Q8y0KTZdtBOwMW1o7CI/Nos2dQtplLU2oy00mthEtEkxGJvZdnWjjaXBdbLYuCKtWidxdnFqNWb/wGUgVEW0TN1ZGZnAuBjU2OwG/PSP++AeLvfX3HvPuee55/VKbu693+f73PM98/Ccc958f1V3BwAAgPH6hZ1uAAAAAFsj2AEAAIycYAcAADBygh0AAMDICXYAAAAjJ9gBAACM3Jk73YC1nH322b1nz56dbgYAAMCOuO+++37Q3eesV2+ug92ePXty9OjRnW4GAADAjqiqP91IPUMxAQAARk6wAwAAGDnBDgAAYOQEOwAAgJET7AAAAEZOsAMAABg5wQ4AAGDkBDsAAICRE+wAAABGbt1gV1UXVtWXq+qhqnqwqt4xlL+vqr5fVQ8MX1dPnPOeqjpWVd+pqtdNlF85lB2rqpun85QAAAAWy5kbqPNUknd19/1V9YtJ7quqI8Ox3+nu356sXFUXJ7kuySuS/FKS/1lVf2M4/NEkv5rkRJJ7q+pQd39rO54IAADAolo32HX3ySQnh59/XFUPJTl/jVOuTXJnd/95kj+pqmNJLh2OHevu7yZJVd051F24YLfn5i+seuz4gWtm2BIAAGA3OK05dlW1J8krk3xtKLqpqr5eVQer6qyh7Pwkj0ycdmIoW60cAACALdhwsKuqFyX5TJJ3dvePktyW5OVJLslSj94Hn6m6wum9Rvnyx9lfVUer6uipU6c22jwAAICFtaFgV1XPy1Ko+2R3fzZJuvux7n66u3+W5GP5y+GWJ5JcOHH6BUkeXaP8Wbr79u7e1937zjnnnNN9PgAAAAtnI6tiVpKPJ3mouz80UX7eRLU3JPnm8POhJNdV1Quq6qIke5P8cZJ7k+ytqouq6vlZWmDl0PY8DQAAgMW1kVUxX5PkTUm+UVUPDGXvTXJ9VV2SpeGUx5O8LUm6+8GquitLi6I8leTG7n46SarqpiR3JzkjycHufnAbnwsAAMBC2siqmH+UlefHHV7jnA8k+cAK5YfXOg8AAIDTd1qrYgIAADB/BDsAAICRE+wAAABGTrADAAAYOcEOAABg5Day3QGnac/NX9ixv338wDVTe2wAAGA+6bEDAAAYOcEOAABg5AzFnDPTHMYJAADsTnrsAAAARk6wAwAAGDnBDgAAYOTMsVsgtkoAAIDdSY8dAADAyAl2AAAAIyfYAQAAjJw5dvzcWnPwzL8DAID5JdjtMjY4BwCAxWMoJgAAwMjpsWNDbJUAAADzS48dAADAyAl2AAAAIyfYAQAAjJxgBwAAMHKCHQAAwMhZFZNtYdVMAADYOXrsAAAARk6wAwAAGDnBDgAAYOQEOwAAgJET7AAAAEZOsAMAABg52x2w42yVAAAAWyPYMRPrhTcAAGDzBDtGTW8fAACYYwcAADB6gh0AAMDICXYAAAAjJ9gBAACMnGAHAAAwclbF3ARL9wMAAPNEjx0AAMDI6bFj7ukhBQCAta3bY1dVF1bVl6vqoap6sKreMZS/pKqOVNXDw/ezhvKqqo9U1bGq+npVvWrib90w1H+4qm6Y3tMCAABYHBsZivlUknd1999K8uokN1bVxUluTvKl7t6b5EvD70lyVZK9w9f+JLclS0EwyS1JLktyaZJbngmDAAAAbN66wa67T3b3/cPPP07yUJLzk1yb5I6h2h1JXj/8fG2ST/SSe5K8uKrOS/K6JEe6+4nu/mGSI0mu3NZnAwAAsIBOa/GUqtqT5JVJvpbk3O4+mSyFvyQvHaqdn+SRidNODGWrlS9/jP1VdbSqjp46dep0mgcAALCQNhzsqupFST6T5J3d/aO1qq5Q1muUP7ug+/bu3tfd+84555yNNg8AAGBhbWhVzKp6XpZC3Se7+7ND8WNVdV53nxyGWj4+lJ9IcuHE6RckeXQov3xZ+Vc233RY33orah4/cM2MWgIAANOzkVUxK8nHkzzU3R+aOHQoyTMrW96Q5PMT5W8eVsd8dZInh6Gadye5oqrOGhZNuWIoAwAAYAs20mP3miRvSvKNqnpgKHtvkgNJ7qqqtyb5XpI3DscOJ7k6ybEkP03yliTp7ieq6v1J7h3q3drdT2zLswAAAFhg6wa77v6jrDw/Lkleu0L9TnLjKn/rYJKDp9NAAAAA1rahOXawW601B8/8OwAAxuK0tjsAAABg/gh2AAAAIyfYAQAAjJw5drAKe+ABADAWeuwAAABGTrADAAAYOcEOAABg5Myxg02yBx4AAPNCjx0AAMDICXYAAAAjJ9gBAACMnGAHAAAwcoIdAADAyFkVE6ZgrRUzE6tmAgCwvfTYAQAAjJxgBwAAMHKGYsIOMFQTAIDtpMcOAABg5AQ7AACAkTMUE0bGME4AAJbTYwcAADBygh0AAMDICXYAAAAjZ44dzKH15tEBAMAkPXYAAAAjJ9gBAACMnGAHAAAwcoIdAADAyAl2AAAAIyfYAQAAjJxgBwAAMHL2sYNdZr098I4fuGZGLQEAYFb02AEAAIycYAcAADBygh0AAMDICXYAAAAjJ9gBAACMnFUxYcGstWrmeitmbuVcAACmR48dAADAyAl2AAAAI2coJvBz621uDgDAfFq3x66qDlbV41X1zYmy91XV96vqgeHr6olj76mqY1X1nap63UT5lUPZsaq6efufCgAAwGLayFDM30ty5Qrlv9Pdlwxfh5Okqi5Ocl2SVwzn/OeqOqOqzkjy0SRXJbk4yfVDXQAAALZo3aGY3f3Vqtqzwb93bZI7u/vPk/xJVR1Lculw7Fh3fzdJqurOoe63TrvFAAAAPMtW5tjdVFVvTnI0ybu6+4dJzk9yz0SdE0NZkjyyrPyyLTw2MGfWm59nOwQAgOnZ7KqYtyV5eZJLkpxM8sGhvFao22uUP0dV7a+qo1V19NSpU5tsHgAAwOLYVLDr7se6++nu/lmSj+Uvh1ueSHLhRNULkjy6RvlKf/v27t7X3fvOOeeczTQPAABgoWwq2FXVeRO/viHJMytmHkpyXVW9oKouSrI3yR8nuTfJ3qq6qKqen6UFVg5tvtkAAAA8Y905dlX1qSSXJzm7qk4kuSXJ5VV1SZaGUx5P8rYk6e4Hq+quLC2K8lSSG7v76eHv3JTk7iRnJDnY3Q9u+7MBAABYQBtZFfP6FYo/vkb9DyT5wArlh5McPq3WAQAAsK7NLp4CAADAnBDsAAAARm4r+9gBbJh97gAApkewA+bCWsFP6AMAWJtgB4yeUAgALDpz7AAAAEZOjx0w99abnwcAsOj02AEAAIycYAcAADBygh0AAMDImWMH7Gpb3T/PipsAwBjosQMAABg5PXbAQrPiJgCwG+ixAwAAGDnBDgAAYOQEOwAAgJET7AAAAEZOsAMAABg5wQ4AAGDkBDsAAICRE+wAAABGzgblADtgvY3Rjx+4ZkYtAQB2A8EOYErWC28AANvFUEwAAICRE+wAAABGzlBMgE0y1BIAmBd67AAAAEZOsAMAABg5wQ4AAGDkBDsAAICRE+wAAABGTrADAAAYOcEOAABg5AQ7AACAkRPsAAAARu7MnW4AANtrz81f2PS5xw9cs40tAQBmRbADmEPrhTMBDACYZCgmAADAyAl2AAAAIyfYAQAAjJw5dgAjtJUFUgCA3UePHQAAwMitG+yq6mBVPV5V35woe0lVHamqh4fvZw3lVVUfqapjVfX1qnrVxDk3DPUfrqobpvN0AAAAFs9Geux+L8mVy8puTvKl7t6b5EvD70lyVZK9w9f+JLclS0EwyS1JLktyaZJbngmDAAAAbM26c+y6+6tVtWdZ8bVJLh9+viPJV5K8eyj/RHd3knuq6sVVdd5Q90h3P5EkVXUkS2HxU1t+BgDMjbXm/tl7DwCmZ7Nz7M7t7pNJMnx/6VB+fpJHJuqdGMpWKwcAAGCLtntVzFqhrNcof+4fqNqfpWGcednLXrZ9LQNgy6zGCQDzabPB7rGqOq+7Tw5DLR8fyk8kuXCi3gVJHh3KL19W/pWV/nB3357k9iTZt2/fiuEPgOkQ3ABgnDYb7A4luSHJgeH75yfKb6qqO7O0UMqTQ/i7O8m/n1gw5Yok79l8swEYm/VCozl4ALB56wa7qvpUlnrbzq6qE1la3fJAkruq6q1JvpfkjUP1w0muTnIsyU+TvCVJuvuJqnp/knuHerc+s5AKACSCHwBsxUZWxbx+lUOvXaFuJ7lxlb9zMMnB02odAAAA69rsqpgAAADMCcEOAABg5AQ7AACAkRPsAAAARk6wAwAAGLnN7mMHADO11nYItkIAYNEJdgCMnj3wAFh0hmICAACMnB47AHY9wzgB2O302AEAAIycYAcAADBygh0AAMDImWMHAGuw4iYAYyDYAbDQ1gtuADAGhmICAACMnB47ANgCWykAMA/02AEAAIycYAcAADBygh0AAMDICXYAAAAjZ/EUAJgSe+ABMCt67AAAAEZOjx0A7JCtbI6utw+ASXrsAAAARk6PHQAsGHP/AHYfwQ4ARkg4A2CSoZgAAAAjJ9gBAACMnGAHAAAwcubYAcAutJWtFAAYHz12AAAAIyfYAQAAjJyhmADAs6w1jNM2CgDzSbADADZsq/vnCY0A02EoJgAAwMgJdgAAACMn2AEAAIycOXYAwLbZyv55W52/B7DI9NgBAACMnGAHAAAwcoIdAADAyJljBwCMnvl5wKLTYwcAADByW+qxq6rjSX6c5OkkT3X3vqp6SZJPJ9mT5HiS3+juH1ZVJflwkquT/DTJb3b3/Vt5fABgcWxlxU2A3W47hmL+w+7+wcTvNyf5UncfqKqbh9/fneSqJHuHr8uS3DZ8BwCYKkM1gd1uGkMxr01yx/DzHUleP1H+iV5yT5IXV9V5U3h8AACAhbLVYNdJvlhV91XV/qHs3O4+mSTD95cO5ecneWTi3BND2bNU1f6qOlpVR0+dOrXF5gEAAOx+Wx2K+ZrufrSqXprkSFV9e426tUJZP6eg+/YktyfJvn37nnMcAACAZ9tSsOvuR4fvj1fV55JcmuSxqjqvu08OQy0fH6qfSHLhxOkXJHl0K48PADBt5ucBY7DpYFdVL0zyC9394+HnK5LcmuRQkhuSHBi+f3445VCSm6rqziwtmvLkM0M2AQB2khU3gbHbSo/duUk+t7SLQc5M8vvd/YdVdW+Su6rqrUm+l+SNQ/3DWdrq4FiWtjt4yxYeGwAAgMGmg113fzfJr6xQ/mdJXrtCeSe5cbOPBwAAwMq2Yx87AICFZQ4eMA+msY8dAAAAM6THDgBgh+jtA7aLYAcAMEVbWXFzrXPXC31CIywWQzEBAABGTrADAAAYOUMxAQBGyKbqwCQ9dgAAACMn2AEAAIycoZgAADyLFTVhfAQ7AIAFNK1tGNYjFMJ0GIoJAAAwcoIdAADAyBmKCQDAzJi/B9Mh2AEAMDfWCn5CH6xOsAMAYBT09sHqzLEDAAAYOcEOAABg5AzFBABg1zOMk91OsAMAYOHZdJ2xE+wAANgVthLOpvm4gh+zYI4dAADAyAl2AAAAIyfYAQAAjJxgBwAAMHIWTwEAgB1i4RW2i2AHAABTZLVOZsFQTAAAgJHTYwcAAHNqp3r7GB89dgAAACOnxw4AABaQOXi7i2AHAAA8x1rBT+ibP4ZiAgAAjJweOwAA4LQYxjl/9NgBAACMnB47AABgW21lmwa9fZsj2AEAAHPDMM/NMRQTAABg5AQ7AACAkTMUEwAAGA3z91amxw4AAGDkBDsAAICRm/lQzKq6MsmHk5yR5He7+8Cs2wAAACyetYZxjn2Y5kx77KrqjCQfTXJVkouTXF9VF8+yDQAAALvNrIdiXprkWHd/t7v/IsmdSa6dcRsAAAB2lVkHu/OTPDLx+4mhDAAAgE2a9Ry7WqGsn1Whan+S/cOvP6mq70y9Vafn7CQ/2OlGkMS1mCeuxfxwLeaHazE/XIv54VrMD9dimfqtHXvo9a7FL2/kj8w62J1IcuHE7xckeXSyQnffnuT2WTbqdFTV0e7et9PtwLWYJ67F/HAt5odrMT9ci/nhWswP12J+bNe1mPVQzHuT7K2qi6rq+UmuS3Joxm0AAADYVWbaY9fdT1XVTUnuztJ2Bwe7+8FZtgEAAGC3mfk+dt19OMnhWT/uNprbYaILyLWYH67F/HAt5odrMT9ci/nhWswP12J+bMu1qO5evxYAAABza9Zz7AAAANhmgt0KqurKqvpOVR2rqptXOP6Cqvr0cPxrVbVn9q1cDFV1YVV9uaoeqqoHq+odK9S5vKqerKoHhq9/txNtXQRVdbyqvjH8Ox9d4XhV1UeGe+PrVfWqnWjnbldVf3Piv/cHqupHVfXOZXXcF1NSVQer6vGq+uZE2Uuq6khVPTx8P2uVc28Y6jxcVTfMrtW70yrX4j9W1beH16DPVdWLVzl3zdczTs8q1+J9VfX9idehq1c5d83PXZyeVa7Fpyeuw/GqemCVc90X22i1z7HTes8wFHOZqjojyf9J8qtZ2p7h3iTXd/e3Jur8qyR/p7vfXlXXJXlDd//THWnwLldV5yU5r7vvr6pfTHJfktcvux6XJ/k33f1rO9TMhVFVx5Ps6+4V91oZ3rT/dZKrk1yW5MPdfdnsWrh4htes7ye5rLv/dKL88rgvpqKq/kGSnyT5RHf/7aHsPyR5orsPDB9Mz+rudy877yVJjibZl6U9XO9L8ne7+4czfQK7yCrX4ook/2tYsO23kmT5tRjqHc8ar2ecnlWuxfuS/KS7f3uN89b93MXpWelaLDv+wSRPdvetKxw7HvfFtlntc2yS38wU3jP02D3XpUmOdfd3u/svktyZ5Nplda5Ncsfw839P8tqqWmnzdbaou0929/3Dzz9O8lCS83e2Vazh2iy9kXR335PkxcOLGtPz2iT/dzLUMV3d/dUkTywrnnxfuCNLb9zLvS7Jke5+YnhjPpLkyqk1dAGsdC26+4vd/dTw6z1Z2jOXKVvlvtiIjXzu4jSsdS2Gz6u/keRTM23Uglrjc+xU3jMEu+c6P8kjE7+fyHODxM/rDG8eTyb5azNp3QKrpSGvr0zytRUO/72q+t9V9QdV9YqZNmyxdJIvVtV9VbV/heMbuX/YXtdl9Tdo98XsnNvdJ5OlN/IkL12hjvtj9v55kj9Y5dh6r2dsj5uGYbEHVxlu5r6Yrb+f5LHufniV4+6LKVn2OXYq7xmC3XOt1PO2fLzqRuqwjarqRUk+k+Sd3f2jZYfvT/LL3f0rSf5Tkv8x6/YtkNd096uSXJXkxmG4xyT3xgxV1fOT/HqS/7bCYffF/HF/zFBV/dskTyX55CpV1ns9Y+tuS/LyJJckOZnkgyvUcV/M1vVZu7fOfTEF63yOXfW0FcrWvDcEu+c6keTCid8vSPLoanWq6swkfzWbG37ABlTV87J0M3yyuz+7/Hh3/6i7fzL8fDjJ86rq7Bk3cyF096PD98eTfC5LQ2gmbeT+YftcleT+7n5s+QH3xcw99syw4+H74yvUcX/MyLDIwK8l+We9ymICG3g9Y4u6+7Hufrq7f5bkY1n539h9MSPDZ9Z/nOTTq9VxX2y/VT7HTuU9Q7B7rnuT7K2qi4b/G35dkkPL6hxK8szKNP8kS5O0/d+lKRjGgn88yUPd/aFV6vz1Z+Y4VtWlWfrv+s9m18rFUFUvHCb+pqpemOSKJN9cVu1QkjfXkldnaXL2yRk3dZGs+n9e3RczN/m+cEOSz69Q5+4kV1TVWcOQtCuGMrZRVV2Z5N1Jfr27f7pKnY28nrFFy+ZYvyEr/xtv5HMX2+MfJfl2d59Y6aD7Yvut8Tl2Ku8ZZ269ybvLsIrWTVn6hzsjycHufrCqbk1ytLsPZekC/deqOpalnrrrdq7Fu95rkrwpyTcmluZ9b5KXJUl3/5cshet/WVVPJfl/Sa4TtKfi3CSfG7LCmUl+v7v/sKrenvz8WhzO0oqYx5L8NMlbdqitu15V/ZUsrSL3tomyyWvhvpiSqvpUksuTnF1VJ5LckuRAkruq6q1JvpfkjUPdfUne3t3/orufqKr3Z+mDbJLc2t1Ge2zBKtfiPUlekOTI8Hp1z7CK9S8l+d3uvjqrvJ7twFPYNVa5FpdX1SVZGj52PMPr1eS1WO1z1w48hV1jpWvR3R/PCnOy3RdTt9rn2Km8Z9juAAAAYOQMxQQAABg5wQ4AAGDkBDsAAICRE+wAAABGTrADAAAYOcEOAABg5AQ7AACAkRPsAAAARu7/A4YW1otnD+yvAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 1080x360 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.hist(q2vals,bins=100);\n",
"fig = plt.gcf()\n",
"fig.set_size_inches(15,5)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"def momenta_map(q2_, cos_theta_l_, m_l):\n",
" n_events = q2_.shape[0]\n",
" \n",
" q2 = q2_.reshape(n_events,1)\n",
" cos_theta_l = cos_theta_l_.reshape(n_events,1)\n",
"\n",
" kin_f = (np.square(mB_mass)+np.square(mKst_mass)-q2)/(2*mB_mass)\n",
"\n",
" sin_theta_l = np.sqrt(1-np.square(cos_theta_l))\n",
" \n",
" lambd = lambda_phase_np(mB_mass,mKst_mass,q2)\n",
" \n",
" beta = beta_l_np(q2, m_l) \n",
" \n",
" phi_l = np.random.uniform(low=-pi,high=pi, size=(n_events,1))\n",
" phi_Kst = np.random.uniform(low=-pi,high=pi, size=(n_events,1))\n",
" theta_Kst = np.random.uniform(low=-pi,high=pi, size=(n_events,1))\n",
" \n",
" cos_phi_l = np.cos(phi_l)\n",
" sin_phi_l = np.sin(phi_l)\n",
" \n",
" cos_phi_Kst = np.cos(phi_Kst)\n",
" sin_phi_Kst = np.sin(phi_Kst)\n",
" \n",
" cos_theta_Kst = np.cos(theta_Kst)\n",
" sin_theta_Kst = np.sin(theta_Kst)\n",
"\n",
" p1_B_0 = 0.5*(mB_mass-kin_f) - (beta*np.sqrt(lambd)*cos_theta_l)/(4*mB_mass)\n",
" p1_B_1 = 0.5*beta*np.sqrt(q2)*sin_theta_l*cos_phi_l\n",
" p1_B_2 = 0.5*beta*np.sqrt(q2)*sin_theta_l*sin_phi_l\n",
" p1_B_3 = -np.sqrt(lambd)/(4*mB_mass) + 0.5*(mB_mass-kin_f)*beta*cos_theta_l\n",
"\n",
" p2_B_0 = 0.5*(mB_mass-kin_f) + (beta*np.sqrt(lambd)*cos_theta_l)/(4*mB_mass)\n",
" p2_B_1 = -0.5*beta*np.sqrt(q2)*sin_theta_l*cos_phi_l\n",
" p2_B_2 = -0.5*beta*np.sqrt(q2)*sin_theta_l*sin_phi_l\n",
" p2_B_3 = -np.sqrt(lambd)/(4*mB_mass) - 0.5*(mB_mass-kin_f)*beta*cos_theta_l\n",
" \n",
" pKst_B_0 = kin_f\n",
" pKst_B_1 = (np.sqrt(lambd)/(2*mB_mass))*(sin_theta_Kst*cos_phi_Kst)\n",
" pKst_B_2 = (np.sqrt(lambd)/(2*mB_mass))*(sin_theta_Kst*sin_phi_Kst)\n",
" pKst_B_3 = (np.sqrt(lambd)/(2*mB_mass))*cos_theta_Kst\n",
"\n",
" \n",
" pKst_B = np.concatenate([pKst_B_0,\n",
" pKst_B_1,\n",
" pKst_B_2,\n",
" pKst_B_3], axis =1)\n",
" \n",
" p1_B = np.concatenate([p1_B_0,\n",
" p1_B_1,\n",
" p1_B_2,\n",
" p1_B_3], axis =1)\n",
" \n",
" p2_B = np.concatenate([p2_B_0,\n",
" p2_B_1,\n",
" p2_B_2,\n",
" p2_B_3], axis =1)\n",
" \n",
" return pKst_B, p1_B, p2_B"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"pKst_B, p1_B, p2_B = momenta_map(q2vals_MC, cos_theta_vals_MC, mmu_mass)"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(array([0.011025, 0.011025, 0.011025, ..., 0.011025, 0.011025, 0.011025]),\n",
" array([0.795664, 0.795664, 0.795664, ..., 0.795664, 0.795664, 0.795664]))"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"scalar_product_4_np(p2_B,p2_B),scalar_product_4_np(pKst_B,pKst_B)"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(5.280000000000001, -0.00010350998435309136)"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"(pKst_B+p1_B+p2_B)[:,0].mean(), (pKst_B+p1_B+p2_B)[:,1].mean()\n"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA3MAAAEvCAYAAADvmpjfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAASp0lEQVR4nO3dX6jk51kH8O9j0nphC0lJGsNm4wZZxFg0LUsS6U01mm6yxbRCJLlo1hpYLxJooRdurRBpLayIFYs1utKlKdTGQFuymNV0DZUgmJptCW3SNHaJa3PcJZua0FYCSvT14vw2OUnOzpx/OzPvzOcDhzPzzu/MPAOzh/Pd5/1TrbUAAADQlx+bdgEAAACsnzAHAADQIWEOAACgQ8IcAABAh4Q5AACADglzAAAAHTp/2gWMctFFF7UdO3ZMuwwAAICp+PrXv/791trFqz0202Fux44dOXbs2LTLAAAAmIqq+vezPWaaJQAAQIeEOQAAgA4JcwAAAB0S5gAAADokzAEAAHRImAMAAOiQMAcAANAhYQ4AAKBDwhwAAECHhDkAAIAOCXMAAAAdOn/aBQCwGHbsf2Dk4ycO7JlQJZOziO8ZgMnRmQMAAOiQMAcAANAhYQ4AAKBDwhwAAECHhDkAAIAO2c0SADZh3I6V5+p57YQJgDAHwNwbFYyEIgB6JcwB0IVz1QEDgF5ZMwcAANAhnTkAFpq1aQD0SpgDgBFM7wRgVplmCQAA0CGdOQCYEl0/ADZDmANgywgnADA5whwAdMjGLQBYMwcAANAhYQ4AAKBDplkCwBwaNQ3TFEyA+SDMATATbJ4CAOtjmiUAAECHhDkAAIAOmWYJwJqZCgkAs0NnDgAAoEPCHAAAQIeEOQAAgA4JcwAAAB2yAQoAL7PBCQD0Q5gDAF5lVKg/cWDPBCsBYJSx0yyrantVfbWqnqyqJ6rqg8P4W6rqaFV9d/h+4TBeVfWpqjpeVd+sqneseK69w/Xfraq95+5tAQAAzLe1rJl7KcmHW2s/m+TaJHdU1ZVJ9id5qLW2M8lDw/0kuSHJzuFrX5K7k+Xwl+SuJNckuTrJXWcCIAAAAOszNsy11k611r4x3P5RkieTbEtyU5J7hsvuSfLe4fZNST7Xlj2S5IKqujTJu5Mcba0931p7IcnRJLu39N0AAAAsiHWtmauqHUnenuRrSS5prZ1KlgNfVb11uGxbkmdW/NjSMHa28de+xr4sd/Ry+eWXr6c8ANbAJicAMB/WHOaq6k1JvpjkQ621H1bVWS9dZayNGH/1QGsHkxxMkl27dr3ucQBgcwR6gPmwpnPmquoNWQ5yn2+tfWkYfnaYPpnh++lhfCnJ9hU/flmSkyPGAQAAWKe17GZZST6T5MnW2idXPHQ4yZkdKfcmuX/F+G3DrpbXJvnBMB3zwSTXV9WFw8Yn1w9jAAAArNNaplm+M8n7k3yrqh4bxn43yYEk91XV7Um+l+Tm4bEjSW5McjzJi0k+kCStteer6uNJHh2u+1hr7fkteRcAwESMm6LpHDqAyRkb5lpr/5TV17slyXWrXN+S3HGW5zqU5NB6CgQAAOD11rRmDgAAgNkizAEAAHRImAMAAOiQMAcAANChNR8aDsDsGLWjoN0EAWAx6MwBAAB0SGcOANgyzqEDmBydOQAAgA7pzAHMmXGdEQBgPujMAQAAdEiYAwAA6JAwBwAA0CFr5gCAmWAnTID10ZkDAADokDAHAADQIWEOAACgQ8IcAABAh4Q5AACADglzAAAAHRLmAAAAOiTMAQAAdEiYAwAA6JAwBwAA0CFhDgAAoEPCHAAAQIfOn3YBAMDi2LH/gWmXADA3hDmAKRj3B+2JA3smVAkA0CthDmAG6V4AAONYMwcAANAhYQ4AAKBDwhwAAECHhDkAAIAOCXMAAAAdEuYAAAA6JMwBAAB0SJgDAADokEPDAc4RB3/D5Iz793biwJ4JVQIwOTpzAAAAHRLmAAAAOiTMAQAAdMiaOQCgC5tZh2pNHTCPdOYAAAA6JMwBAAB0SJgDAADokDVzABvkHDkAYJp05gAAADokzAEAAHRobJirqkNVdbqqHl8x9vtV9R9V9djwdeOKxz5SVcer6qmqeveK8d3D2PGq2r/1bwUAAGBxrKUz99kku1cZ/5PW2lXD15Ekqaork9yS5OeGn/nzqjqvqs5L8ukkNyS5Msmtw7UAAABswNgNUFprD1fVjjU+301J7m2t/XeSf6uq40muHh473lp7Okmq6t7h2m+vu2IAAAA2tZvlnVV1W5JjST7cWnshybYkj6y4ZmkYS5JnXjN+zSZeGwBgy4zanfbEgT0TrARg7Ta6AcrdSX46yVVJTiX542G8Vrm2jRh/naraV1XHqurYc889t8HyAAAA5tuGOnOttWfP3K6qv0ryt8PdpSTbV1x6WZKTw+2zjb/2uQ8mOZgku3btWjXwAUyKs+QAgFm1oc5cVV264u77kpzZ6fJwkluq6ser6ookO5P8S5JHk+ysqiuq6o1Z3iTl8MbLBgAAWGxjO3NV9YUk70pyUVUtJbkrybuq6qosT5U8keS3k6S19kRV3ZfljU1eSnJHa+1/h+e5M8mDSc5Lcqi19sSWvxsAAIAFsZbdLG9dZfgzI67/RJJPrDJ+JMmRdVUHAADAqjazmyVA96yJAwB6tdHdLAEAAJgiYQ4AAKBDwhwAAECHhDkAAIAOCXMAAAAdEuYAAAA65GgCAIARxh1hcuLAnglVAvBqOnMAAAAdEuYAAAA6JMwBAAB0yJo5oHuj1rNYywIAzCudOQAAgA7pzAFzbdwudAAAvdKZAwAA6JAwBwAA0CFhDgAAoEPWzAEAbMK4tbl21QXOFZ05AACADglzAAAAHRLmAAAAOiTMAQAAdMgGKAAA59CoDVJsjgJshs4cAABAh4Q5AACADglzAAAAHRLmAAAAOmQDFGDmjdo8AKBn436/2SAFGEVnDgAAoEPCHAAAQIeEOQAAgA4JcwAAAB0S5gAAADokzAEAAHRImAMAAOiQMAcAANAhYQ4AAKBDwhwAAECHhDkAAIAOCXMAAAAdOn/aBQAAsLod+x8462MnDuyZYCXALNKZAwAA6JAwBwAA0CHTLIGZMGoqEQAAr6czBwAA0CGdOWAidN4AALaWzhwAAECHhDkAAIAOjZ1mWVWHkrwnyenW2tuGsbck+ZskO5KcSPIbrbUXqqqS/GmSG5O8mOQ3W2vfGH5mb5LfG572D1pr92ztWwEAWBzjpq87hw7m31o6c59Nsvs1Y/uTPNRa25nkoeF+ktyQZOfwtS/J3cnL4e+uJNckuTrJXVV14WaLBwAAWFRjw1xr7eEkz79m+KYkZzpr9yR574rxz7VljyS5oKouTfLuJEdba8+31l5IcjSvD4gAAACs0UbXzF3SWjuVJMP3tw7j25I8s+K6pWHsbOMAAABswFZvgFKrjLUR469/gqp9VXWsqo4999xzW1ocAADAvNjoOXPPVtWlrbVTwzTK08P4UpLtK667LMnJYfxdrxn/x9WeuLV2MMnBJNm1a9eqgQ8AgNFskALzb6OducNJ9g639ya5f8X4bbXs2iQ/GKZhPpjk+qq6cNj45PphDAAAgA1Yy9EEX8hyV+2iqlrK8q6UB5LcV1W3J/lekpuHy49k+ViC41k+muADSdJae76qPp7k0eG6j7XWXrupCtCxcf8DDADA1hob5lprt57loetWubYlueMsz3MoyaF1VQcAAMCqtnoDFAAAACZgoxugAAvIVEoAgNmhMwcAANAhYQ4AAKBDplkCACygUVPnnUEHfdCZAwAA6JAwBwAA0CFhDgAAoEPCHAAAQIdsgAK8zDlyAAD90JkDAADokDAHAADQIWEOAACgQ8IcAABAh2yAAgDAq4zbEOvEgT0TqgQYRWcOAACgQzpzsGAcPwAAMB905gAAADqkMwcAwLqMmuVhPR1Mjs4cAABAh4Q5AACADglzAAAAHRLmAAAAOiTMAQAAdEiYAwAA6JAwBwAA0CHnzMGcGXX2DwAA80NnDgAAoEPCHAAAQIdMswQAYMuMm+5/4sCeCVUC809nDgAAoEPCHAAAQIeEOQAAgA4JcwAAAB0S5gAAADpkN0sAACbGbpewdXTmAAAAOiTMAQAAdEiYAwAA6JAwBwAA0CEboEBnxi0cBwBgMQhzMIMENgAAxhHmAACYGaP+Q9OxBfBq1swBAAB0SJgDAADokGmWMAXWxAEAsFk6cwAAAB3aVJirqhNV9a2qeqyqjg1jb6mqo1X13eH7hcN4VdWnqup4VX2zqt6xFW8AAABgEW1FZ+6XWmtXtdZ2Dff3J3motbYzyUPD/SS5IcnO4Wtfkru34LUBAAAW0rmYZnlTknuG2/ckee+K8c+1ZY8kuaCqLj0Hrw8AADD3NrsBSkvylapqSf6ytXYwySWttVNJ0lo7VVVvHa7dluSZFT+7NIyd2mQNAAAsgHEbiDmHjkWz2TD3ztbaySGwHa2q74y4tlYZa6+7qGpflqdh5vLLL99keQAAAPNpU9MsW2snh++nk3w5ydVJnj0zfXL4fnq4fCnJ9hU/flmSk6s858HW2q7W2q6LL754M+UBAADMrQ2Huar6iap685nbSa5P8niSw0n2DpftTXL/cPtwktuGXS2vTfKDM9MxAQAAWJ/NTLO8JMmXq+rM8/x1a+3vq+rRJPdV1e1Jvpfk5uH6I0luTHI8yYtJPrCJ14aZ5lBwAADOtQ2Hudba00l+YZXx/0xy3SrjLckdG309AAAAXrHZDVAAAGAmjJoZY6dL5tG5OGcOAACAc0yYAwAA6JAwBwAA0CFr5uAsxu1Iae49AADTJMzBBjl+AACAaTLNEgAAoEM6cwAAzD3LJ5hHwhwAAAtP2KNHplkCAAB0SJgDAADokDAHAADQIWEOAACgQzZAYaE5Kw4AgF4JcwAAMMao/wC20yXTYpolAABAh4Q5AACADplmCQAAm+DAcaZFZw4AAKBDOnPMNbtVAgAwr3TmAAAAOiTMAQAAdEiYAwAA6JAwBwAA0CEboNA1G5wAALCodOYAAAA6pDMHAADn0KiZRA4UZzOEOWaeqZQAwLwa93eOsMcoplkCAAB0SJgDAADokGmWTJ1plAAAsH46cwAAAB3SmQMAgE7ZKXOxCXMAADCjLEdhFGGOifCLCAAAtpY1cwAAAB3SmQMAgDnkQPL5J8yxJUyjBACAyRLmAABgAenc9c+aOQAAgA7pzLEmplECAMBs0ZkDAADokM4cL9N9AwDgjFF/G1pPNxt05gAAADqkM7dAdN4AAGB+6MwBAAB0SGeuM7prAABMmzPqZoMwBwAAbKnNNCAEwbWb+DTLqtpdVU9V1fGq2j/p1wcAAJgHE+3MVdV5ST6d5FeTLCV5tKoOt9a+Pck6ps1USQAAWJ0pnGs36WmWVyc53lp7Okmq6t4kNyVZqDAHAABsjPPvXjHpMLctyTMr7i8luWbCNWyazhoAAMyeRVurN+kwV6uMtVddULUvyb7h7n9V1VPnvKrZc1GS70+7CNgkn2Pmgc8x88DnmHlwzj/H9Yfn8tk35afO9sCkw9xSku0r7l+W5OTKC1prB5McnGRRs6aqjrXWdk27DtgMn2Pmgc8x88DnmHngc7y6Se9m+WiSnVV1RVW9McktSQ5PuAYAAIDuTbQz11p7qaruTPJgkvOSHGqtPTHJGgAAAObBxA8Nb60dSXJk0q/bmYWeZsrc8DlmHvgcMw98jpkHPserqNba+KsAAACYKZNeMwcAAMAWEOZmVFX9UVV9p6q+WVVfrqoLpl0TrFdV3VxVT1TV/1WVHajoSlXtrqqnqup4Ve2fdj2wXlV1qKpOV9Xj064FNqqqtlfVV6vqyeFvig9Ou6ZZIszNrqNJ3tZa+/kk/5rkI1OuBzbi8SS/nuThaRcC61FV5yX5dJIbklyZ5NaqunK6VcG6fTbJ7mkXAZv0UpIPt9Z+Nsm1Se7w+/gVwtyMaq19pbX20nD3kSyfyQddaa092Vp7atp1wAZcneR4a+3p1tr/JLk3yU1TrgnWpbX2cJLnp10HbEZr7VRr7RvD7R8leTLJtulWNTuEuT78VpK/m3YRAAtkW5JnVtxfij8eAKaqqnYkeXuSr023ktkx8aMJeEVV/UOSn1zloY+21u4frvloltvLn59kbbBWa/kcQ4dqlTHbPwNMSVW9KckXk3yotfbDadczK4S5KWqt/cqox6tqb5L3JLmuOUOCGTXucwydWkqyfcX9y5KcnFItAAutqt6Q5SD3+dbal6ZdzywxzXJGVdXuJL+T5Ndaay9Oux6ABfNokp1VdUVVvTHJLUkOT7kmgIVTVZXkM0mebK19ctr1zBphbnb9WZI3JzlaVY9V1V9MuyBYr6p6X1UtJfnFJA9U1YPTrgnWYtiA6s4kD2Z5sf19rbUnplsVrE9VfSHJPyf5mapaqqrbp10TbMA7k7w/yS8PfxM/VlU3TruoWVFm7wEAAPRHZw4AAKBDwhwAAECHhDkAAIAOCXMAAAAdEuYAAAA6JMwBAAB0SJgDAADokDAHAADQof8HWgq4ecLqlakAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 1080x360 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.hist(p1_B[:,2],bins=100);\n",
"fig = plt.gcf()\n",
"fig.set_size_inches(15,5)"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA3MAAAEvCAYAAADvmpjfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAXvElEQVR4nO3df6yeZ3kf8O9FEmg10JIshqWOmaM13RrQGjovyYQ0MWgTk1QLTKMKf0CgkdJNyQQSmjDtNFpoJHdrQVRlTK7ikUyMEBUQFqRLDQUhpAFxsjTFhAwveMRNRMwcfgktU9Jrf5wn9CQ5Pj98js977nM+H+nVed/ruZ/3vR/ptX2+vn881d0BAABgLM+bdQcAAABYOWEOAABgQMIcAADAgIQ5AACAAQlzAAAAAxLmAAAABnTmrDuwmPPOO6937tw5624AAADMxD333PPd7t620LENHeZ27tyZQ4cOzbobAAAAM1FV//tkx0yzBAAAGJAwBwAAMCBhDgAAYEDCHAAAwICEOQAAgAEJcwAAAAMS5gAAAAYkzAEAAAxoyTBXVT9VVV+tqj+vqsNV9dtT/cNV9a2qum96XDLVq6r+oKqOVNX9VfWL897ruqr65vS47vRdFgAAwOZ25jLaPJHk1d39o6o6K8mXqupPpmP/prv/+FntX5vkoulxWZIPJbmsqs5N8u4ku5J0knuq6kB3P74WFwIAALCVLDky13N+NL08a3r0Iqdck+S26bwvJzm7qs5PcmWSg919YgpwB5PsXl33AQAAtqbljMylqs5Ick+Sn03ywe7+SlX9qyQ3V9W/S/K5JHu6+4kk25M8PO/0Y1PtZHUAmJmdez6z6PGje69ep54AwMosawOU7n6quy9JckGSS6vq5UneleTvJ/lHSc5N8s6peS30FovUn6GqbqiqQ1V16Pjx48vpHgAAwJazot0su/t7Sb6QZHd3PzpNpXwiyX9OcunU7FiSHfNOuyDJI4vUn/0Z+7p7V3fv2rZt20q6BwAAsGUsZzfLbVV19vT8p5P8UpJvTOvgUlWV5HVJvjadciDJm6ddLS9P8v3ufjTJXUmuqKpzquqcJFdMNQAAAFZoOWvmzk9y67Ru7nlJ7ujuT1fVn1XVtsxNn7wvyb+c2t+Z5KokR5L8OMlbk6S7T1TVe5PcPbV7T3efWLtLAQAA2DqWDHPdfX+SVyxQf/VJ2neSG09ybH+S/SvsIwAAAM+yojVzAAAAbAzCHAAAwICEOQAAgAEJcwAAAAMS5gAAAAYkzAEAAAxImAMAABiQMAcAADAgYQ4AAGBAwhwAAMCAzpx1BwBgI9u55zOLHj+69+p16gkAPJOROQAAgAEJcwAAAAMS5gAAAAYkzAEAAAxImAMAABiQMAcAADAgYQ4AAGBAwhwAAMCAhDkAAIABCXMAAAADEuYAAAAGJMwBAAAMSJgDAAAYkDAHAAAwIGEOAABgQMIcAADAgIQ5AACAAS0Z5qrqp6rqq1X151V1uKp+e6pfWFVfqapvVtXHqur5U/0F0+sj0/Gd897rXVP9waq68nRdFAAAwGa3nJG5J5K8urt/IcklSXZX1eVJfjfJ+7v7oiSPJ7l+an99kse7+2eTvH9ql6q6OMm1SV6WZHeS/1hVZ6zlxQAAAGwVS4a5nvOj6eVZ06OTvDrJH0/1W5O8bnp+zfQ60/HXVFVN9du7+4nu/laSI0kuXZOrAAAA2GKWtWauqs6oqvuSPJbkYJL/leR73f3k1ORYku3T8+1JHk6S6fj3k/yt+fUFzgEAAGAFlhXmuvup7r4kyQWZG037+YWaTT/rJMdOVn+Gqrqhqg5V1aHjx48vp3sAAABbzop2s+zu7yX5QpLLk5xdVWdOhy5I8sj0/FiSHUkyHf+bSU7Mry9wzvzP2Nfdu7p717Zt21bSPQAAgC1jObtZbquqs6fnP53kl5I8kOTzSf7F1Oy6JJ+anh+YXmc6/mfd3VP92mm3ywuTXJTkq2t1IQAAAFvJmUs3yflJbp12nnxekju6+9NV9fUkt1fV7yT5H0lumdrfkuS/VNWRzI3IXZsk3X24qu5I8vUkTya5sbufWtvLAQAA2BqWDHPdfX+SVyxQfygL7EbZ3f83yRtO8l43J7l55d0EAABgvhWtmQMAAGBjEOYAAAAGJMwBAAAMSJgDAAAYkDAHAAAwIGEOAABgQMIcAADAgIQ5AACAAQlzAAAAAxLmAAAABiTMAQAADEiYAwAAGJAwBwAAMCBhDgAAYEDCHAAAwICEOQAAgAEJcwAAAAMS5gAAAAYkzAEAAAxImAMAABiQMAcAADCgM2fdAQA43Xbu+cysuwAAa87IHAAAwICEOQAAgAGZZgkAq7DYFM6je69ex54AsNUYmQMAABiQMAcAADAgYQ4AAGBAS4a5qtpRVZ+vqgeq6nBVvW2q/1ZV/WVV3Tc9rpp3zruq6khVPVhVV86r755qR6pqz+m5JAAAgM1vORugPJnkHd19b1W9KMk9VXVwOvb+7v69+Y2r6uIk1yZ5WZKfSfLZqvq56fAHk/xykmNJ7q6qA9399bW4EAAAgK1kyTDX3Y8meXR6/sOqeiDJ9kVOuSbJ7d39RJJvVdWRJJdOx45090NJUlW3T22FOQAAgBVa0Zq5qtqZ5BVJvjKVbqqq+6tqf1WdM9W2J3l43mnHptrJ6gAAAKzQssNcVb0wyceTvL27f5DkQ0n+bpJLMjdy9/tPN13g9F6k/uzPuaGqDlXVoePHjy+3ewAAAFvKssJcVZ2VuSD3ke7+RJJ093e6+6nu/qskf5S/nkp5LMmOeadfkOSRRerP0N37untXd+/atm3bSq8HAABgS1jObpaV5JYkD3T3++bVz5/X7PVJvjY9P5Dk2qp6QVVdmOSiJF9NcneSi6rqwqp6fuY2STmwNpcBAACwtSxnN8tXJnlTkr+oqvum2m8keWNVXZK5qZJHk/x6knT34aq6I3MbmzyZ5MbufipJquqmJHclOSPJ/u4+vIbXAgAAsGUsZzfLL2Xh9W53LnLOzUluXqB+52LnAQAAsDwr2s0SAACAjUGYAwAAGJAwBwAAMCBhDgAAYEDCHAAAwICEOQAAgAEJcwAAAAMS5gAAAAYkzAEAAAzozFl3AAA2q517PrPo8aN7r16nngCwGRmZAwAAGJAwBwAAMCBhDgAAYEDCHAAAwICEOQAAgAHZzRKA4S21ayQAbEZG5gAAAAYkzAEAAAxImAMAABiQMAcAADAgYQ4AAGBAdrMEYAh2rASAZzIyBwAAMCBhDgAAYEDCHAAAwICEOQAAgAEJcwAAAAMS5gAAAAa0ZJirqh1V9fmqeqCqDlfV26b6uVV1sKq+Of08Z6pXVf1BVR2pqvur6hfnvdd1U/tvVtV1p++yAAAANrfl3GfuySTv6O57q+pFSe6pqoNJ3pLkc929t6r2JNmT5J1JXpvkoulxWZIPJbmsqs5N8u4ku5L09D4Huvvxtb4oABjBYvfOO7r36nXsCQAjWnJkrrsf7e57p+c/TPJAku1Jrkly69Ts1iSvm55fk+S2nvPlJGdX1flJrkxysLtPTAHuYJLda3o1AAAAW8RyRuZ+oqp2JnlFkq8keUl3P5rMBb6qevHUbHuSh+eddmyqnaz+7M+4IckNSfLSl750Jd0DgE1jsVG7xMgdACsIc1X1wiQfT/L27v5BVZ206QK1XqT+zEL3viT7kmTXrl3POQ7A5rRUeAEAnmlZu1lW1VmZC3If6e5PTOXvTNMnM/18bKofS7Jj3ukXJHlkkToAAAArtJzdLCvJLUke6O73zTt0IMnTO1Jel+RT8+pvnna1vDzJ96fpmHcluaKqzpl2vrxiqgEAALBCy5lm+cokb0ryF1V131T7jSR7k9xRVdcn+XaSN0zH7kxyVZIjSX6c5K1J0t0nquq9Se6e2r2nu0+syVUAAABsMUuGue7+UhZe75Ykr1mgfSe58STvtT/J/pV0EAAAgOda1po5AAAANhZhDgAAYEDCHAAAwIBWdNNwADhV7iMHAGvLyBwAAMCAhDkAAIABCXMAAAADEuYAAAAGZAMUABjQUhvKHN179Tr1BIBZMTIHAAAwIGEOAABgQMIcAADAgKyZA4BNaLE1ddbTAWwORuYAAAAGZGQOgDWz1A6LAMDaMTIHAAAwIGEOAABgQMIcAADAgIQ5AACAAdkABYBls8EJAGwcwhwAbDFLhXL3oQMYg2mWAAAAAxLmAAAABmSaJQA/YU0cAIzDyBwAAMCAhDkAAIABmWYJADzDYtNt7XQJsHEYmQMAABjQkmGuqvZX1WNV9bV5td+qqr+sqvumx1Xzjr2rqo5U1YNVdeW8+u6pdqSq9qz9pQAAAGwdy5lm+eEkf5jktmfV39/dvze/UFUXJ7k2ycuS/EySz1bVz02HP5jkl5McS3J3VR3o7q+vou8AnAI7VgLA5rBkmOvuL1bVzmW+3zVJbu/uJ5J8q6qOJLl0Onakux9Kkqq6fWorzAEAAJyC1ayZu6mq7p+mYZ4z1bYneXhem2NT7WT156iqG6rqUFUdOn78+Cq6BwAAsHmd6m6WH0ry3iQ9/fz9JL+WpBZo21k4NPZCb9zd+5LsS5Jdu3Yt2AYAmI3VTtO1GybA2jmlMNfd33n6eVX9UZJPTy+PJdkxr+kFSR6Znp+sDgAAwAqd0jTLqjp/3svXJ3l6p8sDSa6tqhdU1YVJLkry1SR3J7moqi6squdnbpOUA6febQAAgK1tyZG5qvpoklclOa+qjiV5d5JXVdUlmZsqeTTJrydJdx+uqjsyt7HJk0lu7O6npve5KcldSc5Isr+7D6/51QAAAGwRy9nN8o0LlG9ZpP3NSW5eoH5nkjtX1DsAAAAWdKoboAAArNhiG6jYHAVgZVZzawIAAABmxMgcwIBWuz08ADA+I3MAAAADEuYAAAAGZJolwAZkGiUAsBQjcwAAAAMS5gAAAAYkzAEAAAzImjkAYENYzVpRNxwHtiJhDmAGbHACAKyWMAdwipYKZEYKYP348whsRcIcwGli9A0AOJ1sgAIAADAgYQ4AAGBAplkCLMJUSQBgozIyBwAAMCAjcwDApme3S2AzEuYAgC1vsbAn6AEblWmWAAAAAxLmAAAABmSaJQDAIqy3AzYqI3MAAAADEuYAAAAGZJolsKW5KTgAMCphDhieQAYAbEXCHADAKtggBZiVJcNcVe1P8itJHuvul0+1c5N8LMnOJEeT/Gp3P15VleQDSa5K8uMkb+nue6dzrkvyb6e3/Z3uvnVtLwUAYONxQ3LgdKnuXrxB1T9J8qMkt80Lc/8+yYnu3ltVe5Kc093vrKqrkvzrzIW5y5J8oLsvm8LfoSS7knSSe5L8w+5+fLHP3rVrVx86dGh1VwgMzzRKYKtaLOwZEYStoaru6e5dCx1bcjfL7v5ikhPPKl+T5OmRtVuTvG5e/bae8+UkZ1fV+UmuTHKwu09MAe5gkt0rvxQAAACSU18z95LufjRJuvvRqnrxVN+e5OF57Y5NtZPVgS3C/yADAKyttb7PXC1Q60Xqz32Dqhuq6lBVHTp+/Piadg4AAGCzONWRue9U1fnTqNz5SR6b6seS7JjX7oIkj0z1Vz2r/oWF3ri79yXZl8ytmTvF/gGDsS4O4Ln83Qgs5lRH5g4kuW56fl2ST82rv7nmXJ7k+9N0zLuSXFFV51TVOUmumGoAAACcguXcmuCjmRtVO6+qjiV5d5K9Se6oquuTfDvJG6bmd2ZuJ8sjmbs1wVuTpLtPVNV7k9w9tXtPdz97UxVgYP73GABgfS0Z5rr7jSc59JoF2naSG0/yPvuT7F9R7wAAAFjQWm+AAgAAwDo41Q1QAAAY2GLT490uBsYgzAHLZl0cwDj8nQ2bnzAHAMAzLBUEjdzBxmDNHAAAwICMzAE/YUoOAMA4hDkAAFbkdG6eYmMWWD7TLAEAAAZkZA62GFMpAQA2ByNzAAAAAzIyB4MxsgbAyPw7BmvHyBwAAMCAjMwBALBmjLzB+hHmYAPyDyEAAEsxzRIAAGBARuZgBoy8AQCwWsIcnAbCGgCsvaX+fT269+p16glsDKZZAgAADMjIHAAAm8JiI3dG7diMjMwBAAAMyMgcnCLr4gAAmCVhDk5CWAMAYCMT5gAA2PTshMlmJMwBALDlCXuMyAYoAAAAAxLmAAAABmSaJVuaTU4AgOVYze8MpmhyuhiZAwAAGNCqRuaq6miSHyZ5KsmT3b2rqs5N8rEkO5McTfKr3f14VVWSDyS5KsmPk7ylu+9dzeeDxcoAAGxVazHN8p9293fnvd6T5HPdvbeq9kyv35nktUkumh6XJfnQ9BNOG9MoAQDYrE7HNMtrktw6Pb81yevm1W/rOV9OcnZVnX8aPh8AAGDTW22Y6yR/WlX3VNUNU+0l3f1okkw/XzzVtyd5eN65x6YaAAAAK7TaaZav7O5HqurFSQ5W1TcWaVsL1Po5jeZC4Q1J8tKXvnSV3QMAgNlabNmH9f2sxqrCXHc/Mv18rKo+meTSJN+pqvO7+9FpGuVjU/NjSXbMO/2CJI8s8J77kuxLkl27dj0n7AEAwGZhMzdW45TDXFX9jSTP6+4fTs+vSPKeJAeSXJdk7/TzU9MpB5LcVFW3Z27jk+8/PR0TFmMTEwAAeK7VjMy9JMkn5+44kDOT/Nfu/m9VdXeSO6rq+iTfTvKGqf2dmbstwZHM3Zrgrav4bAAAgC3tlMNcdz+U5BcWqP+fJK9ZoN5JbjzVz2PzMvIGAAArtxb3mQMAAE6D1fynt/V2m58wBwAAm5DNVTY/YY41YaokAACsr9XeNBwAAIAZEOYAAAAGJMwBAAAMSJgDAAAYkDAHAAAwILtZAgDAFuTWBeMT5vgJtxcAAOBpi/1uKOhtDMLcFiKsAQDA5mHNHAAAwICEOQAAgAGZZgkAAKyIzVM2BmFuk7EuDgCAWbN5yvowzRIAAGBAwhwAAMCAhDkAAIABWTMHAACsG5unrB1hbjA2OAEAABLTLAEAAIZkZA4AANgwTMNcPmFugzGNEgAATs497P6aMDcDAhsAALBa1swBAAAMSJgDAAAYkGmWp4FplAAAsP622uYpRuYAAAAGtO5hrqp2V9WDVXWkqvas9+cDAABsBusa5qrqjCQfTPLaJBcneWNVXbyefQAAANgM1nvN3KVJjnT3Q0lSVbcnuSbJ19e5H6tiTRwAAIxns92jbr2nWW5P8vC818emGgAAACuw3iNztUCtn9Gg6oYkN0wvf1RVD572XjEr5yX57qw7wabje8Va851irflOsdZ8p9ZA/e6se3BSf+dkB9Y7zB1LsmPe6wuSPDK/QXfvS7JvPTvFbFTVoe7eNet+sLn4XrHWfKdYa75TrDXfqa1rvadZ3p3koqq6sKqen+TaJAfWuQ8AAADDW9eRue5+sqpuSnJXkjOS7O/uw+vZBwAAgM1gvadZprvvTHLnen8uG5LptJwOvlesNd8p1prvFGvNd2qLqu5euhUAAAAbynqvmQMAAGANCHPMVFX9h6r6RlXdX1WfrKqzZ90nxlZVb6iqw1X1V1VlZy9OWVXtrqoHq+pIVe2ZdX8YX1Xtr6rHquprs+4L46uqHVX1+ap6YPp3722z7hPrT5hj1g4meXl3/4Mk/zPJu2bcH8b3tST/PMkXZ90RxlVVZyT5YJLXJrk4yRur6uLZ9opN4MNJds+6E2waTyZ5R3f/fJLLk9zo76mtR5hjprr7T7v7yenllzN370E4Zd39QHc/OOt+MLxLkxzp7oe6+/8luT3JNTPuE4Pr7i8mOTHrfrA5dPej3X3v9PyHSR5Isn22vWK9CXNsJL+W5E9m3QmAzP1C9PC818filyRgg6qqnUlekeQrs+0J623db03A1lNVn03ytxc49Jvd/ampzW9mbrrAR9azb4xpOd8pWKVaoGb7Z2DDqaoXJvl4krd39w9m3R/WlzDHadfdv7TY8aq6LsmvJHlNu1cGy7DUdwrWwLEkO+a9viDJIzPqC8CCquqszAW5j3T3J2bdH9afaZbMVFXtTvLOJP+su3886/4ATO5OclFVXVhVz09ybZIDM+4TwE9UVSW5JckD3f2+WfeH2RDmmLU/TPKiJAer6r6q+k+z7hBjq6rXV9WxJP84yWeq6q5Z94nxTBsz3ZTkrsxtKnBHdx+eba8YXVV9NMl/T/L3qupYVV0/6z4xtFcmeVOSV0+/Q91XVVfNulOsrzKrDQAAYDxG5gAAAAYkzAEAAAxImAMAABiQMAcAADAgYQ4AAGBAwhwAAMCAhDkAAIABCXMAAAAD+v8Qha88bHPqsAAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 1080x360 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.hist(pKst_B[:,2],bins=100);\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": []
}
],
"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