{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Import"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"c:\\users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\util\\execution.py:53: 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 numpy as np\n",
"from pdg_const import pdg\n",
"import matplotlib\n",
"import matplotlib.pyplot as plt\n",
"import pickle as pkl\n",
"import sys\n",
"import time\n",
"from helperfunctions import display_time, prepare_plot\n",
"import cmath as c\n",
"import scipy.integrate as integrate\n",
"from scipy.optimize import fminbound\n",
"from array import array as arr\n",
"import collections\n",
"from itertools import compress\n",
"import tensorflow as tf\n",
"import zfit\n",
"from zfit import ztf\n",
"from IPython.display import clear_output"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Build model and graphs\n",
"## Create graphs"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def formfactor( q2, subscript): #returns real value\n",
" #check if subscript is viable\n",
"\n",
" if subscript != \"0\" and subscript != \"+\" and subscript != \"T\":\n",
" raise ValueError('Wrong subscript entered, choose either 0, + or T')\n",
"\n",
" #get constants\n",
"\n",
" mK = ztf.constant(pdg['Ks_M'])\n",
" mbstar0 = ztf.constant(pdg[\"mbstar0\"])\n",
" mbstar = ztf.constant(pdg[\"mbstar\"])\n",
" b0 = ztf.constant(pdg[\"b0\"])\n",
" bplus = ztf.constant(pdg[\"bplus\"])\n",
" bT = ztf.constant(pdg[\"bT\"])\n",
"\n",
" mmu = ztf.constant(pdg['muon_M'])\n",
" mb = ztf.constant(pdg['bquark_M'])\n",
" ms = ztf.constant(pdg['squark_M'])\n",
" mB = ztf.constant(pdg['Bplus_M'])\n",
"\n",
" #N comes from derivation in paper\n",
"\n",
" N = 3\n",
"\n",
" #some helperfunctions\n",
"\n",
" tpos = (mB - mK)**2\n",
" tzero = (mB + mK)*(ztf.sqrt(mB)-ztf.sqrt(mK))**2\n",
"\n",
" z_oben = ztf.sqrt(tpos - q2) - ztf.sqrt(tpos - tzero)\n",
" z_unten = ztf.sqrt(tpos - q2) + ztf.sqrt(tpos - tzero)\n",
" z = tf.divide(z_oben, z_unten)\n",
"\n",
" #calculate f0\n",
"\n",
" if subscript == \"0\":\n",
" prefactor = 1/(1 - q2/(mbstar0**2))\n",
" _sum = 0\n",
"\n",
" for i in range(N):\n",
" _sum += b0[i]*(tf.pow(z,i))\n",
"\n",
" return tf.complex(prefactor * _sum, ztf.constant(0.0))\n",
"\n",
" #calculate f+ or fT\n",
"\n",
" else:\n",
" prefactor = 1/(1 - q2/(mbstar**2))\n",
" _sum = 0\n",
"\n",
" if subscript == \"T\":\n",
" b = bT\n",
" else:\n",
" b = bplus\n",
"\n",
" for i in range(N):\n",
" _sum += b[i] * (tf.pow(z, i) - ((-1)**(i-N)) * (i/N) * tf.pow(z, N))\n",
"\n",
" return tf.complex(prefactor * _sum, ztf.constant(0.0))\n",
"\n",
"def resonance(q, _mass, width, phase, scale):\n",
"\n",
" q2 = tf.pow(q, 2)\n",
"\n",
" mmu = ztf.constant(pdg['muon_M'])\n",
"\n",
" p = 0.5 * ztf.sqrt(q2 - 4*(mmu**2))\n",
"\n",
" p0 = 0.5 * ztf.sqrt(_mass**2 - 4*mmu**2)\n",
"\n",
" gamma_j = tf.divide(p, q2) * _mass * width / p0\n",
"\n",
" #Calculate the resonance\n",
"\n",
" _top = tf.complex(_mass * width, ztf.constant(0.0))\n",
"\n",
" _bottom = tf.complex(_mass**2 - q2, -_mass*gamma_j)\n",
"\n",
" com = _top/_bottom\n",
"\n",
" #Rotate by the phase\n",
"\n",
" r = tf.abs(com)\n",
"\n",
" _phase = tf.angle(com)\n",
"\n",
" _phase += phase\n",
"\n",
" x = tf.cos(phase)*r\n",
" y = tf.sin(phase)*r\n",
"\n",
" com = tf.complex(scale* x, scale * y)\n",
"\n",
" return com\n",
"\n",
"def bifur_gauss(q, mean, sigma_L, sigma_R, scale):\n",
"\n",
" _exp = tf.where(q < mean, ztf.exp(- tf.pow((q-mean),2) / (2 * sigma_L**2)), ztf.exp(- tf.pow((q-mean),2) / (2 * sigma_R**2)))\n",
"\n",
" #Scale so the total area under curve is 1 and the top of the cusp is continuous\n",
"\n",
" dgamma = scale*_exp/(ztf.sqrt(2*np.pi))*2*(sigma_L*sigma_R)/(sigma_L+sigma_R)\n",
"\n",
" com = ztf.complex(dgamma, ztf.constant(0.0))\n",
"\n",
" return com\n",
"\n",
"def axiv_nonres(q):\n",
"\n",
" GF = ztf.constant(pdg['GF'])\n",
" alpha_ew = ztf.constant(pdg['alpha_ew'])\n",
" Vtb = ztf.constant(pdg['Vtb'])\n",
" Vts = ztf.constant(pdg['Vts'])\n",
" C10eff = ztf.constant(pdg['C10eff'])\n",
"\n",
" mmu = ztf.constant(pdg['muon_M'])\n",
" mb = ztf.constant(pdg['bquark_M'])\n",
" ms = ztf.constant(pdg['squark_M'])\n",
" mK = ztf.constant(pdg['Ks_M'])\n",
" mB = ztf.constant(pdg['Bplus_M'])\n",
"\n",
" q2 = tf.pow(q, 2)\n",
"\n",
" #Some helperfunctions\n",
"\n",
" beta = ztf.sqrt(tf.abs(1. - 4. * mmu**2. / q2))\n",
"\n",
" kabs = ztf.sqrt(mB**2. +tf.pow(q2, 2)/mB**2. + mK**4./mB**2. - 2. * (mB**2. * mK**2. + mK**2. * q2 + mB**2. * q2) / mB**2.)\n",
"\n",
" #prefactor in front of whole bracket\n",
"\n",
" prefactor1 = GF**2. *alpha_ew**2. * (tf.abs(Vtb*Vts))**2. * kabs * beta / (128. * np.pi**5.)\n",
"\n",
" #left term in bracket\n",
"\n",
" bracket_left = 2./3. * kabs**2. * beta**2. *tf.abs(tf.complex(C10eff, ztf.constant(0.0))*formfactor(q2, \"+\"))**2.\n",
"\n",
" #middle term in bracket\n",
"\n",
" _top = 4. * mmu**2. * (mB**2. - mK**2.) * (mB**2. - mK**2.)\n",
"\n",
" _under = q2 * mB**2.\n",
"\n",
" bracket_middle = _top/_under *tf.pow(tf.abs(tf.complex(C10eff, ztf.constant(0.0)) * formfactor(q2, \"0\")), 2)\n",
"\n",
" #Note sqrt(q2) comes from derivation as we use q2 and plot q\n",
"\n",
" return prefactor1 * (bracket_left + bracket_middle) * 2 *ztf.sqrt(q2)\n",
"\n",
"def vec(q, funcs):\n",
" \n",
" q2 = tf.pow(q, 2)\n",
"\n",
" GF = ztf.constant(pdg['GF'])\n",
" alpha_ew = ztf.constant(pdg['alpha_ew'])\n",
" Vtb = ztf.constant(pdg['Vtb'])\n",
" Vts = ztf.constant(pdg['Vts'])\n",
" C7eff = ztf.constant(pdg['C7eff'])\n",
"\n",
" mmu = ztf.constant(pdg['muon_M'])\n",
" mb = ztf.constant(pdg['bquark_M'])\n",
" ms = ztf.constant(pdg['squark_M'])\n",
" mK = ztf.constant(pdg['Ks_M'])\n",
" mB = ztf.constant(pdg['Bplus_M'])\n",
"\n",
" #Some helperfunctions\n",
"\n",
" beta = ztf.sqrt(tf.abs(1. - 4. * mmu**2. / q2))\n",
"\n",
" kabs = ztf.sqrt(mB**2. + tf.pow(q2, 2)/mB**2. + mK**4./mB**2. - 2 * (mB**2 * mK**2 + mK**2 * q2 + mB**2 * q2) / mB**2)\n",
"\n",
" #prefactor in front of whole bracket\n",
"\n",
" prefactor1 = GF**2. *alpha_ew**2. * (tf.abs(Vtb*Vts))**2 * kabs * beta / (128. * np.pi**5.)\n",
"\n",
" #right term in bracket\n",
"\n",
" prefactor2 = kabs**2 * (1. - 1./3. * beta**2)\n",
"\n",
" abs_bracket = tf.abs(c9eff(q, funcs) * formfactor(q2, \"+\") + tf.complex(2.0 * C7eff * (mb + ms)/(mB + mK), ztf.constant(0.0)) * formfactor(q2, \"T\"))**2\n",
"\n",
" bracket_right = prefactor2 * abs_bracket\n",
"\n",
" #Note sqrt(q2) comes from derivation as we use q2 and plot q\n",
"\n",
" return prefactor1 * bracket_right * 2 * ztf.sqrt(q2)\n",
"\n",
"def c9eff(q, funcs):\n",
"\n",
" C9eff_nr = tf.complex(ztf.constant(pdg['C9eff']), ztf.constant(0.0))\n",
"\n",
" c9 = C9eff_nr\n",
"\n",
" c9 = c9 + funcs\n",
"\n",
" return c9"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"def G(y):\n",
" \n",
" def inner_rect_bracket(q):\n",
" return tf.log(ztf.to_complex((1+tf.sqrt(q))/(1-tf.sqrt(q)))-tf.complex(ztf.constant(0), -1*ztf.constant(np.pi))) \n",
" \n",
" def inner_right(q):\n",
" return ztf.to_complex(2 * tf.atan(1/tf.sqrt(-q)))\n",
" \n",
" big_bracket = tf.where(y > ztf.const(0.0), inner_rect_bracket(y), inner_right(y))\n",
" \n",
" return ztf.to_complex(tf.sqrt(tf.abs(y))) * big_bracket\n",
"\n",
"def h_S(m, q):\n",
" \n",
" return tf.constant(2) - G(tf.constant(1) - 4*tf.pow(m, 2) / tf.pow(q, 2))\n",
"\n",
"def h_P(m,q):\n",
" \n",
" return 2/3 + (1 - (tf.constant(1) - 4*tf.pow(m, 2) / tf.pow(q, 2))) * h_S(m,q)\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Build pdf"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"class total_pdf(zfit.pdf.ZPDF):\n",
" _N_OBS = 1 # dimension, can be omitted\n",
" _PARAMS = ['jpsi_mass', 'jpsi_scale', 'jpsi_phase', 'jpsi_width',\n",
" 'psi2s_mass', 'psi2s_scale', 'psi2s_phase', 'psi2s_width',\n",
" 'cusp_mass', 'sigma_L', 'sigma_R', 'cusp_scale'\n",
" ] # the name of the parameters\n",
"\n",
" def _unnormalized_pdf(self, x):\n",
" \n",
" x = x.unstack_x()\n",
"\n",
" def jpsi_res(q):\n",
" return resonance(q, _mass = self.params['jpsi_mass'], scale = self.params['jpsi_scale'], phase = self.params['jpsi_phase'], width = self.params['jpsi_width'])\n",
"\n",
" def psi2s_res(q):\n",
" return resonance(q, _mass = self.params['psi2s_mass'], scale = self.params['psi2s_scale'], phase = self.params['psi2s_phase'], width = self.params['psi2s_width'])\n",
"\n",
" def cusp(q):\n",
" return bifur_gauss(q, mean = self.params['cusp_mass'], sigma_L = self.params['sigma_L'], sigma_R = self.params['sigma_R'], scale = self.params['cusp_scale'])\n",
"\n",
" funcs = jpsi_res(x) + psi2s_res(x) + cusp(x)\n",
"\n",
" vec_f = vec(x, funcs)\n",
"\n",
" axiv_nr = axiv_nonres(x)\n",
"\n",
" tot = vec_f + axiv_nr\n",
"\n",
" return tot"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Load data"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"x_min = 2*pdg['muon_M']\n",
"x_max = (pdg[\"Bplus_M\"]-pdg[\"Ks_M\"]-0.1)\n",
"\n",
"obs = zfit.Space('q', limits = (x_min, x_max))\n",
"\n",
"with open(r\"./data/slim_points/slim_points_toy_0_range({0}-{1}).pkl\".format(int(x_min), int(x_max)), \"rb\") as input_file:\n",
" part_set = pkl.load(input_file)\n",
"\n",
"x_part = part_set['x_part']\n",
"\n",
"x_part = x_part.astype('float64')\n",
"\n",
"data = zfit.data.Data.from_numpy(array=x_part, obs=obs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup parameters"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"WARNING:tensorflow:From c:\\users\\sa_li\\.conda\\envs\\rmd\\lib\\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": [
"#jpsi\n",
"\n",
"jpsi_mass, jpsi_width, jpsi_phase, jpsi_scale = pdg[\"jpsi\"]\n",
"jpsi_scale *= pdg[\"factor_jpsi\"]\n",
"\n",
"jpsi_m = zfit.Parameter(\"jpsi_m\", ztf.constant(jpsi_mass), floating = False)\n",
"jpsi_w = zfit.Parameter(\"jpsi_w\", ztf.constant(jpsi_width), floating = False)\n",
"jpsi_p = zfit.Parameter(\"jpsi_p\", ztf.constant(jpsi_phase))\n",
"jpsi_s = zfit.Parameter(\"jpsi_s\", ztf.constant(jpsi_scale))\n",
"\n",
"#psi2s\n",
"\n",
"psi2s_mass, psi2s_width, psi2s_phase, psi2s_scale = pdg[\"psi2s\"]\n",
"psi2s_scale *= pdg[\"factor_psi2s\"]\n",
"\n",
"psi2s_m = zfit.Parameter(\"psi2s_m\", ztf.constant(psi2s_mass), floating = False)\n",
"psi2s_w = zfit.Parameter(\"psi2s_w\", ztf.constant(psi2s_width), floating = False)\n",
"psi2s_p = zfit.Parameter(\"psi2s_p\", ztf.constant(psi2s_phase))\n",
"psi2s_s = zfit.Parameter(\"psi2s_s\", ztf.constant(psi2s_scale))\n",
"\n",
"#cusp\n",
"\n",
"cusp_mass, sigma_R, sigma_L, cusp_scale = 3550, 3e-7, 200, 0\n",
"\n",
"cusp_m = zfit.Parameter(\"cusp_m\", ztf.constant(cusp_mass), floating = False)\n",
"sig_L = zfit.Parameter(\"sig_L\", ztf.constant(sigma_L), floating = False)\n",
"sig_R = zfit.Parameter(\"sig_R\", ztf.constant(sigma_R), floating = False)\n",
"cusp_s = zfit.Parameter(\"cusp_s\", ztf.constant(cusp_scale), floating = False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup pdf"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"total_f = total_pdf(obs=obs, jpsi_mass = jpsi_m, jpsi_scale = jpsi_s, jpsi_phase = jpsi_p, jpsi_width = jpsi_w,\n",
" psi2s_mass = psi2s_m, psi2s_scale = psi2s_s, psi2s_phase = psi2s_p, psi2s_width = psi2s_w,\n",
" cusp_mass = cusp_m, sigma_L = sig_L, sigma_R = sig_R, cusp_scale = cusp_s)\n",
"\n",
"# print(total_pdf.obs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Test if graphs actually work and compute values"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"def total_test_tf(xq):\n",
"\n",
" def jpsi_res(q):\n",
" return resonance(q, jpsi_m, jpsi_s, jpsi_p, jpsi_w)\n",
"\n",
" def psi2s_res(q):\n",
" return resonance(q, psi2s_m, psi2s_s, psi2s_p, psi2s_w)\n",
"\n",
" def cusp(q):\n",
" return bifur_gauss(q, cusp_m, sig_L, sig_R, cusp_s)\n",
"\n",
" funcs = jpsi_res(xq) + psi2s_res(xq) + cusp(xq)\n",
"\n",
" vec_f = vec(xq, funcs)\n",
"\n",
" axiv_nr = axiv_nonres(xq)\n",
"\n",
" tot = vec_f + axiv_nr\n",
" \n",
" return tot\n",
"\n",
"def jpsi_res(q):\n",
" return resonance(q, jpsi_m, jpsi_s, jpsi_p, jpsi_w)\n",
"\n",
"# calcs = zfit.run(total_test_tf(x_part))\n",
"\n",
"test_q = np.linspace(x_min, x_max, 2000000)\n",
"\n",
"probs = total_f.pdf(test_q)\n",
"\n",
"calcs_test = zfit.run(probs)\n",
"res_y = zfit.run(jpsi_res(test_q))"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.09\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAD8CAYAAACl69mTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8FdX5+PHPk4WEkBAghLAkkABhCbIIEXEtdSm44oISrYoLoha/Xfy1Vtpv1Vrtt7YWq9YNBcUVKHWJFhcUtypbkEXCGvZAgEBIAtmX5/fHHTDEhFwgYe7yvF+vvJg798yZZ4bkPvfMmXNGVBVjjDGmJYS4HYAxxpjAZUnGGGNMi7EkY4wxpsVYkjHGGNNiLMkYY4xpMZZkjDHGtBivkoyIjBaRdSKSIyL3NfB+hIjMct5fJCLJdd6b7KxfJyKjjqHOp0TkoDf7MMYY45uaTDIiEgo8DVwEpAHXiUhavWK3AftVtTfwOPCos20akAEMAEYDz4hIaFN1ikg60M6bfRhjjPFd3rRkhgM5qrpJVSuBmcCYemXGADOc5TnA+SIizvqZqlqhqpuBHKe+Rut0EtDfgHu93IcxxhgfFeZFmW7A9jqvc4HTGyujqtUiUgTEOesX1tu2m7PcWJ13A5mqmlcvhzS2j711C4nIRGAiQJs2bYb169fPi0M0xhyy92AFeUXlxLVpRdd2rd0O5wjVtcqavGIABnaLdTmawLV06dK9qhrfHHV5k2Qaai3Un4umsTKNrW+oBaUi0hW4Bhh5nHGgqlOBqQDp6emalZXVwGbGmMZM++9m/vT+am4+M5kHLx/gdjhH2HuwgvSHPwEg6y+XuBxN4BKRrc1VlzeXy3KBpDqvE4GdjZURkTAgFig4yraNrT8V6A3kiMgWIEpEcprYhzGmGfnyNeham2vR73iTZJYAqSKSIiKt8HTkZ9YrkwmMd5bHAvPVM/NmJpDh3BmWAqQCixurU1X/o6qdVTVZVZOBUqej/2j7MMY0I5/u6bS/eL/T5OUyp//jbuAjIBSYrqrZIvIQkKWqmcA04FWn1VGAJ2nglJsNrAaqgUmqWgPQUJ1NhNLgPowxLcMXv8PV+l5Ipgne9MmgqnOBufXW3V9nuRxPX0pD2z4CPOJNnQ2UifZmH8eiqqqK3NxcysvLT7Qq10VGRpKYmEh4eLjboZgAcqgh44uf5+qTUZmj8SrJBJLc3FxiYmJITk7Gn++AVlX27dtHbm4uKSkpbodjzElhLRn/E3TTypSXlxMXF+fXCQZARIiLiwuIFpkx3qq1LON3gi7JAH6fYA4JlOMwvuXQ75UPdskYPxSUScYY0zhf/u5itzD7H0syPurzzz/n0ksvBaCiooILLriAIUOGMGvWLJcjM8HCFzvZLcf4n6Dr+PdHy5Yto6qqiuXLl7sdigkCh+8u88EPdGvJ+B9rybhgy5Yt9OvXj/HjxzNo0CDGjh1LaWkpH374If369ePss8/mrbfeAmDPnj3ccMMNLF++nCFDhrBx40aXozfGPdbv73+CuiXzx/eyWb2zuFnrTOvalgcua3q+p3Xr1jFt2jTOOussbr31VqZMmcLzzz/P/Pnz6d27N+PGjQOgU6dOvPjiizz22GO8//77zRqrMQ059Dnui30zzdGSKausYc2uYtbmHWDrvhL2l1ZSWllD6/BQYluHkxLfhgFdYxnYLZbQEB88CX4mqJOMm5KSkjjrrLMAuOGGG3jyySdJSUkhNTX18LqpU6e6GaIJUr58Raqyuva4tss/UMF7K3by2bo9LNpccLieVmEhdIhqRetWoVRU1VBQWkl5lee99lHhjBrQmRtG9OAUm/H5uAV1kvGmxdFS6t9+XFRUZLckG5/gi9PJHFJV432SUVW+ztnHy99s4bN1e6ipVXp3iubGET0Y0TOO/l1i6Nau9RF/d7W1ys6iMr7dVshna/fw7vKdzFyynZF94/ndxf3pkxDTEocV0II6ybhp27ZtLFiwgDPOOIM333yTCy64gOeff56NGzfSq1cv3nzzTbdDNEHKd1OM9y2ZFdsLefTDtXyzcR/xMRFMOCeFa4Yl0rvT0ZNESIiQ2D6KxPZRXD64Kw9eXsWbi7fx9Gc5XPTEV9z94978z3m9CQu17mxvWZJxSf/+/ZkxYwZ33HEHqampPPHEEwwbNoxLLrmEjh07cvbZZ7Nq1Sq3wzRByIcbMlTVHD24jfkH+fvH65j73S7i2rTigcvSuP707kSEhR7X/mJbh3Pnj3oxLj2Jh95fzROfbuDrnL08f+Mw4qIjjqvOYGNJxiUhISE899xzR6wbPXo0a9eu/UHZkSNHMnLkyJMUmQl2PpxjqKiuObysqocvde0qKueJT9czOyuXyLAQfnlBKhPO6Ul0RPN8xLVv04rHxw1hZN947p2zkque/YaXbxlOSsc2zVJ/ILMkY4w5wqE+GfHBx5cVllYdXj5YUU1tLTzzRQ4vf72FWlVuHNGDu8/rTccWamWMGdKNxPZR3P5KFtdNXci/7jyDpA5RLbKvQGFJxgXJycl2Kcz4rEOXy3zxPpT9pZWHl//f7BUs3LSPAxXVXDmkG7+6sM9J+cAf1qM9r084nYypC7n+xYX8+84z6dQ2ssX366+CsvfKl++eORaBchzGt/jidDKHbNh9kIiwEHp3imbemt2c2asjc39+DlPGDTmpLYr+Xdryyq3D2XugkjtfW3rEZTxzpKBryURGRrJv3z6/n+7/0PNkIiPtG5RpXr763aWkoppP1+5mZN94nr5+KJU1tUS1cu8jbHBSOx67ZjCT3viWBzOz+b+rBrkWiy8LuiSTmJhIbm4u+fn5bodywg49GdOY5uSLOaayupZ756xkX0klt5/Tk7DQEJ+4jfiSQV1YtbMXz36+kXNT47loYBe3Q/I5XiUZERkNPAGEAi+q6l/qvR8BvAIMA/YB41R1i/PeZOA2oAb4uap+dLQ6RWQakI5nnr71wM2qelBEbgb+BuxwdvtPVX3xWA84PDzcniRpzFH4WkvmYEU1d722lK827OV3F/cjPbmD2yEd4Z4L+/B1zl5+9/Z3DEtuT6cYu7pQV5NfBUQkFHgauAhIA64TkbR6xW4D9qtqb+Bx4FFn2zQgAxgAjAaeEZHQJur8laoOVtVBwDbg7jr7maWqQ5yfY04wxpimHeqT8YWLyXsPVnDd1IV8s3Effx07iInn9nI7pB8IDw1hyrWDKa2s4Q/v2A099XnT3hwO5KjqJlWtBGYCY+qVGQPMcJbnAOeLp8NjDDBTVStUdTOQ49TXaJ2qWgzgbN8a32y9GxOwfKUls3VfCVc/+w0b9hzghZuGcW16ktshNap3pxh+fn4qH2Xv5vN1e9wOx6d4k2S6AdvrvM511jVYRlWrgSIg7ijbHrVOEXkJ2AX0A56qU+5qEVkpInNExHd/44wxJ2TVjiKufvYbisqqeOP2EZzXL8HtkJo04ZwUUjq24Y/vrba7zerwJsk01Gqu/12nsTLHut6zoHoL0BVYA4xzVr8HJDuX0T7h+5bTkYGITBSRLBHJCoTOfWNONrdvjf9qQz7jnl9ARFgoc+48k6Hd27saj7ciwkJ54LI0Nu8t4dUFW90Ox2d4k2RygbqthkRgZ2NlRCQMiAUKjrJtk3Wqag0wC7jaeb1PVSuct1/Ac5PBD6jqVFVNV9X0+Ph4Lw7PGFOXmznmnWU7uOWlJSR1iOKtn51J707R7gVzHEb27cQ5qR155vONlFRUux2OT/AmySwBUkUkRURa4enIz6xXJhMY7yyPBear5+tQJpAhIhEikgKkAosbq1M8esPhPpnLgLXO67r3Bl6Op5VjjGlmbuWYF77cxC9nLSc9uT2z7zyDBD8dRX/PhX0oKKnk5W+2uB2KT2jyFmZVrRaRu4GP8NxuPF1Vs0XkISBLVTOBacCrIpKDpwWT4WybLSKzgdVANTDJaaHQSJ0hwAwRaYvnktoK4C4nlJ+LyOVOPQXAzc1yBowxR/h+WpmTc39Zba3yyNw1TPvvZi4Z2IUp4wYf96zJvuDU7u05v18nnv9iIzee0YO2keFuh+Qqr8bJqOpcYG69dffXWS4Hrmlk20eAR7yssxY4q5F6JgOTvYnXGHP8Tua0MhXVNfz6Xyt5b8VObj4zmfsvTSMkAB55/MsL+nDZP//LrMXbuf3cnm6H4yr3h8waY3zKyeqTKa+qYeIrS3lvxU7uu6gfD1wWGAkGYGBiLCN6dmD615uP6WmegciSjDHmCIdyTEveZVZWWcOEGVl8uSGfv1w1kDt/1Muv5xJsyMRze5JXVM5/Vua5HYqrLMkYY47kJJeWSjElFdXc8vJivtm4l8fGDiZjePcW2pO7RvbpRK/4Nrz4302u3xbuJksyxpgjtOTH4YHyKsZPX8zizQU8Pm4IVw8L3AleQ0KEm89MZtWOYr7bUeR2OK6xJGOMOcKhL93N/eW7qKyKm6YvZtn2Qp66bihjhtSfOCTwjDm1G5HhIby5eHvThQOUJRljzBFqD18ua74sU1hayY3TFrFqRxFPXz+USwYFx5T4bSPDuWRgVzKX7wjawZmWZIwxR/i+47956isoqeT6FxaxNu8Az90wjNGndG6eiv3E9acnUVJZw/sr60+UEhwsyRhjjlDbjB3/ew9WcP0LC9mYf5AXxqdzfn/fn+iyuQ3t3p7enaKZszTX7VBcYUnGGHOEmhonyZxgltlTXE7G1IVs2VfC9JtP40d9gnMuQRHhiiFdWbJlPzsKy9wO56SzJGOMOUKNk11qa48/y+wq8iSYnYVlvHzLcM7q3bG5wvNLlw3uCsD7K4LvkpklGWPMEWqc5HK8I9V3FJYxbuoC9hyo4JVbhzOiZ1xzhueXesS1YXBSOzItyRhjgt3hJHMcLZntBaWMe34BBSWVvHrbcNKTOzR3eH7r8sFdyd5ZzMb8g26HclJZkjHGHOFQkqk+xpbMlr0ljHt+AQfKq3l9wumc6icPGztZLh3UBRF4f0VwTTNjScYYc4TjuVy2Mf8g46YuoKyqhjduP51Bie1aKjy/ldA2klOT2vHJmt1uh3JSWZIxxhzh+yTj3eWyDbsPkDF1ITW1ysyJZzCga2xLhufXLkzrzHc7isgrCp67zCzJGGOOcOjuMm9aMmt3FZMxdSEAMyeOoG/nmBaNzd9dmNYJgE9WB09rxpKMMeYI1Yf7ZI7eklm1o4jrpi4kPDSEWRNH0LuTJZim9IqPJqVjGz62JGOMCVaHOvwrj9KSWbG9kOtfWEhUqzBm3TGCnvHRJys8vyYiXJiWwMJN+ygur3I7nJPCkowx5ghlVZ7kUl3bcJJZunU/N7y4iNiocGZOHEGPuDYnMzy/d0H/BKpqlP9u2Ot2KCeFV0lGREaLyDoRyRGR+xp4P0JEZjnvLxKR5DrvTXbWrxORUU3VKSLTRGSFiKwUkTkiEt3UPowxzae8sgaA0oqaH7y3aNM+bpq2iA7RrZg18QySOkSd7PD83qnd2xEdEcZXlmQ8RCQUeBq4CEgDrhORtHrFbgP2q2pv4HHgUWfbNCADGACMBp4RkdAm6vyVqg5W1UHANuDuo+3DGNO8yqs9yWV/aeUR679cn8/4lxbTOTaSWRPPoGu71m6E5/fCQ0M4o1ccX23ID4onZnrTkhkO5KjqJlWtBGYCY+qVGQPMcJbnAOeL54HdY4CZqlqhqpuBHKe+RutU1WIAZ/vWfD8ZbGP7MMY0ozKnJVNYVnX4duZ5q3czYUYWKR2jmXXHGXSOjXQzRL93TmpHcveXsXVfqduhtDhvkkw3oO5j3XKddQ2WUdVqoAiIO8q2R61TRF4CdgH9gKea2McRRGSiiGSJSFZ+fr4Xh2eMqauozNMhrQp5RWW8smALd762lP5d2zLz9hF0jI5wN8AAcE6qZ0bqr3IC/5KZN0mmodZC/TZeY2WOdb1nQfUWoCuwBhh3DHGgqlNVNV1V0+Pjg3NqcWOOV22tsq+kkvP6ecZzjHt+Ife/m83IPvG8PuF0YqPCXY4wMCTHRdGtXWu+Wh/4X4S9STK5QFKd14lA/alED5cRkTAgFig4yrZN1qmqNcAs4Oom9mGMaSZ7SyqoqVV+1Cee64Z7/kR/O7ofU29KJzoizOXoAoeIcG6fjizYuO+Y54jzN9781iwBUkUkBdiBpyP/+nplMoHxwAJgLDBfVVVEMoE3RGQKnpZJKrAYT6vkB3U6fSy9VDXHWb4MWHu0fRzncRtjGrA27wAAqQnRjD8z2d1gAtyZvTry5uLtZO8sZnBS4M711mSSUdVqEbkb+AgIBaararaIPARkqWomMA14VURy8LQuMpxts0VkNrAaqAYmOS0UGqkzBJghIm3xJKIVwF1OKA3uwxjTfL7euJewELH5x06C01M8j0FYvLkgoJOMBHJjID09XbOystwOwxi/UFpZzY/+9jmndG3LS7cMdzucoDDyb5+RmhDDCzelux3KEURkqao2S1A24t8YA8CTn+aQf6CCu8/r7XYoQWN4SgeWbCk4oUdd+zpLMsYYPlu3h+e+2Mi49CSG9bCnWZ4sw1PiKCytYsOewH1apiUZY4Lc+t0H+OXM5fTrHMMfxwxwO5yg8n2/zD6XI2k5lmSMCWI7Csu4adpiWoWF8MJN6USGh7odUlBJbN+azm0jWbxlv9uhtBhLMsYEqYKSSm6atoiSimpm3DLcJrt0gYgwPKUDizfvC9h5zCzJGBOESiurufXlJWzfX8aL49NJ69rW7ZCC1rAe7dldXEFeUbnbobQISzLGBJnyqhrufO1bVuYW8tR1p3J6zx9MAWhOokNjZJZvL3Q5kpZhScaYIFJZXcvdb3zLl+vz+ctVgxg1oLPbIQW9/l1iaBUaYknGGOPfqmtq+cXMZXyyZg9/uuIUrj0tqemNTIuLCAslrWtbSzLGGP9VU6vcM3sFH6zaxR8uTePGET3cDsnUMSSpHd/lFgXkZJmWZIwJcLW1ym//vZLMFTv57eh+3HZ2itshmXqGJLWjrKqG9bsDb1CmJRljApiq8r/vrmLO0lx+dUEf7hrZy+2QTAOGBHDnvyUZYwKUqvLH91bzxqJt/GxkL35+vs1J5qt6xEXRPiqc5dsDb1CmJRljApCq8pcP1vLyN1uYcHYKvxnVF88jmowvEhEGJbZjxfYit0NpdpZkjAlAj89bz/NfbuKmM3rw+0v6W4LxA6d0a0tO/kHKq2rcDqVZWZIxJsD8c/4GnpyfQ8ZpSTx42QBLMH5iQNdYamqV9bsPuB1Ks7IkY0wAefqzHB77eD1XndqNP185kJAQSzD+YoAztU/2zmKXI2leTT5+2RjjH/45fwOPfbyeK0/txt+uGWwJxs8ktY8iJiKM7J2B1S/jVUtGREaLyDoRyRGR+xp4P0JEZjnvLxKR5DrvTXbWrxORUU3VKSKvO+tXich0EQl31o8UkSIRWe783H8iB25MIKmbYB67ZjChlmD8TkiI0L9r24BryTSZZEQkFHgauAhIA64TkbR6xW4D9qtqb+Bx4FFn2zQgAxgAjAaeEZHQJup8HegHDARaAxPq7OcrVR3i/Dx0PAdsTKB56tMNhy+RWYLxbwO6tmVt3gFqAuhxzN60ZIYDOaq6SVUrgZnAmHplxgAznOU5wPni6W0cA8xU1QpV3QzkOPU1WqeqzlUHsBhIPLFDNCZwPfnpBv4+z5Ng/mYJxu+ldWlLWVUNm/eWuB1Ks/EmyXQDttd5neusa7CMqlYDRUDcUbZtsk7nMtmNwId1Vp8hIitE5AMRafA5sSIyUUSyRCQrPz/fi8Mzxj89+ekGpsxbz1VDLcEEigFdYwECql/GmyTT0G9u/bZcY2WOdX1dzwBfqupXzutvgR6qOhh4CninoWBVdaqqpqtqenx8fENFjPF7T3xSJ8GMtQQTKFITomkVGsLqAOqX8SbJ5AJ15wRPBHY2VkZEwoBYoOAo2x61ThF5AIgH7jm0TlWLVfWgszwXCBeRjl7Eb0xAeeKTDTz+iSWYQBQeGkKfztGszguuJLMESBWRFBFphacjP7NemUxgvLM8Fpjv9KlkAhnO3WcpQCqefpZG6xSRCcAo4DpVPTzvtYh0dvp5EJHhTuz7juegjfFX//hkPY9/sp6rhyZagglQfRJiAmpAZpPjZFS1WkTuBj4CQoHpqpotIg8BWaqaCUwDXhWRHDwtmAxn22wRmQ2sBqqBSapaA9BQnc4unwO2AgucnPKWcyfZWOAuEakGyoAMJ5EZE/BUlSc+3cA/PtnA1UMT+evYQZZgAlTfhBje+nYHRaVVxEaFux3OCZNA/pxOT0/XrKwst8Mw5oSoKo99vI6nP9toCSYIfLZ2D7e8vIR/3XkGpyV3cCUGEVmqqunNUZdNK2OMD1NVHv7PGp7+bCPXDU/ib5ZgAl6fzjEArNsVGJfMbFoZY3xUba1yf+YqXlu4jZvPTOaBy9Jssssg0DU2kuiIMDYESL+MJRljfFBNrTL5rZXMzsrljh/15L7R/SzBBAkRITUhmnUBkmTscpkxPqa6ppZ7Zi9ndlYuvzg/1RJMEOqbEMOG3QfdDqNZWJIxxodUVtfyP28u493lO/nNqL786sI+lmCCUGpCDPtKKtl7sMLtUE6YJRljfER5VQ13vbaUD1bt4g+XpjHpx73dDsm4pG+Cp/N/fQB0/luSMcYHlFXWcPsrWXy6dg9/uuIUbjs7xe2QjIv6JEQDBMSgTOv4N8ZlJRXV3DZjCYs2F/DXsYO4Nj2p6Y1MQIuPiaBdVDjrAqBfxpKMMS4qLK3k5peW8N2OIv4xbghjhtSf4NwEIxGhV3w0m/L9P8nY5TJjXLKnuJxxzy9k9c5inv3pUEsw5gg9O7ZhUwA8V8aSjDEu2F5QytjnFrB9fykv33IaPxnQ2e2QjI/pGR9N/oEKDpRXuR3KCbEkY8xJtn73Aa5+9huKy6t44/YRnNnbnlhhfqhnfBsANuX7d2vGkowxJ9GK7YVc+/wCAGZNPIMhSe1cjsj4ql6Hksxe/+6XsY5/Y06Sbzbu5fYZWXSIbsXrt42ge1yU2yEZH9a9QxtCQ4SNe/y7JWNJxpiTYN7q3Ux641uS46J49bbTSWgb6XZIxse1CgshqX1ra8kYY47u7WW5/PpfKzmla1tevmU47du0cjsk4yd6xkdbn4wxpnEvf72ZX81awfDkDrx++whLMOaY9OzYhs17S6it9d+HS1qSMaYFqCp//XAtD763mgvTEnjpltOIjrALB+bY9IyPpqK6lh2FZW6Hcty8SjIiMlpE1olIjojc18D7ESIyy3l/kYgk13lvsrN+nYiMaqpOEXndWb9KRKaLSLizXkTkSaf8ShEZeiIHbkxLqaqp5TdzVvLM5xu5bnh3nv3pUCLDQ90Oy/ihw7cx+/GgzCaTjIiEAk8DFwFpwHUiklav2G3AflXtDTwOPOpsmwZkAAOA0cAzIhLaRJ2vA/2AgUBrYIKz/iIg1fmZCDx7PAdsTEsqraxm4itZzFnqeRbMn688hbBQu2Bgjs/3Y2X8t/Pfm9/+4UCOqm5S1UpgJjCmXpkxwAxneQ5wvngegjEGmKmqFaq6Gchx6mu0TlWdqw5gMZBYZx+vOG8tBNqJSJfjPG5jml1BSSXXv7CIL9bn8/AVp9izYMwJi4+OICYijM2B3JIBugHb67zOddY1WEZVq4EiIO4o2zZZp3OZ7Ebgw2OIAxGZKCJZIpKVn5/vxeEZc+Jy95cy9rlvWJ1XzDM/HcYNI3q4HZIJACJC97gotu4rdTuU4+ZNkmnoq1j9Wx0aK3Os6+t6BvhSVb86hjhQ1amqmq6q6fHx8Q1sYkzzWpNXzFXPfMPeAxW8dtvpjD7F5iEzzadHXBTbCgI7yeQCdR9wkQjsbKyMiIQBsUDBUbY9ap0i8gAQD9xzjHEYc1It3LSPa59fQIgI/7rzTIandHA7JBNgundoQ+7+Umr89DZmb5LMEiBVRFJEpBWejvzMemUygfHO8lhgvtOnkglkOHefpeDptF98tDpFZAIwCrhOVWvr7eMm5y6zEUCRquYdxzEb0yzmfpfHTdMX0ykmgn//7Ez6do5xOyQTgHrERVFVo+z009uYm7xxX1WrReRu4CMgFJiuqtki8hCQpaqZwDTgVRHJwdOCyXC2zRaR2cBqoBqYpKo1AA3V6ezyOWArsMDpNH1LVR8C5gIX47l5oBS4pTlOgDHHSlV54atN/HnuWob1aM+LN6XbIEvTYnp08Mxxt62glKQO/jffnVejw1R1Lp4P+brr7q+zXA5c08i2jwCPeFOns77BmJyW0SRv4jWmpVTX1PLge9m8tnAblwzswt+vHWxjYEyLOjSR6tZ9pZzV2+VgjoMNQTbGSyUV1dz9xrd8ti6fO37Uk9+O6kdIiN2ibFpWl9jWhIcKWwv88zZmSzLGeGF3cTm3vryENXnFPHzFKXaLsjlpQkOEpPZRbPPT25gtyRjThLW7irn1pSUUllUxbfxp/LhfJ7dDMkHGn8fK2HwXxhzFfzfs5ZpnF1Bdq8y+4wxLMMYVPTp4xsp4uqb9iyUZYxoxO2s7N7+0mK7tWvPOpLM4pVus2yGZINU9rg0HK6opKKl0O5RjZpfLjKmntlb5+7x1PP3ZRs5J7cjTPx1K28hwt8MyQezQbcxbC0qJi45wOZpjYy0ZY+ooqajmrteX8vRnG8k4LYnpN59mCca4rodzG7M/dv5bS8YYx87CMibMyGLtrmL+cGkat56VbLMoG59waBCmP3b+W5IxBli2bT+3v7KU8qoapt18Gj/uax38xndEhoeS0DaC7fstyRjjd95dvoPfzFlJQtsI3rj9dPok2Bxkxvckto9ix37/m7/MkowJWrW1yj8+Wc+T83MYntyB524cRgebg8z4qG7tWrNs+363wzhm1vFvglJpZTWT3viWJ+fncM2wRF6bcLolGOPTEtu3Jq+w3O+m/LeWjAk6OwvLmPhqFtk7i/n9xf2ZcE6KdfAbn9etfWuqa5XdxeV0bdfa7XC8ZknGBJVFm/Yx6Y1vKa+q5cWb0jm/f4LbIRnjlcT2njvMdhSW+VWSsctlJiioKq8u2MJPX1xE28hw3pl0piUY41e6OYkl18/uMLOWjAl4FdU13P9ONrOytnNev078I2OIDbA0ficovkUsAAAX20lEQVSxvSfJ+NsdZpZkTEDbVVTOna8tZfn2Qv7nvN786oI+9gwY45ciw0PpGN2KXEsyxviGpVsLuPO1bymtqOa5G4Yy+pQubodkzAnp1j6KHYX+lWS86pMRkdEisk5EckTkvgbejxCRWc77i0Qkuc57k53160RkVFN1isjdzjoVkY511o8UkSIRWe78HH78szH1vbFoGxlTF9KmVShvTzrLEowJCIntWgdeS0ZEQoGngQuBXGCJiGSq6uo6xW4D9qtqbxHJAB4FxolIGpABDAC6Ap+ISB9nm8bq/Bp4H/i8gXC+UtVLj+M4TZCoqK7hj++t5o1F2xjZN54nxp1KbJT1v5jAkNi+NfNW76a2Vv3msq83l8uGAzmquglARGYCY4C6SWYM8KCzPAf4p3gGHowBZqpqBbBZRHKc+misTlVd5qw7keMyQWhHYRk/e20pK3KLuGtkL379k76E+skfojHe6Na+NZU1tew9WEGntpFuh+MVby6XdQO213md66xrsIyqVgNFQNxRtvWmzoacISIrROQDERnQUAERmSgiWSKSlZ+f70WVJhB8sT6fS5/8ik35JTx3wzB+O7qfJRgTcA7dYbbdjy6ZeZNkGvpLrT+vQWNljnX90XwL9FDVwcBTwDsNFVLVqaqarqrp8fHxTVRp/N2h+cdufmkxCW0jyfyfsxl9Sme3wzKmRXRr9/2ATH/hzeWyXCCpzutEYGcjZXJFJAyIBQqa2LapOo+gqsV1lueKyDMi0lFV93pxDCYAFZRU8stZy/lyfT5XDe3GI1cMpHWrULfDMqbFdGvvfwMyvWnJLAFSRSRFRFrh6cjPrFcmExjvLI8F5quqOusznLvPUoBUYLGXdR5BRDo7/TyIyHAn9n3eHKQJPMu3F3LZU/9l4cZ9/PnKgfz9msGWYEzAi44Io11UuF8NyGyyJaOq1SJyN/AREApMV9VsEXkIyFLVTGAa8KrTsV+AJ2nglJuN5yaBamCSqtaA51bl+nU6638O3At0BlaKyFxVnYAned0lItVAGZDhJDITRFSV1xZt46H3skloG8m/7zqTgYmxbodlzEnTrV1rv7pcJoH8OZ2enq5ZWVluh2GaSUlFNf/7zireXraDH/eN5/FxQ2gXZdPzm+AyYUYWuftL+fCX57bYPkRkqaqmN0ddNuLf+IU1ecVMeuNbNu8t4f9d2IdJP+7tN+MEjGlOXWIjWbzZf3oKLMkYn6aqvLF4G398bzXtWofz+oTTObNXx6Y3NCZAdWkXSXF5NSUV1bSJ8P2PcN+P0ASt4vIqJr/1Hf9Zmce5feKZcu1gOkZHuB2WMa7qEusZhJlXVE7vTtEuR9M0SzLGJ63YXsj/vLmMHYVl/HZ0P+44t6ddHjMG6BLruY05r6jMkowxx0pVmfbfzTz64Vo6xUQy+44RDOvRwe2wjPEZXQ8nmXKXI/GOJRnjM/aXVPKbOSv4ZM0eLkxL4G9jB9ndY8bUkxDruWScV2hJxhivLd5cwC9mLmPvwQoeuCyNm89MtklSjWlARJjn4WW7iv1jrIwlGeOqqppanvhkA898nkNShyj+fdeZDEps53ZYxvi0LrGt2WktGWOObsveEn4xazkrthdyzbBEHrh8ANF+cEumMW7rHBvJ1n0lbofhFfuLNiedqvKvpbk8mJlNWIjw9PVDuWSQPbnSGG91jY1k4Sb/GJBpScacVIWllfzu7e+Y+90uRvTswJRrh9C1XWu3wzLGr3SObc2B8moOVlT7fOvft6MzAeWbjXv5f7NXkH+ggvsu6sft5/S0B4sZcxy6tvMMyNxVVEbvTjEuR3N0lmRMi6usrmXKvPU8/+VGUuLa8PbPzrKZk405AYcGZO4sLLckY4Lb+t0HuGf2clbtKOa64d35w6X9iWplv3bGnIjvp5bx/duY7a/dtIiaWmXafzfx2MfriYkI4/kbhzFqgD0W2ZjmkNA2EhH/GPVvScY0u637Svj1v1awZMt+Rg1I4JErB9rElsY0o1ZhIXSMjvCLUf+WZEyzUVVeX7SNP89dQ2iIMOXawVx5ajcbuW9MC+gSG0lesSUZEyR2FZVz779X8uX6fM5J7cijVw+yW5ONaUFdYiPZlO/7AzJDvCkkIqNFZJ2I5IjIfQ28HyEis5z3F4lIcp33Jjvr14nIqKbqFJG7nXUqIh3rrBcRedJ5b6WIDD3egzbNR1V5e1kuP3n8C5ZsLuBPYwbwyq3DLcEY08K6xLYOjD4ZEQkFngYuBHKBJSKSqaqr6xS7Ddivqr1FJAN4FBgnImlABjAA6Ap8IiJ9nG0aq/Nr4H3g83qhXASkOj+nA886/xqX7DtYwe/fXsWH2bsY1qM9f79mMMkd27gdljFBIaFtJAcrfH9ApjeRDQdyVHUTgIjMBMYAdZPMGOBBZ3kO8E/xXIgfA8xU1Qpgs4jkOPXRWJ2qusxZVz+OMcArqqrAQhFpJyJdVDXvWA7YnDhV5b2VeTyYmc3B8mobWGmMCzo7U/7vLi4nOt53H17mTZLpBmyv8zqXH7YgDpdR1WoRKQLinPUL623bzVluqk5v4ugGHJFkRGQiMBGge/fuTVRpjtWe4nJ+/84q5q3ezeDEWP46djB9O/v2YDBjAlFCW89Ymd3F5fTy8yTT0NdT9bJMY+sb6guqX+fxxIGqTgWmAqSnpzdVp/GSqjJnaS5/en81FdW1/O7iftx6VgphoV516xljmlndJOPLvEkyuUBSndeJwM5GyuSKSBgQCxQ0sW1TdR5PHKYF7CgsY/Jb3/Hl+nyGJ3fg0bGDSLG+F2NcdSjJ7CqqcDmSo/Pma+gSIFVEUkSkFZ6O/Mx6ZTKB8c7yWGC+03eSCWQ4d5+l4Om0X+xlnfVlAjc5d5mNAIqsP6Zl1dYqry3cyk+mfEHWlgIeGjOAmRNHWIIxxgdER4QRHRHm/y0Zp4/lbuAjIBSYrqrZIvIQkKWqmcA04FWnY78AT9LAKTcbz00C1cAkVa0Bz63K9et01v8cuBfoDKwUkbmqOgGYC1wM5AClwC3NdRLMD23dV8Jv/72ShZsKOLt3R/7vqoEkdYhyOyxjTB0JbSN8PsmIp8ERmNLT0zUrK8vtMPxKdU0t07/ezJR56wkPCeF/L+3PtelJNmrfGB90/QsLKa+q4a2fndWs9YrIUlVNb466fPfmanPSrdheyOS3vmN1XjEX9E/g4StOobMz26sxxvd0bhvJos0FbodxVJZkDAcrqnnso3W8smAL8TERPHfDUEYN6GytF2N8XEJsJHsOlFNbq4T46Dg1SzJB7uPsXTyQmc2u4nJuHNGDX4/qS9vIcLfDMsZ4ISEmgqoapaC00mdnOrckE6Tyisp44N1sPl69m36dY3j6p0MZ2r2922EZY47BocvZu4rKLckY31BTq7y6YAuPfbye6tpafju6HxPOSSHcBlUa43c6OWNl9hwoxzM80fdYkgkiq3YU8ft3VrFieyHnpHbkkSsG0j3Obks2xl919oMBmZZkgkBRWRVTPl7Hqwu30qFNK57IGMLlg7tax74xfi4+JgIR355axpJMAFNV3vp2B//3wRoKSiq5cUQP7vlJX2JbW8e+MYEgPDSEuDa+PSDTkkyAWpNXzP3vrmLJlv2c2r0dL98ynFO6+eY1W2PM8esca0nGnEQHyqt4fN4GZizYQmzrcP569SDGDkv02XvojTEnJiEmkp0+/IRMSzIBQlXJXLGTh/+zhr0HK7h+eHd+M6ov7aJauR2aMaYFJcRGsmx7odthNMqSTABYu6uYBzOzWbipgEGJsbx4UzqDk9q5HZYx5iRIiImkoKSSiuoaIsJC3Q7nByzJ+LGCkkoen7ee1xdtJSYynEeuPIWM07rbY5CNCSKHHsO8p7jCJ2dKtyTjh6pqanlt4VYen7eeksoabhzRg19e0If2bezSmDHBpu6ATEsy5oR9sT6fP72/mpw9BzkntSN/uDSNPgkxbodljHGJrw/ItCTjJzbvLeHh91fz6do9JMdF8cJN6VzQv5MNqDQmyB1KMr56G7MlGR9XXF7FP+fn8NLXm4kIC2XyRf24+axkn+zgM8acfO2iwmkVFmJJxhybqppa3ly8jSc/3cC+kkquGZbIr0f1pVOMPUTMGPM9ESGhbQS7fDTJeDX1roiMFpF1IpIjIvc18H6EiMxy3l8kIsl13pvsrF8nIqOaqlNEUpw6Njh1tnLW3ywi+SKy3PmZcCIH7qtUlfdX7uTCKV9w/7vZ9IqPJnPS2fx17GBLMMaYBiXERPpvS0ZEQoGngQuBXGCJiGSq6uo6xW4D9qtqbxHJAB4FxolIGpABDAC6Ap+ISB9nm8bqfBR4XFVnishzTt3POtvMUtW7T/CYfdY3G/fylw/WsjK3iL4JMbx082mM7Btv/S7GmKNKiI1k9c5it8NokDeXy4YDOaq6CUBEZgJjgLpJZgzwoLM8B/ineD4ZxwAzVbUC2CwiOU59NFSniKwBzgOud8rMcOo9lGQC0pq8Yv7ywVq+WJ9P19hIHrtmMFee2s3GuxhjvJIQE8lnxXtQVZ/7UupNkukGbK/zOhc4vbEyqlotIkVAnLN+Yb1tuznLDdUZBxSqanUD5QGuFpFzgfXAr1S1bh1+J3d/KVM+Xs/by3fQNjKc313cj5vOSCYy3Dr1jTHe6xwbQWllDQcqqn3u8eneJJmG0qJ6Waax9Q31BR2tPMB7wJuqWiEid+Jp5Zz3g2BFJgITAbp3795Ade7bX1LJM5/nMOObrSAw8dye/OxHvYmN8q1fDmOMf0g4NCCzuNwvk0wukFTndSKws5EyuSIShuc5oAVNbNvQ+r1AOxEJc1ozh8ur6r465V/A03fzA6o6FZgKkJ6eXj8Zuqq8qoaXvt7CM5/ncLCimrFDE/nVhX3o2q6126EZY/xYQp0Bmb07+dbgbG+SzBIgVURSgB14OvKvr1cmExgPLADGAvNVVUUkE3hDRKbg6fhPBRbjabH8oE5nm8+cOmY6db4LICJdVDXP2d/lwJrjPOaTrqZW+ffSXKbMW8+u4nLO79eJe0f3o29n3/plMMb4p8Oj/n3wDrMmk4zTx3I38BEQCkxX1WwReQjIUtVMYBrwqtOxX4AnaeCUm43nJoFqYJKq1gA0VKezy98CM0XkYWCZUzfAz0XkcqeeAuDmEz76FqaqfLpmD49+uJYNew4yJKkd/8gYwoiecW6HZowJIAk+POpfVH3qilKzSk9P16ysLFf2vXTrfh79YC2LtxSQ0rEN947qy+hTOvvcnR/GmMAw5KGPuXRQFx6+YuAJ1yUiS1U1vRnCshH/zW1j/kH+9uE6PszeRcfoCB6+4hTGnZZEeKhX416NMea4dIltTV6h77VkLMk0k70HK3jikw28sXgbkWEh3HNhH247O4U2EXaKjTEtr0tsJHk++Bhm+wQ8QWWVNUz/ejPPfr6Rsqoarh/enV9ckErH6Ai3QzPGBJEusZEs27bf7TB+wJLMCfgoexd/zMxmZ1E5F6YlcN9F/egVH+12WMaYINQlNpL9pVWUV9X41IBuSzLHYWdhGfe/m80na3bTr3MMU8bZHWPGGHd1ifWMt8srKielYxuXo/meJZlj9FH2Lu6ds5LK6lomX9SPW89OsU59Y4zrusR6bmPOKyqzJOOPVJUp89bz1PwcBnaL5anrTiXZh/4jjTHBrYszc4iv3WFmScYLqsr/vrOK1xdtY1x6Eg9dMcCeTGmM8Sm+OurfkowXnpqfw+uLtnHHj3py3+h+NqDSGONzWrcKpX1UODsLy9wO5QjWmdCEhZv2MWXeeq46tZslGGOMT+sc25pdPjZWxpLMUVTV1DL5re/o3iGKh688xRKMMcandY2NZKclGf/xzrIdbN5bwv2XphHVyq4sGmN8W5d2keQV2eUyvzH96y3079KW8/t3cjsUY4xpUpfY1hSWVlFWWeN2KIdZkmlEzp4DrMkr5tr0RLtMZozxC3XHyvgKSzKN+HTNHgAuHtjF5UiMMcY7h56ym7vfkozPW7p1P8lxUYcfBmSMMb7u0Ej/LftKXI7ke5ZkGrEit5BTu7d3OwxjjPFap5gIolqFsnmvJRmfVl5Vw+7iCp+a/8cYY5oiIiTHtbEk4+sOXc9M6tDa5UiMMebYpHRswxZ/SzIiMlpE1olIjojc18D7ESIyy3l/kYgk13lvsrN+nYiMaqpOEUlx6tjg1NmqqX00t4KSSgDio60/xhjjX5I7RrF9fxlVNbVuhwJ4kWREJBR4GrgISAOuE5G0esVuA/aram/gceBRZ9s0IAMYAIwGnhGR0CbqfBR4XFVTgf1O3Y3uoyWUVlYDnrmAjDHGn/Tv0paaWiV7Z7HboQDetWSGAzmquklVK4GZwJh6ZcYAM5zlOcD54hlcMgaYqaoVqroZyHHqa7BOZ5vznDpw6ryiiX00u/Iqz0CmKEsyxhg/c3qK5wGKH67a5XIkHt7MldIN2F7ndS5wemNlVLVaRIqAOGf9wnrbdnOWG6ozDihU1eoGyje2j711AxGRicBE5+VBEdlXv4y30lqsreSajhznuQhAdi487Dx8L6DOxeRHYfLxbdoR6NFccXiTZBpqLaiXZRpb31AL6mjlvY0DVZ0KTD0cmEiWqqY3sG3QsXPxPTsXHnYevmfnwsM5D8nNVZ83l8tygaQ6rxOBnY2VEZEwIBYoOMq2ja3fC7Rz6qi/r8b2YYwxxkd5k2SWAKnOXV+t8HTkZ9YrkwmMd5bHAvNVVZ31Gc6dYSlAKrC4sTqdbT5z6sCp890m9mGMMcZHNXm5zOn/uBv4CAgFpqtqtog8BGSpaiYwDXhVRHLwtC4ynG2zRWQ2sBqoBiapag1AQ3U6u/wtMFNEHgaWOXXT2D68MLXpIkHDzsX37Fx42Hn4np0Lj2Y9D2KNAWOMMS3FRvwbY4xpMZZkjDHGtJiATjJNTYfj70RkuojsEZFVddZ1EJF5zrQ880SkvbNeRORJ51ysFJGhdbYZ75TfICLjG9qXrxORJBH5TETWiEi2iPzCWR9U50NEIkVksYiscM7DH531xzxdU2NTQvkbZ5aRZSLyvvM6KM+FiGwRke9EZLmIZDnrWv7vQ1UD8gfPDQUbgZ5AK2AFkOZ2XM18jOcCQ4FVddb9FbjPWb4PeNRZvhj4AM94oxHAImd9B2CT8297Z7m928d2HOeiCzDUWY4B1uOZsiiozodzPNHOcjiwyDm+2UCGs/454C5n+WfAc85yBjDLWU5z/mYigBTnbynU7eM7znNyD/AG8L7zOijPBbAF6FhvXYv/fQRyS8ab6XD8mqp+yQ/HCtWdfqf+tDyvqMdCPOORugCjgHmqWqCq+4F5eOaZ8yuqmqeq3zrLB4A1eGaJCKrz4RzPQedluPOjHPt0TY1NCeVXRCQRuAR40Xl9PFNXBcS5aESL/30EcpJpaDqcbo2UDSQJqpoHng9eoJOzvrHzEXDnybnMcSqeb/FBdz6cy0PLgT14PgQ24uV0TUDdKaH8+jw4/gHcCxyaktjrqasIvHOhwMcislQ802/BSfj78GZaGX/l1TQ0QeRYp/7xSyISDfwb+KWqFkvjc6gG7PlQz1i0ISLSDngb6N9QMeffgD0PInIpsEdVl4rIyEOrGyga8OfCcZaq7hSRTsA8EVl7lLLNdi4CuSXjzXQ4gWi306zF+XePs/5Yp/jxOyISjifBvK6qbzmrg/Z8qGoh8Dmea+rHOl1TIJyHs4DLRWQLnsvl5+Fp2QTjuUBVdzr/7sHz5WM4J+HvI5CTjDfT4QSiutPv1J+W5ybnrpERQJHTPP4I+ImItHfuLPmJs86vONfOpwFrVHVKnbeC6nyISLzTgkFEWgMX4OmfOtbpmhqbEspvqOpkVU1Uz2SPGXiO7acE4bkQkTYiEnNoGc/v9SpOxt+H23c8tOQPnjsk1uO5Jv17t+NpgeN7E8gDqvB8w7gNzzXkT4ENzr8dnLKC50FxG4HvgPQ69dyKpzMzB7jF7eM6znNxNp5m+0pgufNzcbCdD2AQnumYVjofIvc763vi+WDMAf4FRDjrI53XOc77PevU9Xvn/KwDLnL72E7wvIzk+7vLgu5cOMe8wvnJPvR5eDL+PmxaGWOMMS0mkC+XGWOMcZklGWOMMS3GkowxxpgWY0nGGGNMi7EkY4wxpsVYkjHGGNNiLMkYY4xpMf8fCbWrKjzpHYAAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.clf()\n",
"# plt.plot(x_part, calcs, '.')\n",
"plt.plot(test_q, calcs_test, label = 'pdf')\n",
"# plt.plot(test_q, res_y, label = 'res')\n",
"plt.legend()\n",
"plt.ylim(0.0, 4e-4)\n",
"# plt.yscale('log')\n",
"# plt.xlim(3080, 3110)\n",
"plt.savefig('test.png')\n",
"print(jpsi_width)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Adjust scaling of different parts"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"# total_f.update_integration_options(draws_per_dim=10000000, mc_sampler=None)\n",
"# inte = total_f.integrate(limits = (3000, 3200), norm_range=False)\n",
"# print(zfit.run(inte))\n",
"# print(pdg[\"jpsi_BR\"]/pdg[\"NR_BR\"], zfit.run(inte)/pdg[\"NR_auc\"])"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"# factor_jpsi = pdg[\"NR_auc\"]*pdg[\"jpsi_BR\"]/(pdg[\"NR_BR\"]*pdg[\"jpsi_auc\"])\n",
"# print(np.sqrt(factor_jpsi))\n",
"# factor_psi2s = pdg[\"NR_auc\"]*pdg[\"psi2s_BR\"]/(pdg[\"NR_BR\"]*pdg[\"psi2s_auc\"])\n",
"# print(np.sqrt(factor_psi2s))"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"# def _t_f(xq):\n",
"\n",
"# def jpsi_res(q):\n",
"# return resonance(q, jpsi_m, jpsi_s, jpsi_p, jpsi_w)\n",
"\n",
"# def psi2s_res(q):\n",
"# return resonance(q, psi2s_m, psi2s_s, psi2s_p, psi2s_w)\n",
"\n",
"# def cusp(q):\n",
"# return bifur_gauss(q, cusp_m, sig_L, sig_R, cusp_s)\n",
"\n",
"# funcs = psi2s_res(xq) + jpsi_res(xq) + cusp(xq)\n",
"\n",
"# vec_f = vec(xq, funcs)\n",
"\n",
"# axiv_nr = axiv_nonres(xq)\n",
"\n",
"# tot = vec_f + axiv_nr\n",
" \n",
"# return tot\n",
"\n",
"# def t_f(x):\n",
"# probs = zfit.run(_t_f(ztf.constant(x)))\n",
"# return probs"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"5404695.652173913\n"
]
}
],
"source": [
"print(36000*(1+ pdg[\"jpsi_BR\"]/pdg[\"NR_BR\"] + pdg[\"psi2s_BR\"]/pdg[\"NR_BR\"]))"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"# start = time.time()\n",
"\n",
"# result, err = integrate.quad(lambda x: t_f(x), x_min, x_max, limit = 50)\n",
"# print(result, \"{0:.2f} %\".format(err/result))\n",
"# print(\"Time:\", time.time()-start)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Sampling\n",
"## One sample"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"nevents = int(pdg[\"number_of_decays\"])\n",
"event_stack = 100000\n",
"\n",
"calls = int(nevents/event_stack + 1)\n",
"\n",
"# total_samp = []\n",
"\n",
"# start = time.time()\n",
"\n",
"# for call in range(calls):\n",
"# samp = total_f.sample(n=event_stack)\n",
"# sam = samp.unstack_x()\n",
"# sam = zfit.run(sam)\n",
"# clear_output(wait=True)\n",
" \n",
"# print(\"{0}/{1}\".format(call + 1, calls))\n",
"# print(\"Time taken: {}\".format(display_time(int(time.time() - start))))\n",
"# c = call + 1\n",
"# print(\"Projected time left: {}\".format(display_time(int((time.time() - start)/c*(calls-c)))))\n",
" \n",
" \n",
"# with open(\"data/zfit_toys/toy_1/{}.pkl\".format(call), \"wb\") as f:\n",
"# pkl.dump(sam, f, pkl.HIGHEST_PROTOCOL)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"# print(\"Time to generate full toy: {} s\".format(int(time.time()-start)))"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(5404696,)\n"
]
}
],
"source": [
"total_samp = []\n",
"\n",
"for call in range(calls):\n",
" with open(r\"data/zfit_toys/toy_1/{}.pkl\".format(call), \"rb\") as input_file:\n",
" sam = pkl.load(input_file)\n",
" total_samp = np.append(total_samp, sam)\n",
"\n",
"total_samp = total_samp.astype('float64')\n",
"\n",
"data2 = zfit.data.Data.from_numpy(array=total_samp[:int(nevents)], obs=obs)\n",
"\n",
"print(total_samp[:nevents].shape)"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAD8CAYAAACcjGjIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAGRpJREFUeJzt3X20ZXV93/H3JxBRTJUBwYUDFNRRfGgkeAOYtWJtiDyZBk3jCroqE0IcYyDGtKvN0HYFqyELWyLV1JA1hlGwBkTy4KwySkdSE1crD4MiDhrKFRAujDjJ4NgUg4Lf/nF+A4e7770znH3nnnNn3q+1zrr7fPdv7/uds86c7/k97H1TVUiSNOxHxp2AJGnyWBwkSR0WB0lSh8VBktRhcZAkdVgcJEkduywOSdYn+XaSLUOxg5NsSnJX+7mixZPkQ0mmk9ye5PihY1a39nclWT0Uf3WSr7ZjPpQki/2PlCQ9PbvTc/gYcNqs2FrghqpaBdzQngOcDqxqjzXAZTAoJsCFwInACcCFOwtKa7Nm6LjZv0uStMR2WRyq6q+B7bPCZwJXtO0rgDcOxa+sgRuBg5IcDpwKbKqq7VX1MLAJOK3te05VfbEGV+NdOXQuSdKY7D/icc+vqq0AVbU1yWEtvhK4f6jdTIstFJ+ZIz6nJGsY9DJ49rOf/epjjz12xPSXxlcf2PHE9j9Z+dwxZiJpkozrs+HWW2/926o6dHfajloc5jPXfEGNEJ9TVa0D1gFMTU3V5s2bR8lxyRy99rontjdf/IYxZiJpkozrsyHJN3e37airlR5qQ0K0n99u8RngyKF2RwAP7iJ+xBxxSdIYjVocNgA7VxytBj49FD+7rVo6CdjRhp+uB05JsqJNRJ8CXN/2/d8kJ7VVSmcPnUuSNCa7HFZKchXwOuB5SWYYrDq6GLgmybnAfcCbW/ONwBnANPAIcA5AVW1P8j7gltbuvVW1c5L7nQxWRD0L+Ex7SJLGaJfFoareMs+uk+doW8B585xnPbB+jvhm4JW7ykOStHS8QlqS1GFxkCR1WBwkSR0WB0lSh8VBktRhcZAkdVgcJEkdFgdJUofFQZLUYXGQJHVYHCRJHRYHSVKHxUGS1GFxkCR1WBwkSR0WB0lSh8VBktRhcZAkdfQqDkl+M8mWJHckeXeLvSfJA0lua48zhtpfkGQ6yZ1JTh2Kn9Zi00nW9slJktTfLv+G9HySvBJ4O3AC8H3gs0mua7svrapLZrV/OXAW8ArgBcDnkryk7f4w8HpgBrglyYaq+tqouUnScnH02uue2L734jeMMZOnGrk4AC8DbqyqRwCS/BXwpgXanwlcXVWPAvckmWZQWACmq+rudp6rW1uLgyQxngLSZ1hpC/DaJIckORA4Aziy7Ts/ye1J1idZ0WIrgfuHjp9psfnikqRFcPTa655SYHbHyMWhqr4OvB/YBHwW+ArwGHAZ8CLgOGAr8PvtkMx1mgXiHUnWJNmcZPO2bdtGTV2StAu9JqSr6vKqOr6qXgtsB+6qqoeq6vGq+iHwEZ4cOprhyZ4FwBHAgwvE5/p966pqqqqmDj300D6pS5IW0GfOgSSHVdW3kxwF/ALwmiSHV9XW1uRNDIafADYAf5LkAwwmpFcBNzPoOaxKcgzwAINJ67f2yUuS9kWLOTfRqzgAf5rkEOAHwHlV9XCSjyc5jsHQ0L3AOwCq6o4k1zCYaH6stX8cIMn5wPXAfsD6qrqjZ16StKw93TmCxdarOFTVT88Re9sC7S8CLpojvhHY2CcXSdLi8QppSVJH32ElSdIu7KkhooXO2/d32nOQJHXYc5Ckp2lPXbG8O9/2l2qi2uIgST1M6r2R+nJYSZLUYXGQJHU4rCRJi2RvGmKy5yBJ6rA4SJI6LA6SpA7nHCTt02ZfN7Dc5woWiz0HSVKHPQdJmsfetPro6bI4SNIeMO6/x9CXxUHSPmG5f1gvNeccJEkd9hwk7VXsISwOi4OkZWOcE8T7WtHpNayU5DeTbElyR5J3t9jBSTYluav9XNHiSfKhJNNJbk9y/NB5Vrf2dyVZ3e+fJEnqa+SeQ5JXAm8HTgC+D3w2yXUtdkNVXZxkLbAW+G3gdGBVe5wIXAacmORg4EJgCijg1iQbqurh0f9ZkpbafN+sl+Ib/r72rX4p9BlWehlwY1U9ApDkr4A3AWcCr2ttrgA+z6A4nAlcWVUF3JjkoCSHt7abqmp7O88m4DTgqh65SZpAuzMstC9fWzBJ+hSHLcBFSQ4BvgecAWwGnl9VWwGqamuSw1r7lcD9Q8fPtNh88Y4ka4A1AEcddVSP1CVpbvZCBkYuDlX19STvBzYBfw98BXhsgUMy12kWiM/1O9cB6wCmpqbmbCNp3+CH+J7Va7VSVV0OXA6Q5PcYfOt/KMnhrddwOPDt1nwGOHLo8COAB1v8dbPin++Tl6TJ54f7ZOu7Wumw9vMo4BcYzBNsAHauOFoNfLptbwDObquWTgJ2tOGn64FTkqxoK5tOaTFJ0pj0vc7hT9ucww+A86rq4SQXA9ckORe4D3hza7uRwbzENPAIcA5AVW1P8j7gltbuvTsnpyVJ49F3WOmn54j9HXDyHPECzpvnPOuB9X1ykSQtHq+QljQS5wz2bhYHSXtUnyJiARofi4Ok3eaH9b7DW3ZLkjosDpKkDoeVJAHe00hPZc9BktRhcZAkdTisJGlBrlDaN1kcpH2YH/yaj8VBUodFQxYHaYLsiRVDsz/oXYmk3WFxkPYSftvXYnK1kiSpw+IgSepwWElaxhxK0p5icZCWQN+JZm9toaVmcZD2MfY2tDt6FYckvwX8KlDAVxn8Xeg/Av4psKM1++Wqui1JgA8y+DvSj7T4l9p5VgP/obX/3aq6ok9e0t5gvg9xP9y1FEYuDklWAu8CXl5V30tyDXBW2/1vquraWYecDqxqjxOBy4ATkxwMXAhMMSgytybZUFUPj5qbJKmfvquV9geelWR/4EDgwQXanglcWQM3AgclORw4FdhUVdtbQdgEnNYzL0lSDyMXh6p6ALgEuA/YCuyoqv/Rdl+U5PYklyY5oMVWAvcPnWKmxeaLdyRZk2Rzks3btm0bNXVpXkevve6Jh7QvG7k4JFnBoDdwDPAC4NlJ/iVwAXAs8JPAwcBv7zxkjtPUAvFusGpdVU1V1dShhx46auqSpF3oM6z0s8A9VbWtqn4A/BnwU1W1tQ0dPQp8FDihtZ8Bjhw6/ggGw1DzxSVJY9JntdJ9wElJDgS+B5wMbE5yeFVtbauT3ghsae03AOcnuZrBhPSO1u564PdaTwTgFAa9D2lZmO8ahN0ZmnL4SpNq5OJQVTcluRb4EvAY8GVgHfCZJIcyGC66Dfi1dshGBstYpxksZT2nnWd7kvcBt7R2762q7aPmJUnqr9d1DlV1IYNlqMN+Zp62BZw3z771wPo+uUhLqc83fnsLWg688Z4kqcPiIEnqsDhIkjq88Z72akt9N1PnE7S3sOcgSeqwOEiSOhxWWiL+sZblzyEj7UssDton7U6xthhoX2Zx0LI0+4Pb3pi0uCwO2is4bCctLouD9joOB0n9uVpJktRhcZAkdTispH3GfMNNDkNJXfYcJEkd9hw0cVx5JI2fPQdJUoc9B43N0+0hODcgLZ1ePYckv5XkjiRbklyV5JlJjklyU5K7knwyyTNa2wPa8+m2/+ih81zQ4ncmObXfP0mS1NfIxSHJSuBdwFRVvRLYDzgLeD9waVWtAh4Gzm2HnAs8XFUvBi5t7Ujy8nbcK4DTgD9Mst+oeUmS+us7rLQ/8KwkPwAOBLYCPwO8te2/AngPcBlwZtsGuBb4r0nS4ldX1aPAPUmmgROAL/bMTUtoviEi74EkLU8j9xyq6gHgEuA+BkVhB3Ar8J2qeqw1mwFWtu2VwP3t2Mda+0OG43Mc8xRJ1iTZnGTztm3bRk1dkrQLI/cckqxg8K3/GOA7wKeA0+doWjsPmWfffPFusGodsA5gampqzjbas/ouM3VSWVoe+gwr/SxwT1VtA0jyZ8BPAQcl2b/1Do4AHmztZ4AjgZkk+wPPBbYPxXcaPkYTwA90ad/TpzjcB5yU5EDge8DJwGbgfwK/CFwNrAY+3dpvaM+/2Pb/ZVVVkg3AnyT5APACYBVwc4+8tAxZgKTJMnJxqKqbklwLfAl4DPgygyGf64Crk/xui13eDrkc+HibcN7OYIUSVXVHkmuAr7XznFdVj4+al8ZvMT/oLRrSePRarVRVFwIXzgrfzWC10ey2/wC8eZ7zXARc1CeXSeGHmaS9gVdIa04WOWnfZnHQEywIknayOCxD3rVU0p5mcdDI7GlIey+Lwxgs9Td/exqSni6Lw17EIiBpsVgc9lIO+Ujqw78EJ0nqsOewD7N3IWk+Focxc55A0iRyWEmS1GFxkCR1WBwkSR3OOSwT800eO6ksaU+w5yBJ6rDnsI+xpyFpd9hzkCR1WBwkSR0jDysleSnwyaHQC4HfAQ4C3g5sa/F/V1Ub2zEXAOcCjwPvqqrrW/w04IPAfsAfV9XFo+a1nHlBnKRJMXJxqKo7geMAkuwHPAD8OXAOcGlVXTLcPsnLgbOAVwAvAD6X5CVt94eB1wMzwC1JNlTV10bNTZLUz2JNSJ8MfKOqvplkvjZnAldX1aPAPUmmgRPavumquhsgydWt7diLw3yTt36rl7S3W6zicBZw1dDz85OcDWwG/nVVPQysBG4cajPTYgD3z4qfONcvSbIGWANw1FFHPe0kF+vD3hU/kvZ2vSekkzwD+HngUy10GfAiBkNOW4Hf39l0jsNrgXg3WLWuqqaqaurQQw/tlbckaX6L0XM4HfhSVT0EsPMnQJKPAP+9PZ0Bjhw67gjgwbY9X3yfZe9E0jgtRnF4C0NDSkkOr6qt7embgC1tewPwJ0k+wGBCehVwM4Oew6okxzCY1D4LeOsi5DUSP5QlqWdxSHIgg1VG7xgK/6ckxzEYGrp3576quiPJNQwmmh8Dzquqx9t5zgeuZ7CUdX1V3dEnL0lSP72KQ1U9AhwyK/a2BdpfBFw0R3wjsLFPLpKkxeO9lXAoSZJmW7a3z/jqAzs4eu11frBL0h6wV/ccLBySNJpl23OQJO05FgdJUofFQZLUYXGQJHVYHCRJHRYHSVKHxUGS1GFxkCR1WBwkSR17xRXSw1dC+yc8Jak/ew6SpA6LgySpY68YVhrmzfYkqT97DpKkDouDJKlj5OKQ5KVJbht6fDfJu5McnGRTkrvazxWtfZJ8KMl0ktuTHD90rtWt/V1JVi/GP0ySNLqRi0NV3VlVx1XVccCrgUeAPwfWAjdU1SrghvYc4HRgVXusAS4DSHIwcCFwInACcOHOgiJJGo/FGlY6GfhGVX0TOBO4osWvAN7Yts8ErqyBG4GDkhwOnApsqqrtVfUwsAk4bZHykiSNYLGKw1nAVW37+VW1FaD9PKzFVwL3Dx0z02LzxTuSrEmyOcnmxx/ZsUipS5Jm610ckjwD+HngU7tqOkesFoh3g1Xrqmqqqqb2O/C5Ty9RSdJuW4yew+nAl6rqofb8oTZcRPv57RafAY4cOu4I4MEF4pKkMVmM4vAWnhxSAtgA7FxxtBr49FD87LZq6SRgRxt2uh44JcmKNhF9SotJksak1xXSSQ4EXg+8Yyh8MXBNknOB+4A3t/hG4AxgmsHKpnMAqmp7kvcBt7R2762q7X3ykiT106s4VNUjwCGzYn/HYPXS7LYFnDfPedYD6/vkIklaPF4hLUnqsDhIkjosDpKkDouDJKnD4iBJ6rA4SJI6LA6SpA6LgySpw+IgSeqwOEiSOiwOkqQOi4MkqcPiIEnqsDhIkjosDpKkDouDJKnD4iBJ6rA4SJI6ehWHJAcluTbJ3yT5epLXJHlPkgeS3NYeZwy1vyDJdJI7k5w6FD+txaaTrO2TkySpv15/Qxr4IPDZqvrFJM8ADgROBS6tqkuGGyZ5OXAW8ArgBcDnkryk7f4w8HpgBrglyYaq+lrP3CRJIxq5OCR5DvBa4JcBqur7wPeTzHfImcDVVfUocE+SaeCEtm+6qu5u5726tbU4SNKY9BlWeiGwDfhoki8n+eMkz277zk9ye5L1SVa02Erg/qHjZ1psvrgkaUz6FIf9geOBy6rqJ4D/B6wFLgNeBBwHbAV+v7Wfq0tRC8Q7kqxJsjnJ5scf2dEjdUnSQvoUhxlgpqpuas+vBY6vqoeq6vGq+iHwEZ4cOpoBjhw6/gjgwQXiHVW1rqqmqmpqvwOf2yN1SdJCRi4OVfUt4P4kL22hk4GvJTl8qNmbgC1tewNwVpIDkhwDrAJuBm4BViU5pk1qn9XaSpLGpO9qpd8APtE+1O8GzgE+lOQ4BkND9wLvAKiqO5Jcw2Ci+THgvKp6HCDJ+cD1wH7A+qq6o2dekqQeehWHqroNmJoVftsC7S8CLpojvhHY2CcXSdLi8QppSVKHxUGS1GFxkCR1WBwkSR0WB0lSh8VBktRhcZAkdVgcJEkdFgdJUofFQZLUYXGQJHVYHCRJHRYHSVKHxUGS1GFxkCR1WBwkSR0WB0lSh8VBktRhcZAkdfQqDkkOSnJtkr9J8vUkr0lycJJNSe5qP1e0tknyoSTTSW5PcvzQeVa39nclWd33HyVJ6qdvz+GDwGer6ljgVcDXgbXADVW1CrihPQc4HVjVHmuAywCSHAxcCJwInABcuLOgSJLGY+TikOQ5wGuBywGq6vtV9R3gTOCK1uwK4I1t+0zgyhq4ETgoyeHAqcCmqtpeVQ8Dm4DTRs1LktTf/j2OfSGwDfhoklcBtwK/CTy/qrYCVNXWJIe19iuB+4eOn2mx+eIdSdYw6HUAPPrN9//clh75L7XnAX877iSepuWW83LLF8x5KSy3fGHP5fyPd7dhn+KwP3A88BtVdVOSD/LkENJcMkesFoh3g1XrgHUASTZX1dTTS3l8llu+sPxyXm75gjkvheWWL0xGzn3mHGaAmaq6qT2/lkGxeKgNF9F+fnuo/ZFDxx8BPLhAXJI0JiMXh6r6FnB/kpe20MnA14ANwM4VR6uBT7ftDcDZbdXSScCONvx0PXBKkhVtIvqUFpMkjUmfYSWA3wA+keQZwN3AOQwKzjVJzgXuA97c2m4EzgCmgUdaW6pqe5L3Abe0du+tqu278bvX9cx9qS23fGH55bzc8gVzXgrLLV+YgJxTNefwviRpH+YV0pKkDouDJKljYopDkmcmuTnJV5LckeQ/tvgxSW5qt9b4ZJvfIMkB7fl023/00LkuaPE7k5y6xPme3353JXneUPux3z5kgZw/0V6rLUnWJ/nRZZDz5S12e7uFy4+1+ES+L4b2/0GSvx96PtZ8F8o5yceS3JPktvY4rsUn+X2RJBcl+T8Z3NLnXZOQ8wL5fmHo9X0wyV9MQr4AVNVEPBhc7/BjbftHgZuAk4BrgLNa/I+Ad7btXwf+qG2fBXyybb8c+ApwAHAM8A1gvyXM9yeAo4F7gecNtT8D+Ew77iTgphY/mMFk/sHAira9Yolf4zPavgBXDb3Gk5zzc4bafABYO8nvi/Z8Cvg48PdD7cea7y5e448BvzhH+0l+X5wDXAn8SNt32CTkvND7YqjNnwJnT0K+VTU5PYca2PmN6kfbo4CfYXANBXRvx7HzNh3XAicnSYtfXVWPVtU9DFZHnbBU+VbVl6vq3jkOGfvtQxbIeWPbV8DNDK41mfScvwuDb1jAs3jywsmJfF8k2Q/4z8C/nXXIWPNdKOcFDpnY9wXwTgYrHn/Y2u28zmqsOe/qNU7yjxh81v3FJOQLEzSsBJBkvyS3MbhwbhODb0vfqarHWpPhW2s8cduNtn8HcAhP43Yci51vPXlB4Fx63z5kMSyUcwbDSW8DPrscck7yUeBbwLHAH8zOecLeF+cDG6rdWmbI2PNdIGeAi9qwxqVJDpid86zcJiHnFwG/lGRzks8kWTUpOe/i8+JNDG5Y+t1JyXeiikNVPV5VxzH45noC8LK5mrWfvW/H0dfsfJO8coHmY88XdpnzHwJ/XVVfaM8nOueqOgd4AYO7Af9Saz72nOfI97UMrvf5gzmajz1fmPc1voBB4f1JBsMYv92aT3LOBwD/UINbT3wEWN+ajz3nXfzfewuDId2dxp7vRBWHnWpwd9fPMxhrOyjJzov1hm+t8cRtN9r+5wLbGcPtOIbyXah7N1G3D5mdc5ILgUOBfzXUbKJzbrHHgU8C/6KFJvF98c+AFwPTSe4FDkwyPWn5zsr5tKra2oY1HgU+ypPDWpP8vphhMHYP8OfAj7fticl5jv97hzB4ba8bajb+fGsPTGSM8mDwwXRQ234W8AXg54BP8dQJ6V9v2+fx1Im8a9r2K3jqRN7d7JmJxznzHdp/L0+dkH4DT51gurmenGC6h8Hk0oq2ffASv8a/Cvxv4Fmz2k9qzv8ceHGLBbgEuGQ5vC9afHhCeqz57uJ9cfjQa/xfgIsn/H3xc8DFwK+0+OuAWyYh54XeF8CvAVfMaj/+13hPnHTEF+/HgS8DtwNbgN9p8RcymCSdZlAoDmjxZ7bn023/C4fO9e8ZzFfcCZy+xPm+i0F1f4xBRf/jFg/w4ZbXV4GpoXP9Svt3TAPnjOE1fqzldVt7/M4k58ygx/u/Wk5bgE/QVi9N6vtiVpvh4jDWfHfxvvjLodf4v/HkapuJfF+0+EEMvoF/Ffgi8KpJyHmh9wVP9tSG24/9Nfb2GZKkjomcc5AkjZfFQZLUYXGQJHVYHCRJHRYHSVKHxUGS1GFxkCR1/H/Hzgfuht0LLwAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"bins = int((x_max-x_min)/7)\n",
"\n",
"# calcs = zfit.run(total_test_tf(samp))\n",
"\n",
"plt.hist(total_samp, bins = bins, range = (x_min,x_max))\n",
"\n",
"# plt.plot(sam, calcs, '.')\n",
"# plt.plot(test_q, calcs_test)\n",
"plt.ylim(6000, 10000)\n",
"plt.xlim(3000, 3750)\n",
"\n",
"plt.savefig('test.png')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Toys"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"# sampler = total_f.create_sampler(n=nevents)\n",
"# nll = zfit.loss.UnbinnedNLL(model=total_f, data=sampler, fit_range = (x_min, x_max))\n",
"\n",
"# # for param in pdf.get_dependents():\n",
"# # param.set_value(initial_value)\n",
"\n",
"# sampler.resample(n=nevents)\n",
"\n",
"# # Randomise initial values\n",
"# # for param in pdf.get_dependents():\n",
"# # param.set_value(random value here)\n",
"\n",
"# # Minimise the NLL\n",
"# minimizer = zfit.minimize.MinuitMinimizer(verbosity = 10)\n",
"# minimum = minimizer.minimize(nll)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Fitting"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<hr>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
"<table>\n",
" <tr>\n",
" <td title=\"Minimum value of function\">FCN = 19224270.58100143</td>\n",
" <td title=\"Total number of call to FCN so far\">TOTAL NCALL = 58</td>\n",
" <td title=\"Number of call in last migrad\">NCALLS = 58</td>\n",
" </tr>\n",
" <tr>\n",
" <td title=\"Estimated distance to minimum\">EDM = 1.523018659051133e-07</td>\n",
" <td title=\"Maximum EDM definition of convergence\">GOAL EDM = 5e-06</td>\n",
" <td title=\"Error def. Amount of increase in FCN to be defined as 1 standard deviation\">\n",
" UP = 0.5</td>\n",
" </tr>\n",
"</table>\n",
"<table>\n",
" <tr>\n",
" <td align=\"center\" title=\"Validity of the migrad call\">Valid</td>\n",
" <td align=\"center\" title=\"Validity of parameters\">Valid Param</td>\n",
" <td align=\"center\" title=\"Is Covariance matrix accurate?\">Accurate Covar</td>\n",
" <td align=\"center\" title=\"Positive definiteness of covariance matrix\">PosDef</td>\n",
" <td align=\"center\" title=\"Was covariance matrix made posdef by adding diagonal element\">Made PosDef</td>\n",
" </tr>\n",
" <tr>\n",
" <td align=\"center\" style=\"background-color:#92CCA6\">True</td>\n",
" <td align=\"center\" style=\"background-color:#92CCA6\">True</td>\n",
" <td align=\"center\" style=\"background-color:#92CCA6\">True</td>\n",
" <td align=\"center\" style=\"background-color:#92CCA6\">True</td>\n",
" <td align=\"center\" style=\"background-color:#92CCA6\">False</td>\n",
" </tr>\n",
" <tr>\n",
" <td align=\"center\" title=\"Was last hesse call fail?\">Hesse Fail</td>\n",
" <td align=\"center\" title=\"Validity of covariance\">HasCov</td>\n",
" <td align=\"center\" title=\"Is EDM above goal EDM?\">Above EDM</td>\n",
" <td align=\"center\"></td>\n",
" <td align=\"center\" title=\"Did last migrad call reach max call limit?\">Reach calllim</td>\n",
" </tr>\n",
" <tr>\n",
" <td align=\"center\" style=\"background-color:#92CCA6\">False</td>\n",
" <td align=\"center\" style=\"background-color:#92CCA6\">True</td>\n",
" <td align=\"center\" style=\"background-color:#92CCA6\">False</td>\n",
" <td align=\"center\"></td>\n",
" <td align=\"center\" style=\"background-color:#92CCA6\">False</td>\n",
" </tr>\n",
"</table>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
"<table>\n",
" <tr>\n",
" <td><a href=\"#\" onclick=\"$('#SHoMzlGpgn').toggle()\">+</a></td>\n",
" <td title=\"Variable name\">Name</td>\n",
" <td title=\"Value of parameter\">Value</td>\n",
" <td title=\"Hesse error\">Hesse Error</td>\n",
" <td title=\"Minos lower error\">Minos Error-</td>\n",
" <td title=\"Minos upper error\">Minos Error+</td>\n",
" <td title=\"Lower limit of the parameter\">Limit-</td>\n",
" <td title=\"Upper limit of the parameter\">Limit+</td>\n",
" <td title=\"Is the parameter fixed in the fit\">Fixed?</td>\n",
" </tr>\n",
" <tr>\n",
" <td>0</td>\n",
" <td>psi2s_p</td>\n",
" <td>-9.42478</td>\n",
" <td>0.00651368</td>\n",
" <td></td>\n",
" <td></td>\n",
" <td></td>\n",
" <td></td>\n",
" <td>No</td>\n",
" </tr>\n",
" <tr>\n",
" <td>1</td>\n",
" <td>psi2s_s</td>\n",
" <td>74.7349</td>\n",
" <td>0.048852</td>\n",
" <td></td>\n",
" <td></td>\n",
" <td></td>\n",
" <td></td>\n",
" <td>No</td>\n",
" </tr>\n",
" <tr>\n",
" <td>2</td>\n",
" <td>jpsi_p</td>\n",
" <td>-9.42478</td>\n",
" <td>0.00631191</td>\n",
" <td></td>\n",
" <td></td>\n",
" <td></td>\n",
" <td></td>\n",
" <td>No</td>\n",
" </tr>\n",
" <tr>\n",
" <td>3</td>\n",
" <td>jpsi_s</td>\n",
" <td>444.495</td>\n",
" <td>0.220516</td>\n",
" <td></td>\n",
" <td></td>\n",
" <td></td>\n",
" <td></td>\n",
" <td>No</td>\n",
" </tr>\n",
"</table>\n",
"<pre id=\"SHoMzlGpgn\" style=\"display:none;\">\n",
"<textarea rows=\"14\" cols=\"50\" onclick=\"this.select()\" readonly>\n",
"\\begin{tabular}{|c|r|r|r|r|r|r|r|c|}\n",
"\\hline\n",
" & Name & Value & Hesse Error & Minos Error- & Minos Error+ & Limit- & Limit+ & Fixed?\\\\\n",
"\\hline\n",
"0 & $psi2s_{p}$ & -9.42478 & 0.00651368 & & & & & No\\\\\n",
"\\hline\n",
"1 & $psi2s_{s}$ & 74.7349 & 0.048852 & & & & & No\\\\\n",
"\\hline\n",
"2 & $jpsi_{p}$ & -9.42478 & 0.00631191 & & & & & No\\\\\n",
"\\hline\n",
"3 & $jpsi_{s}$ & 444.495 & 0.220516 & & & & & No\\\\\n",
"\\hline\n",
"\\end{tabular}\n",
"</textarea>\n",
"</pre>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
"<hr>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
"<span>Minos status for psi2s_p: <span style=\"background-color:#92CCA6\">VALID</span></span>\n",
"<table>\n",
" <tr>\n",
" <td title=\"lower and upper minos error of the parameter\">Error</td>\n",
" <td>-0.006515140889355781</td>\n",
" <td>0.006517797883262716</td>\n",
" </tr>\n",
" <tr>\n",
" <td title=\"Validity of minos error\">Valid</td>\n",
" <td style=\"background-color:#92CCA6\">True</td>\n",
" <td style=\"background-color:#92CCA6\">True</td>\n",
" </tr>\n",
" <tr>\n",
" <td title=\"Did minos error search hit limit of any parameter?\">At Limit</td>\n",
" <td style=\"background-color:#92CCA6\">False</td>\n",
" <td style=\"background-color:#92CCA6\">False</td>\n",
" </tr>\n",
" <tr>\n",
" <td title=\"I don't really know what this one means... Post it in issue if you know\">Max FCN</td>\n",
" <td style=\"background-color:#92CCA6\">False</td>\n",
" <td style=\"background-color:#92CCA6\">False</td>\n",
" </tr>\n",
" <tr>\n",
" <td title=\"New minimum found when doing minos scan.\">New Min</td>\n",
" <td style=\"background-color:#92CCA6\">False</td>\n",
" <td style=\"background-color:#92CCA6\">False</td>\n",
" </tr>\n",
"</table>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
"<span>Minos status for psi2s_s: <span style=\"background-color:#92CCA6\">VALID</span></span>\n",
"<table>\n",
" <tr>\n",
" <td title=\"lower and upper minos error of the parameter\">Error</td>\n",
" <td>-0.04922765388108916</td>\n",
" <td>0.0492610170079898</td>\n",
" </tr>\n",
" <tr>\n",
" <td title=\"Validity of minos error\">Valid</td>\n",
" <td style=\"background-color:#92CCA6\">True</td>\n",
" <td style=\"background-color:#92CCA6\">True</td>\n",
" </tr>\n",
" <tr>\n",
" <td title=\"Did minos error search hit limit of any parameter?\">At Limit</td>\n",
" <td style=\"background-color:#92CCA6\">False</td>\n",
" <td style=\"background-color:#92CCA6\">False</td>\n",
" </tr>\n",
" <tr>\n",
" <td title=\"I don't really know what this one means... Post it in issue if you know\">Max FCN</td>\n",
" <td style=\"background-color:#92CCA6\">False</td>\n",
" <td style=\"background-color:#92CCA6\">False</td>\n",
" </tr>\n",
" <tr>\n",
" <td title=\"New minimum found when doing minos scan.\">New Min</td>\n",
" <td style=\"background-color:#92CCA6\">False</td>\n",
" <td style=\"background-color:#92CCA6\">False</td>\n",
" </tr>\n",
"</table>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
"<span>Minos status for jpsi_p: <span style=\"background-color:#92CCA6\">VALID</span></span>\n",
"<table>\n",
" <tr>\n",
" <td title=\"lower and upper minos error of the parameter\">Error</td>\n",
" <td>-0.00631100168189077</td>\n",
" <td>0.006310920173057535</td>\n",
" </tr>\n",
" <tr>\n",
" <td title=\"Validity of minos error\">Valid</td>\n",
" <td style=\"background-color:#92CCA6\">True</td>\n",
" <td style=\"background-color:#92CCA6\">True</td>\n",
" </tr>\n",
" <tr>\n",
" <td title=\"Did minos error search hit limit of any parameter?\">At Limit</td>\n",
" <td style=\"background-color:#92CCA6\">False</td>\n",
" <td style=\"background-color:#92CCA6\">False</td>\n",
" </tr>\n",
" <tr>\n",
" <td title=\"I don't really know what this one means... Post it in issue if you know\">Max FCN</td>\n",
" <td style=\"background-color:#92CCA6\">False</td>\n",
" <td style=\"background-color:#92CCA6\">False</td>\n",
" </tr>\n",
" <tr>\n",
" <td title=\"New minimum found when doing minos scan.\">New Min</td>\n",
" <td style=\"background-color:#92CCA6\">False</td>\n",
" <td style=\"background-color:#92CCA6\">False</td>\n",
" </tr>\n",
"</table>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
"<span>Minos status for jpsi_s: <span style=\"background-color:#92CCA6\">VALID</span></span>\n",
"<table>\n",
" <tr>\n",
" <td title=\"lower and upper minos error of the parameter\">Error</td>\n",
" <td>-0.22065928968672088</td>\n",
" <td>0.22053487734368604</td>\n",
" </tr>\n",
" <tr>\n",
" <td title=\"Validity of minos error\">Valid</td>\n",
" <td style=\"background-color:#92CCA6\">True</td>\n",
" <td style=\"background-color:#92CCA6\">True</td>\n",
" </tr>\n",
" <tr>\n",
" <td title=\"Did minos error search hit limit of any parameter?\">At Limit</td>\n",
" <td style=\"background-color:#92CCA6\">False</td>\n",
" <td style=\"background-color:#92CCA6\">False</td>\n",
" </tr>\n",
" <tr>\n",
" <td title=\"I don't really know what this one means... Post it in issue if you know\">Max FCN</td>\n",
" <td style=\"background-color:#92CCA6\">False</td>\n",
" <td style=\"background-color:#92CCA6\">False</td>\n",
" </tr>\n",
" <tr>\n",
" <td title=\"New minimum found when doing minos scan.\">New Min</td>\n",
" <td style=\"background-color:#92CCA6\">False</td>\n",
" <td style=\"background-color:#92CCA6\">False</td>\n",
" </tr>\n",
"</table>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"psi2s_p: ^{+0.006517797883262716}_{-0.006515140889355781}\n",
"psi2s_s: ^{+0.0492610170079898}_{-0.04922765388108916}\n",
"jpsi_p: ^{+0.006310920173057535}_{-0.00631100168189077}\n",
"jpsi_s: ^{+0.22053487734368604}_{-0.22065928968672088}\n",
"Function minimum: 19224270.58100143\n"
]
}
],
"source": [
"nll = zfit.loss.UnbinnedNLL(model=total_f, data=data2, fit_range = (x_min, x_max))\n",
"\n",
"minimizer = zfit.minimize.MinuitMinimizer()\n",
"# minimizer._use_tfgrad = False\n",
"result = minimizer.minimize(nll)\n",
"\n",
"param_errors = result.error()\n",
"\n",
"for var, errors in param_errors.items():\n",
" print('{}: ^{{+{}}}_{{{}}}'.format(var.name, errors['upper'], errors['lower']))\n",
"\n",
"print(\"Function minimum:\", result.fmin)"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-3.1420346928204133"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"(-9.42522+2*np.pi)#/np.pi"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'5 h, 55 min'"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"display_time(int(395*pdg[\"number_of_decays\"]/100000))"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"6 h, 12 min\n"
]
}
],
"source": [
"print(display_time(22376))"
]
},
{
"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
}
