{
"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\n",
"import phasespace\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def beta_l(q2, m):\n",
" return ztf.sqrt(1-(4*ztf.square(m))/(q2))"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"def spatial_sp(p1, p2):\n",
" \n",
" val=p1[0]*p2[0]+p1[1]*p2[1]+p1[2]*p2[2]\n",
" \n",
" return val"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def sp(p1, p2):\n",
" \n",
" val0 = p1[3]*p2[3]\n",
" \n",
" val = spatial_sp(p1, p2)\n",
" \n",
" sp = val0 - val\n",
" \n",
" return sp"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"def get_costheta_l(p1, p2): \n",
" \n",
" num = spatial_sp(p1, p2)\n",
" \n",
" den1= ztf.sqrt(spatial_sp(p1, p1))\n",
" den2= ztf.sqrt(spatial_sp(p2, p2))\n",
" \n",
" costheta_l = num/(den1*den2)\n",
" \n",
" return costheta_l"
]
},
{
"cell_type": "code",
"execution_count": 6,
"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": [
"from phasespace import Particle\n",
"mZ = 80000.00\n",
"mmu = 105.3\n",
"\n",
"\n",
"el1 = Particle('l1', mmu)\n",
"el2 = Particle('l2', mmu)\n",
"Z = Particle('Z', mZ).set_children(el1, el2)\n",
"\n",
"#weights, particles = B.generate(n_events=10000)\n",
"weights_np, particles_np = Z.generate(n_events=10000)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAD8CAYAAABU4IIeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAGulJREFUeJzt3X9sXed93/H3p2Qkt8li2TTbOZI8MhBbjEKyLLlRGqz1DHmx5Li1HExe6QSovBrQtlj7AS+rpXlLFscezHSLgiJ2YgF2qzrxaFWtIS52orqzXQxFJonyz1AOI0ZWLUZBLU2Si6yIHNrf/XEexVfX98fh5UPySvq8AELnPud5vs/3Hh7yy/PjHikiMDMzy+XnFjoBMzM7v7iwmJlZVi4sZmaWlQuLmZll5cJiZmZZubCYmVlWLixmZpaVC4uZmWXlwmJmZll1L3QCC+Gyyy6Lvr6+hU7DzOycsn///uMR0duq3wVZWPr6+hgbG1voNMzMzimS/qpMP58KMzOzrFxYzMwsKxcWMzPLyoXFzMyycmExM7OsXFjMzCwrFxYzM8vKhcXMzLIqVVgkrZU0IWlS0uY66xdLeiSt3yOpr2rdltQ+IWlNq5iS+lOMgynmotR+paRnJE1LWl8nh3dL+qGkr8xsE5iZWU4tP3kvqQu4F/gYMAXskzQaEQequt0CnIyIFZKGgGHgtyQNAkPASuA9wJ9L+uU0plHMYWBrRIxI+lqK/VXgFeBm4DMNUv0C8Bfl37pd6Po2P9Zw3eF7rpvHTMzOL2WOWFYBkxFxKCJeB0aAdTV91gHb0/JO4GpJSu0jEXE6Il4GJlO8ujHTmNUpBinmDQARcTgiXgDerE1Q0oeAXwL+rOT7NjOzOVKmsCwFjlS9nkptdftExDTwGtDTZGyj9h7gVIrRaK6zSPo54L8D/6HEezEzszlWprCoTluU7JOrvZlPA49HxJFmnSRtlDQmaezYsWMtQpqZWbvKPN14Clhe9XoZcLRBnylJ3cDFwIkWY+u1HweWSOpORy315qr1UeDXJX0aeBewSNKPI+KsmwwiYhuwDaBSqbQqVmZm1qYyhWUfMCCpH/ghxcX4T9b0GQU2AN8B1gNPRkRIGgUelvQliov3A8BeiiOTt8VMY55KMUZSzF3NkouIT51ZlnQzUKktKma2sJrdKAG+WeJ807KwRMS0pE3AbqALeDAixiXdCYxFxCjwAPCQpEmKI5WhNHZc0g7gADAN3BoRbwDUi5mmvB0YkXQX8GyKjaQPA48ClwC/KenzEbEyy1Yw/+CbWTal/qOviHgceLym7bNVyz8Bbmww9m7g7jIxU/shirvGatv3UZwaa5bnHwJ/2KyPmZnNrQvyf5A0s3OHj6bPPX6ki5mZZeXCYmZmWflUmFkdsz390qmPi+nUvOz84sJi561WxcHM5oZPhZmZWVYuLGZmlpULi5mZZeXCYmZmWbmwmJlZVi4sZmaWlW83Po906mcUOjWvc5Efb/J23iadx4Uls9l8dqKTfwBcHMysLBcWs/OIPxRqncCFxczOaz7ann8uLHZO81/oZp3Hd4WZmVlWPmIxm2e+i+ntOvXI06fR2uPCYmZZdGpxsPnnwnIO8Q/uhWGhvs8+kuoc5/r3otQ1FklrJU1ImpS0uc76xZIeSev3SOqrWrcltU9IWtMqpqT+FONgirkotV8p6RlJ05LWV/X/gKTvSBqX9IKk32pvU5iZWQ4tC4ukLuBe4FpgELhJ0mBNt1uAkxGxAtgKDKexg8AQsBJYC9wnqatFzGFga0QMACdTbIBXgJuBh2vm/lvgtyPizBxflrSk3Ns3M7PcypwKWwVMRsQhAEkjwDrgQFWfdcB/Scs7ga9IUmofiYjTwMuSJlM86sWU9BKwGvhk6rM9xf1qRBxOfd+sTi4ivl+1fFTSq0AvcKrEe7MO16mn/zo1L5s/nXq6qhPyKnMqbClwpOr1VGqr2ycipoHXgJ4mYxu19wCnUoxGczUkaRWwCPhB2TFmZpZXmSMW1WmLkn0atdcraM36tyTpcuAhYENEvFln/UZgI8AVV1xRJuS881/BZnY+KFNYpoDlVa+XAUcb9JmS1A1cDJxoMbZe+3FgiaTudNRSb663kfRu4DHgP0XE/6nXJyK2AdsAKpVKqWJlc8/F1Oz8U6aw7AMGJPUDP6S4GP/Jmj6jwAbgO8B64MmICEmjwMOSvgS8BxgA9lIcmbwtZhrzVIoxkmLuapZcumvsUeCPIuKPS7wfM7NzWqf/QdbyGks6ctgE7AZeAnZExLikOyVdn7o9APSki/O3AZvT2HFgB8WF/m8Dt0bEG41ipli3A7elWD0pNpI+LGkKuBG4X9KZ/v8MuBK4WdJz6esDs9gmZmY2C6U+IBkRjwOP17R9tmr5JxS/8OuNvRu4u0zM1H6It+4cq27fR3FqrLb968DXW74JMzObF34IpZmZZeVHutisdfr5XjObXy4sFwj/8rfZ8j5kZbmwtME/YGZmjfkai5mZZeXCYmZmWflUmJnZHOiEh0EuFB+xmJlZVj5iMbML1kLeiHM+3wTkIxYzM8vKhcXMzLJyYTEzs6xcWMzMLCsXFjMzy8qFxczMsnJhMTOzrFxYzMwsKxcWMzPLyoXFzMyycmExM7OsShUWSWslTUialLS5zvrFkh5J6/dI6qtatyW1T0ha0yqmpP4U42CKuSi1XynpGUnTktbXzL8h9T8oacPMN4OZmeXSsrBI6gLuBa4FBoGbJA3WdLsFOBkRK4CtwHAaOwgMASuBtcB9krpaxBwGtkbEAHAyxQZ4BbgZeLgmv0uBzwEfAVYBn5N0SdkNYGZmeZU5YlkFTEbEoYh4HRgB1tX0WQdsT8s7gaslKbWPRMTpiHgZmEzx6sZMY1anGKSYNwBExOGIeAF4s2buNcATEXEiIk4CT1AUMTMzWwBlCstS4EjV66nUVrdPREwDrwE9TcY2au8BTqUYjeZqJz8zM5snZQqL6rRFyT652pspNUbSRkljksaOHTvWIqSZmbWrTGGZApZXvV4GHG3UR1I3cDFwosnYRu3HgSUpRqO52smPiNgWEZWIqPT29rYIaWZm7SpTWPYBA+lurUUUF+NHa/qMAmfuxloPPBkRkdqH0l1j/cAAsLdRzDTmqRSDFHNXi/x2A9dIuiRdtL8mtZmZ2QJoWVjS9Y5NFL+sXwJ2RMS4pDslXZ+6PQD0SJoEbgM2p7HjwA7gAPBt4NaIeKNRzBTrduC2FKsnxUbShyVNATcC90saT3OcAL5AUaz2AXemNjMzWwAqDhIuLJVKJcbGxtoefz7/X9Vmdn47fM91bY+VtD8iKq36+ZP3ZmaWlQuLmZll5cJiZmZZubCYmVlWLixmZpaVC4uZmWXlwmJmZlm5sJiZWVYuLGZmlpULi5mZZeXCYmZmWbmwmJlZVi4sZmaWlQuLmZll5cJiZmZZubCYmVlWLixmZpaVC4uZmWXlwmJmZlm5sJiZWVYuLGZmllWpwiJpraQJSZOSNtdZv1jSI2n9Hkl9Veu2pPYJSWtaxZTUn2IcTDEXNZtD0jskbZf0oqSXJG1pd2OYmdnstSwskrqAe4FrgUHgJkmDNd1uAU5GxApgKzCcxg4CQ8BKYC1wn6SuFjGHga0RMQCcTLEbzgHcCCyOiPcBHwL+RXVhMzOz+VXmiGUVMBkRhyLidWAEWFfTZx2wPS3vBK6WpNQ+EhGnI+JlYDLFqxszjVmdYpBi3tBijgDeKakb+HngdeBvSm8BMzPLqkxhWQocqXo9ldrq9omIaeA1oKfJ2EbtPcCpFKN2rkZz7AT+H/Aj4BXgv0XEido3IWmjpDFJY8eOHSvxts3MrB1lCovqtEXJPrnam82xCngDeA/QD/x7Se99W8eIbRFRiYhKb29vnVBmZpZDmcIyBSyver0MONqoTzoldTFwosnYRu3HgSUpRu1cjeb4JPDtiPhpRLwK/CVQKfG+zMxsDpQpLPuAgXS31iKKi/GjNX1GgQ1peT3wZEREah9Kd3T1AwPA3kYx05inUgxSzF0t5ngFWK3CO4FfBb5XfhOYmVlO3a06RMS0pE3AbqALeDAixiXdCYxFxCjwAPCQpEmKo4ihNHZc0g7gADAN3BoRbwDUi5mmvB0YkXQX8GyKTaM5KO4u+wPguxSny/4gIl5oe4uYmdmsqPij/8JSqVRibGys7fF9mx/LmI2Z2fw5fM91bY+VtD8iWl5q8CfvzcwsKxcWMzPLyoXFzMyycmExM7OsXFjMzCwrFxYzM8vKhcXMzLJyYTEzs6xcWMzMLCsXFjMzy8qFxczMsnJhMTOzrFxYzMwsKxcWMzPLyoXFzMyycmExM7OsXFjMzCwrFxYzM8vKhcXMzLJyYTEzs6xKFRZJayVNSJqUtLnO+sWSHknr90jqq1q3JbVPSFrTKqak/hTjYIq5qMQc75f0HUnjkl6UdFE7G8PMzGavZWGR1AXcC1wLDAI3SRqs6XYLcDIiVgBbgeE0dhAYAlYCa4H7JHW1iDkMbI2IAeBkit1sjm7g68C/jIiVwFXAT2e4HczMLJMyRyyrgMmIOBQRrwMjwLqaPuuA7Wl5J3C1JKX2kYg4HREvA5MpXt2YaczqFIMU84YWc1wDvBARzwNExP+NiDfKbwIzM8upTGFZChypej2V2ur2iYhp4DWgp8nYRu09wKkUo3auRnP8MhCSdkt6RtLv1nsTkjZKGpM0duzYsRJv28zM2lGmsKhOW5Tsk6u92RzdwK8Bn0r/fkLS1W/rGLEtIioRUent7a0TyszMcihTWKaA5VWvlwFHG/VJ1zwuBk40Gduo/TiwJMWonavZHH8REccj4m+Bx4EPlnhfZmY2B8oUln3AQLpbaxHFxfjRmj6jwIa0vB54MiIitQ+lO7r6gQFgb6OYacxTKQYp5q4Wc+wG3i/pF1LB+cfAgfKbwMzMcupu1SEipiVtovgF3gU8GBHjku4ExiJiFHgAeEjSJMVRxFAaOy5pB8Uv+mng1jMX1uvFTFPeDoxIugt4NsWmyRwnJX2JolgF8HhEPDarrWJmZm1T8Uf/haVSqcTY2Fjb4/s2u26Z2bnp8D3XtT1W0v6IqLTq50/em5lZVi4sZmaWlQuLmZll5cJiZmZZubCYmVlWLixmZpaVC4uZmWXlwmJmZlm5sJiZWVYuLGZmlpULi5mZZeXCYmZmWbmwmJlZVi4sZmaWlQuLmZll5cJiZmZZubCYmVlWLixmZpaVC4uZmWXlwmJmZlmVKiyS1kqakDQpaXOd9YslPZLW75HUV7VuS2qfkLSmVUxJ/SnGwRRzUas50vorJP1Y0mdmuhHMzCyfloVFUhdwL3AtMAjcJGmwptstwMmIWAFsBYbT2EFgCFgJrAXuk9TVIuYwsDUiBoCTKXbDOapsBb5V9o2bmdncKHPEsgqYjIhDEfE6MAKsq+mzDtielncCV0tSah+JiNMR8TIwmeLVjZnGrE4xSDFvaDEHkm4ADgHj5d+6mZnNhTKFZSlwpOr1VGqr2ycipoHXgJ4mYxu19wCnUozauerOIemdwO3A55u9CUkbJY1JGjt27FiLt2xmZu0qU1hUpy1K9snV3myOz1OcOvtxnfVvdYzYFhGViKj09vY262pmZrPQXaLPFLC86vUy4GiDPlOSuoGLgRMtxtZrPw4skdSdjkqq+zea4yPAeklfBJYAb0r6SUR8pcR7MzOzzMocsewDBtLdWosoLsaP1vQZBTak5fXAkxERqX0o3dHVDwwAexvFTGOeSjFIMXc1myMifj0i+iKiD/gy8F9dVMzMFk7LI5aImJa0CdgNdAEPRsS4pDuBsYgYBR4AHpI0SXEUMZTGjkvaARwApoFbI+INgHox05S3AyOS7gKeTbFpNIeZmXUWFQcJF5ZKpRJjY2Ntj+/b/FjGbMzM5s/he65re6yk/RFRadXPn7w3M7OsXFjMzCwrFxYzM8vKhcXMzLJyYTEzs6xcWMzMLCsXFjMzy8qFxczMsnJhMTOzrFxYzMwsKxcWMzPLyoXFzMyycmExM7OsXFjMzCwrFxYzM8vKhcXMzLJyYTEzs6xcWMzMLCsXFjMzy8qFxczMsipVWCStlTQhaVLS5jrrF0t6JK3fI6mvat2W1D4haU2rmJL6U4yDKeaiZnNI+pik/ZJeTP+ubndjmJnZ7LUsLJK6gHuBa4FB4CZJgzXdbgFORsQKYCswnMYOAkPASmAtcJ+krhYxh4GtETEAnEyxG84BHAd+MyLeB2wAHprZJjAzs5zKHLGsAiYj4lBEvA6MAOtq+qwDtqflncDVkpTaRyLidES8DEymeHVjpjGrUwxSzBuazRERz0bE0dQ+DlwkaXHZDWBmZnmVKSxLgSNVr6dSW90+ETENvAb0NBnbqL0HOJVi1M7VaI5q/xR4NiJOl3hfZmY2B7pL9FGdtijZp1F7vYLWrH/LPCStpDg9dk2dfkjaCGwEuOKKK+p1MTOzDMocsUwBy6teLwOONuojqRu4GDjRZGyj9uPAkhSjdq5GcyBpGfAo8NsR8YN6byIitkVEJSIqvb29Jd62mZm1o0xh2QcMpLu1FlFcjB+t6TNKceEcYD3wZEREah9Kd3T1AwPA3kYx05inUgxSzF3N5pC0BHgM2BIRfzmTN29mZvm1LCzpesYmYDfwErAjIsYl3Snp+tTtAaBH0iRwG7A5jR0HdgAHgG8Dt0bEG41ipli3A7elWD0pdsM5UpwVwH+W9Fz6+sU2t4eZmc2SioOEC0ulUomxsbG2x/dtfixjNmZm8+fwPde1PVbS/oiotOrnT96bmVlWLixmZpaVC4uZmWXlwmJmZlm5sJiZWVYuLGZmlpULi5mZZeXCYmZmWbmwmJlZVi4sZmaWlQuLmZll5cJiZmZZubCYmVlWLixmZpaVC4uZmWXlwmJmZlm5sJiZWVYuLGZmlpULi5mZZeXCYmZmWZUqLJLWSpqQNClpc531iyU9ktbvkdRXtW5Lap+QtKZVTEn9KcbBFHNRu3OYmdn8a1lYJHUB9wLXAoPATZIGa7rdApyMiBXAVmA4jR0EhoCVwFrgPkldLWIOA1sjYgA4mWLPeI6ZbggzM8ujzBHLKmAyIg5FxOvACLCups86YHta3glcLUmpfSQiTkfEy8Bkilc3ZhqzOsUgxbyhzTnMzGwBlCksS4EjVa+nUlvdPhExDbwG9DQZ26i9BziVYtTONdM5zMxsAXSX6KM6bVGyT6P2egWtWf925jg7QWkjsDG9/LGkiTrjyroMOD6L8XPFec2M85oZ5zUzHZmXhmeV198r06lMYZkClle9XgYcbdBnSlI3cDFwosXYeu3HgSWSutNRSXX/dub4mYjYBmwr8X5bkjQWEZUcsXJyXjPjvGbGec3MhZxXmVNh+4CBdLfWIooL5aM1fUaBDWl5PfBkRERqH0p3dPUDA8DeRjHTmKdSDFLMXW3OYWZmC6DlEUtETEvaBOwGuoAHI2Jc0p3AWESMAg8AD0mapDiKGEpjxyXtAA4A08CtEfEGQL2YacrbgRFJdwHPpti0M4eZmc0/FX/020xI2phOrXUU5zUzzmtmnNfMXMh5ubCYmVlWfqSLmZll5cJSRdJnJIWky9JrSfr99LiYFyR9sKrvhvTYmYOSNlS1f0jSi2nM76cPcSLpUklPpP5PSLqkRD5fSPM+J+nPJL2nQ/L6PUnfS3M/KmlJ1bo5f4RPk7xulDQu6U1JlZp1C5ZXWY1yyUnSg5JelfTdqra6+0DO/axEXsslPSXppfQ9/LedkJukiyTtlfR8yuvzqX3G+8dM98GS261L0rOSvtlJeRER/ipOBy6nuJngr4DLUtvHgW9RfFbmV4E9qf1S4FD695K0fElatxf4aBrzLeDa1P5FYHNa3gwMl8jp3VXL/wb4WofkdQ3QnZaHz4yheDzP88BioB/4AcXNGV1p+b3AotRnMI3ZAQyl5a8B/yotf7rq/Q4Bj5TI6+8DvwI8DVSq2hc0r5L7X8NcMu/nVwIfBL5b1VZ3H8i5n5XI63Lgg2n57wDfT9+3Bc0t9X1XWn4HsCfNN6P9o519sOR2uw14GPhmO/vtnOWVe8c9V78oHhPzD4DDvFVY7gduquozkX4AbgLur2q/P7VdDnyvqv1n/c6MrfohmphhfluAr3ZgXp8AvlGV45aqdbvTD/JHgd0172VL+qE9zltF6mf9zoxNy92pn0rm9DRnF5aOyKtFznVzmaN9vY+zC0vdfSDnftZGjruAj3VSbsAvAM8AH5np/jHTfbBkPsuA/0XxGKxvtrPfzkVeEeFTYQCSrgd+GBHP16ya6SNplqbl2naAX4qIHwGkf3+xZG53SzoCfAr4bKfkVeV3KP76ayevdh7h045OzatMjvOh0T6Qcz8rLZ2m+YcURwcLnls63fQc8CrwBMVf8rkePTWb7/uXgd8F3kyvcz4Sa1b7Y5lP3p8XJP058HfrrLoD+I8Up3feNqxOW7PHyJR6vEzZvCJiV0TcAdwhaQuwCfhcJ+SV+txB8dmhb5wZ1mCebI/wKZNXHfPxaKHZmqu4szHn+9nbJpTeBfwJ8O8i4m+aXAaZt9yi+FzcB1RcS3yU4pRro1i5Hm/VlKTfAF6NiP2Srmox97zldcYFU1gi4p/Ua5f0Popzi8+nnXgZ8IykVTR+XMwUcFVN+9OpfVmd/gB/LenyiPiRpMsp/vppmFcdDwOPURSWBc8rXRT9DeDqSMfKTfKiQfuMH+Ezg+1VbT4eLTRbpR5NNEfq7gNNcmpnP2tJ0jsoiso3IuJPOyk3gIg4JelpimssOR891c73/R8B10v6OHAR8G6KI5iFzqvQzvnP8/mLs6+xXMfZFwj3pvZLgZcpLg5ekpYvTev2pb5nLhB+PLX/HmdfhPxiiVwGqpb/NbCzQ/JaS/Gkg96a9pWcfSHwEMVFwO603M9bFwJXpjF/zNkXGz+dlm/l7IuNO2bwPXyas6+xdEReLXJumMsc7ON9nH2Npe4+kHM/K5GTgD8CvlzTvqC5Ab3AkrT888D/pviDakb7Rzv74Ay+n1fx1sX7jshr3n9xd/oXZxcWUfyHZD8AXuTsX1a/Q/F/v0wC/7yqvQJ8N435Cm99CLWH4kLbwfTvpSVy+ZMU6wXgfwJLOySvSYrzr8+lr69VrbsjzTFB1V03FHfxfD+tu6Oq/b0Ud+tMph+Kxan9ovR6Mq1/b4m8PkHxF9hp4K85++LjguU1g32vbi6Z9+//AfwI+GnaVrc02gdy7mcl8vo1ilMtL1TtVx9f6NyA91M8WuqFNPaz7e4fM90HZ/A9vYq3CktH5OVP3puZWVa+K8zMzLJyYTEzs6xcWMzMLCsXFjMzy8qFxczMsnJhMTOzrFxYzMwsKxcWMzPL6v8DtSa0I/aKUyMAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"l1_X=particles_np['l1'][:,0];\n",
"l1_Y=particles_np['l1'][:,1];\n",
"l1_Z=particles_np['l1'][:,2];\n",
"l1_E=particles_np['l1'][:,3];\n",
"\n",
"l2_X=particles_np['l2'][:,0];\n",
"l2_Y=particles_np['l2'][:,1];\n",
"l2_Z=particles_np['l2'][:,2];\n",
"l2_E=particles_np['l2'][:,3];\n",
"\n",
"plt.hist(l1_X,bins=40, density=True);"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0.0, 0.0, 0.0, 80000.0)"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"((l1_X[:]+l2_X[:]).mean(), (l1_Y[:]+l2_Y[:]).mean(), (l1_Z[:]+l2_Z[:]).mean(), (l1_E[:]+l2_E[:]).mean())"
]
},
{
"cell_type": "code",
"execution_count": 9,
"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": 10,
"metadata": {},
"outputs": [],
"source": [
"class GaussianSampleAndWeights():\n",
"\n",
" def __call__(self, n_to_produce, limits, dtype):\n",
" \n",
" n_to_produce = tf.cast(n_to_produce, dtype=tf.int32)\n",
"\n",
" weights, phase_space_sample = Z.generate_tensor(n_events=n_to_produce) #check order \n",
" \n",
" phase_space_sample_tensor = tf.concat([phase_space_sample[\"l1\"],phase_space_sample[\"l2\"]], axis=1)\n",
" \n",
" weights_max = None\n",
" \n",
" thresholds = tf.random_uniform(shape=(n_to_produce,), dtype=dtype)\n",
" \n",
" return phase_space_sample_tensor, thresholds, weights, weights_max, n_to_produce\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"class dGamma(zfit.pdf.ZPDF):\n",
" \n",
" #_PARAMS = ['ml']\n",
" _PARAMS = []\n",
" _N_OBS = 8\n",
"\n",
" def _unnormalized_pdf(self, x):\n",
" \n",
" #ml = self.params['ml']\n",
" \n",
" p1x, p1y, p1z, p1E, p2x, p2y, p2z, p2E =x.unstack_x()\n",
" \n",
" p1 = [p1x, p1y, p1z, p1E]\n",
" p2 = [p2x, p2y, p2z, p2E]\n",
" z = [0., 0., 1.]\n",
" \n",
" costheta_l= get_costheta_l(p1, z)\n",
" #q2 = sp(p1+p2,p1+p2)\n",
"\n",
" pdf = (1+costheta_l+ztf.square(costheta_l))#*beta_l(q2, 0.)\n",
" \n",
" return pdf\n"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"#m_lepton = zfit.Parameter('ml', mmu)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"upper_vec=np.sqrt(np.square(mZ)/4-np.square(mmu))\n",
"lower_vec=-np.sqrt(np.square(mZ)/4-np.square(mmu))"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"lower = ((lower_vec, lower_vec, lower_vec, 0., lower_vec, lower_vec, lower_vec, 0., ),)\n",
"upper = ((upper_vec, upper_vec, upper_vec, mZ/2, upper_vec, upper_vec, upper_vec, mZ/2, ),)\n",
"\n",
"obs = zfit.Space([\"p1x\",\"p1y\",\"p1z\",\"p1E\", \"p2x\",\"p2y\", \"p2z\",\"p2E\"], limits=(lower,upper))\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"pdf = dGamma(obs=obs)\n",
"#pdf = dGamma(obs=obs, ml = m_lepton)\n",
"pdf._sample_and_weights=GaussianSampleAndWeights"
]
},
{
"cell_type": "code",
"execution_count": 20,
"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=100000)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"for i in range(1):\n",
" sampler.resample()\n",
"a=sampler.to_pandas() "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.hist(a['p1z']/(mZ/2),bins=40);"
]
},
{
"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.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Davide Lancierini