{
"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\n",
"import os"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# chunksize = 1000000\n",
"# zfit.run.chunking.active = True\n",
"# zfit.run.chunking.max_n_points = chunksize"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Build model and graphs\n",
"## Create graphs"
]
},
{
"cell_type": "code",
"execution_count": 3,
"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 = ztf.to_complex(scale*tf.abs(com))\n",
"\n",
" _phase = tf.angle(com)\n",
"\n",
" _phase += phase\n",
"\n",
" com = r * tf.exp(tf.complex(ztf.constant(0.0), _phase))\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": 4,
"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 ztf.to_complex(2) - G(ztf.to_complex(1) - 4*tf.pow(m, 2) / ztf.to_complex(tf.pow(q, 2)))\n",
"\n",
"def h_P(m,q):\n",
" \n",
" return ztf.to_complex(2/3) + (ztf.to_complex(1) - 4*tf.pow(m, 2) / ztf.to_complex(tf.pow(q, 2))) * h_S(m,q)\n",
"\n",
"def two_p_ccbar(mD, m_D_bar, m_D_star, q):\n",
" \n",
" \n",
" #Load constants\n",
" nu_D_bar = ztf.to_complex(pdg[\"nu_D_bar\"])\n",
" nu_D = ztf.to_complex(pdg[\"nu_D\"])\n",
" nu_D_star = ztf.to_complex(pdg[\"nu_D_star\"])\n",
" \n",
" phase_D_bar = ztf.to_complex(pdg[\"phase_D_bar\"])\n",
" phase_D = ztf.to_complex(pdg[\"phase_D\"])\n",
" phase_D_star = ztf.to_complex(pdg[\"phase_D_star\"])\n",
" \n",
" #Calculation\n",
" left_part = nu_D_bar * tf.exp(tf.complex(ztf.constant(0.0), phase_D_bar)) * h_S(m_D_bar, q) \n",
" \n",
" right_part_D = nu_D * tf.exp(tf.complex(ztf.constant(0.0), phase_D)) * h_P(m_D, q) \n",
" \n",
" right_part_D_star = nu_D_star * tf.exp(tf.complex(ztf.constant(0.0), phase_D_star)) * h_P(m_D_star, q) \n",
"\n",
" return left_part + right_part_D + right_part_D_star"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Build pdf"
]
},
{
"cell_type": "code",
"execution_count": 5,
"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": 6,
"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": 7,
"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",
"\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": 8,
"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": 9,
"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": 10,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAD8CAYAAACl69mTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl4VdW5+PHvS+YQSCAJJJAAAQIYBBHCJNaLVQuOtBYrjsggXivXtvZXh3tvbeutvdp6tdWqiKKiBYGitqkTdawTUxgljGEOUyCEhBAynOT9/XE2EGJCDpCwz/B+nieP+6yz9trv2ZLzZu219tqiqhhjjDEtoZXbARhjjAlelmSMMca0GEsyxhhjWowlGWOMMS3GkowxxpgWY0nGGGNMi/EpyYjIaBHZICL5IvJgA+9Hichc5/3FItKtznsPOeUbRGTUabT5jIiU+XIMY4wx/qnJJCMiYcCzwJVAFnCTiGTVqzYJKFbVnsBTwOPOvlnAOKAvMBp4TkTCmmpTRLKBBF+OYYwxxn/50pMZAuSr6hZVrQLmAGPq1RkDzHS25wOXiYg45XNUtVJVtwL5TnuNtukkoD8A9/t4DGOMMX4q3Ic6nYGddV4XAEMbq6OqHhEpARKd8kX19u3sbDfW5lQgR1X31MshjR3jQN1KIjIFmALQunXrQX369PHhIxpjjtl/uJK9pRUkxUWRGh/tdjgn8dQq6/aUAtCvc7zL0QSvZcuWHVDV5OZoy5ck01Bvof5aNI3Vaay8oR6Uikgn4AZg5BnGgapOB6YDZGdna25ubgO7GWMaM+1fm3ns/fVMuaQ7/3nVeW6Hc5L9hysZ/OhHAOQ+drXL0QQvEdneXG35crmsAEiv8zoN2N1YHREJB+KBg6fYt7HyC4GeQL6IbANiRSS/iWMYY5qRP1+D1m//XWn8nC9JZimQKSIZIhKJdyA/p16dHGC8sz0W+ES9K2/mAOOcmWEZQCawpLE2VfVdVU1R1W6q2g0odwb6T3UMY0yosN/4gNPk5TJn/GMqsAAIA15W1TwReQTIVdUcYAbwutPrOIg3aeDUmwesBTzAPapaA9BQm02E0uAxjDGhw3JM4PFlTAZVfQ94r17Zw3W2K/COpTS076PAo7602UCdOF+OcTqqq6spKCigoqLibJtyXXR0NGlpaURERLgdiglC/nihwA9DMk3wKckEk4KCAtq0aUO3bt0I5BnQqkpRUREFBQVkZGS4HY4x54SNyQSekFtWpqKigsTExIBOMAAiQmJiYlD0yIzxlfVkAk/IJRkg4BPMMcHyOYx/8ed/VpZjAk9IJhljTGDyx3Eic2qWZPzUZ599xjXXXANAZWUll19+OQMGDGDu3LkuR2ZChT9+n/tjTObUQm7gPxCtWLGC6upqVq5c6XYoJgSIH9+OaUkm8FhPxgXbtm2jT58+jB8/nv79+zN27FjKy8v54IMP6NOnDxdffDFvvfUWAIWFhdx6662sXLmSAQMGsHnzZpejN6HCH7/PPbW1bodgTlNI92R+84881u4ubdY2szq15VfX9m2y3oYNG5gxYwYjRoxg4sSJPPnkk7zwwgt88skn9OzZkxtvvBGADh068NJLL/HEE0/wzjvvNGusxgSamlp/TH3mVEI6ybgpPT2dESNGAHDrrbfy9NNPk5GRQWZm5vGy6dOnuxmiMX7H0wxJ5nBFNSt2HGLjvsNs3n+EkqNVlFXWEBXeirbREXRPbs15qW0YmpFI6yj7ijxbIX0GfelxtJT6049LSkpsSrLxC8f+Gfrj+MeZ9mQKSyt4e8UuPl5XyPIdxceTVfvWkbSLjSAuKpxKTy2Hyqt5c3kBAJFhrfhOZhK3De/KJZnJtGplv59nIqSTjJt27NjBwoULGT58OG+88QaXX345L7zwAps3b6ZHjx688cYbbodoQtSx5OKPf/OcTk9GVfkqv4gZX27hXxv3U6twfue2TLmkOxf1SOK81DYkxkV9a7/DFdWsLijhk/WF5KzazR2vLCUrtS2/vCaL4T0Sm/PjhARLMi4577zzmDlzJnfddReZmZn86U9/YtCgQVx99dUkJSVx8cUXs2bNGrfDNCHMD3MMNT4M/KsqC/L28dxn+awuKCG5TRR3j+zBDwem0T05rsn920RHMKJnEiN6JvHA6D68s3o3//fPjdz04iJ+lJ3Gr6/rS2ykfXX6ys6US1q1asW0adNOKhs9ejTr16//Vt2RI0cycuTIcxSZMf6r0nPqJLNy5yF++85acrcXk5HUmseu78cPBnYmKjzsjI4XGd6K6wemcVW/VJ7+eBPP/2szuduLeXn8YLoltT6jNkONJRljzEmOLULpj5fLyio8ANQfHik+UsX/vLuWt5bvIikuiseu78cN2emENdM4SnREGPeP7sN3MpO5Z/Zyrn/+a2aMz+bCLu2apf1gZvfJuKBbt252Kcz4rRNjMv6XZcoqvUmmdZ3LVe99s4crnvoXOSt3c/fIHnz2i5GMG9Kl2RJMXcN7JPLm3RcRFxXObTOW8E1BSbMfI9iEZJIJlvWPguVzGOOrA2WVAERFhFFe5eG+uSv58azlpMRHkzP1Yh4Y3Ye4Fp52nJHUmrl3DSM+JoLbX15MfuHhFj1eoAu5JBMdHU1RUVHAf0Efe55MdHS026GYIOPPvxkb95UB3mTz/We/4u2Vu/jJZZn87ccjyOrU9pzFkRofw+w7hxLWqhWTZ+ZScrT6nB070ITcmExaWhoFBQXs37/f7VDO2rEnYxrTnPz1769D5VV8tG7f8df7D1fy2sQhfCcz2ZV4uia25vlbB3LT9EX8bO5KXro92+6laYBPSUZERgN/AsKAl1T1sXrvRwGvAYOAIuBGVd3mvPcQMAmoAe5V1QWnalNEZgDZeGdQbgTuUNUyEbkD+AOwyznsn1X1pdP9wBEREfYkSWNO4fjAv8tx1HXwSBUTXl3K4QoPsyYPZfHWg9w4OJ3OCTGuxjW4W3sevjaLh/+ex6tfb2PixfbdUl+Tl8tEJAx4FrgSyAJuEpGsetUmAcWq2hN4Cnjc2TcLGAf0BUYDz4lIWBNt/kxVL1DV/sAOYGqd48xV1QHOz2knGGNM0/ytJ7PzYDljn/+a9XtKef6WgYzomcR9V/RyPcEcc9uwrlzWpwOPf7CezfvL3A7H7/gyJjMEyFfVLapaBcwBxtSrMwaY6WzPBy4T79SUMcAcVa1U1a1AvtNeo22qaimAs38M/n2J2BjTgtbuLuX657/mQFklsyYP5Xt9U9wO6VtEhP+9vh8xkWH8fN4qam0Rz5P4kmQ6AzvrvC5wyhqso6oeoARIPMW+p2xTRF4B9gJ9gGfq1PuhiKwWkfkiku5D7MaY03R8UozL18sWbi7ixhcWEt5KmH/3RWR3a+9uQKfQoW00D1+Txcqdh5i/rMDtcPyKL0mmoX9q9VN1Y3VOt9y7oToB6ASsA250iv8BdHMuo33EiZ7TyYGITBGRXBHJDYbBfWPOtdrjOca9LPPeN3sY//ISOsZH8+bdF9GrYxvXYvHVDy7szKCu7fj9gvWUVthss2N8STIFQN1eQxqwu7E6IhIOxAMHT7Fvk22qag0wF/ih87pIVSudt1/EO8ngW1R1uqpmq2p2crI7s06MCWRuj8m8vnAb98xeTr+0eOb/+3A6+cnYS1NEhF9f25eiI1X8+ZN8t8PxG74kmaVApohkiEgk3oH8nHp1coDxzvZY4BP19rlzgHEiEiUiGUAmsKSxNsWrJxwfk7kWWO+8Tq1zvOvw9nKMMc3MrWVlVJUn/7mBX/49j8v6dOAvk4aSEBt5boM4S/3S4rn+wjRmfr2NwtIKt8PxC00mGWeMZSqwAO8X+zxVzRORR0TkOqfaDCBRRPKB+4AHnX3zgHnAWuAD4B5VrWmsTbyX0WaKyDfAN0Aq8IhzjHtFJE9EVgH3Anec9ac3xnyLGz2ZmlrlP99ew9Of5POj7DSm3TqImMgzW9TSbT+5LBNPrfLcZ/aodPDxPhlVfQ94r17Zw3W2K4AbGtn3UeBRH9usBUY00s5DwEO+xGuMOXPHcsy5SjYV1TX8ZM4KFuTt455Le/D/vtfbL9dN81WXxFhuGJTG7MU7uOvfupMaHxiX+1pKyC0rY4xpwjnsypQcreb2l5ewIG8fv7o2i1+M6hPQCeaYey7tSY0qM77Y6nYorrMkY4w5yfGeTAvfolZYWsGNLyxkxY5i/jRuABNGBM/d8untY7m6Xypzlu7kcIjPNLMkY4w5Se2xnkwL5pidB8sZO20hOw6WM2P8YMYMqH/rXeCb/J0Myio9zF26s+nKQcySjDHmJC19tWzrgSPc+MJCDpVXMWvyUC7pFZy3GvRPS2BIt/a88tU2PDVNPzY6WFmSMcacROv9tzlt3HeYH72wkApPLW9MGRb0T5aceHE3dh06ymcbQvfGcEsyxpiTHL9a1sxdmjW7Shg3fREAc6cMo2+n+GZt3x9ddl5HkuKimJsbupfMLMkYY07SEgP+K3YUc/OLi4gOb8W8u4aTGQDLxDSHiLBWjB2UxifrC0P25kxLMsaYkx3vyTRPc4u3FHHrS4tJiI1k3r8PJyOpdfM0HCBuHJxOTa0yf3loLpxpScYYc5Jjs8uaI8d8uekA419ZQkp8NPPuGk5au9hmaDWwZCS1ZmhGe/6aWxDwj30/E5ZkjDEn8TjLMJ/t9+GXmw4waeZSuiW2Zu5dw0mJj26G6ALTDwemsfXAEdbsKnU7lHPOkowx5iQ1zfDQra/yvQkmI6k1s+8cRlJcVDNEFrhG9U0hIkz4x+r6C9gHP0syxpiTHO/JnOEFs6/rJJhZk4fSvnVgraTcEuJjI7gkM5l3Vu0OuSdnWpIxxpykpubML5d9nX+AiTOX0rW9N8EkhngPpq5rL+jE7pIKlu8odjuUc8qSjDHmJDVnOBizcHMRE2cupUv7WGbdaQmmvsuzOhIV3oqcVaF1ycySjDHmJGcyJrNoSxETX11KertYG4NpRFxUOCN7J/PPvH0hNcvMkowx5iTHxmRqffwiXLSliAmvLCWtXYwlmCZcfl5H9pZWkLc7dGaZWZIxxpykpta7mKPHhx7NYifBdHYSTHIbSzCncmmfDojAR+v2uR3KOWNJxhhzEo8z8F/tOfXKwct3FDPh1aV0Sohm9p1DLcH4ICkuioFd2lmSMcaErqqapnsyebtLuOPlJSS3ieKNO4fRoU3o3mh5ui4/ryNrdpWyp+So26GcEz4lGREZLSIbRCRfRB5s4P0oEZnrvL9YRLrVee8hp3yDiIxqqk0RmSEiq0RktYjMF5G4po5hjGk+R6tqgBPJpr78wsPcPmMJcVHhzJo8lA5tLcGcjsvP6wDAx+sKXY7k3GgyyYhIGPAscCWQBdwkIln1qk0CilW1J/AU8LizbxYwDugLjAaeE5GwJtr8mapeoKr9gR3A1FMdwxjTvI5We5NMQ5fLdhSVc8tLixERZt05LCTXIjtbPTvEkdYuhs83hsYzZnzpyQwB8lV1i6pWAXOAMfXqjAFmOtvzgctERJzyOapaqapbgXynvUbbVNVSAGf/GE6s09fYMYwxzehYT6a6Xk9m96Gj3PzSIio9tcyaPDTkVlNuLiLCdzKTWLilKCSemOlLkukM1H3iToFT1mAdVfUAJUDiKfY9ZZsi8gqwF+gDPNPEMU4iIlNEJFdEcvfvD42/FIxpTuXHk8yJMZn9hyu59aXFlJRX8/rEofROCY3nwbSUET2TOFzhYfWuErdDaXG+JJmGegv1RwQbq3O65d4N1QlAJ2AdcONpxIGqTlfVbFXNTk4OzmeHG9NSVJWiI5XAiTGZQ+VV3DZjMXtKKnh5wmD6pQX/Ey1b2ogeSYh4V6oOdr4kmQIgvc7rNKD+ugjH64hIOBAPHDzFvk22qao1wFzgh00cwxjTTMoqPVRUe5NLeZWHwxXVjH95CVv2H+HF27MZ3K29yxEGh3atIzm/U7wlGcdSIFNEMkQkEu9Afk69OjnAeGd7LPCJetdNyAHGOTPDMoBMYEljbYpXTzg+JnMtsL6JYxhjmsmOg+WAdwmULfuPMOnVXPJ2l/LcLQO5ODPJ5eiCy4ieSSzfUcyRSo/bobSoJpOMM/4xFViA9/LVPFXNE5FHROQ6p9oMIFFE8oH7gAedffOAecBa4APgHlWtaaxNvJfEZorIN8A3QCrwyKmOYYxpPmucMYKr+6VSXlXDkm0H+eO4AVye1dHlyILPdzKT8NQqi7cWuR1Kiwr3pZKqvge8V6/s4TrbFcANjez7KPCoj23WAiMaaafRYxhjmseCvH2kxkfz0FV9iIkMY2TvZEb27uB2WEFpUNd2RIQJi7ce5Lt9gjeJ+5RkjDHBb82uEj7dUMg9I3uSEBvJr6/r63ZIQS06Ioz+aQks3RrcQ8u2rIwxhorqGh54czUJMRHceUl3t8MJGYO7teebXSVUODfABiNLMsaEOFXlN/9YS97uUp644QLiYyLcDilkDO7WjuoaZcWOQ26H0mIsyRgT4l74fAtvLNnB3SN7cNl5wTs24I+yu7ZHBJZuC95LZpZkjAlh83J38tj767n2gk784nu93Q4n5MTHRtC7YxtLMsaY4PPh2n089NY3fCczif+74QJatbKlAN0wuFt7lm8vDtp1zCzJGBOClmw9yNTZyzm/U1um3TqIyHD7KnDL4Iz2HKmqYd2ew26H0iLsX5YxIWblzkNMfHUpnRNiePmOwbSOsjsZ3HRhegIAK3cWuxxJy7AkY0wIWbOrhNtnLKZd6whm3TmUxDh7ZLLb0trFkBQXycqdwbkisyUZY0LE+r2l3DpjMW2iI5g9eRip8TFuh2TwPl/mgrQEVhUE5zRmSzLGhID8wsPc8uJiosJbMfvOoaS3tyda+pML0hPYvL+M0opqt0NpdpZkjAlyWw8c4eYXvY9Mnn3nMLom2hMt/c0F6QmowjcFwXfJzJKMMUFsR1E5N7+4CE+tMvvOofRIjnM7JNOAC5wHwa3cGXyXzCzJGBOkCorLuenFRRytruEvk4bSq6M9MtlfJcRGkpHUmlWWZIwxgWBvSQU3v7iY0opq/jJpKFmd2rodkmnCBWnxQTn4b0nGmCBTeLiCm19cxMEjVbw2cQjnd453OyTjgwvSE9hXWsm+0gq3Q2lWlmSMCSL7D1dyy4uL2VtawSsTBnNhl3Zuh2R81LeT94+BtbtLXY6keVmSMSZIFB6u4KYXF1FQfJQZ4wczuFt7t0Myp+G8VO+YWd7u4Jph5lOSEZHRIrJBRPJF5MEG3o8SkbnO+4tFpFud9x5yyjeIyKim2hSRWU75GhF5WUQinPKRIlIiIiudn4cxxgBQWFrBTdMXsav4KK9MGMzwHoluh2ROU5voCLomxrJ2T4j1ZEQkDHgWuBLIAm4Skax61SYBxaraE3gKeNzZNwsYB/QFRgPPiUhYE23OAvoA/YAYYHKd43yhqgOcn0fO5AMbE2wKSysY9+Ii9pRU8OqEwQzrbgkmUGWltg3Jy2VDgHxV3aKqVcAcYEy9OmOAmc72fOAyERGnfI6qVqrqViDfaa/RNlX1PXUAS4C0s/uIxgSvfaUVjJu+iL0lFbw6YQhDLcEEtKzUtmwrKqes0uN2KM3GlyTTGdhZ53WBU9ZgHVX1ACVA4in2bbJN5zLZbcAHdYqHi8gqEXlfRPo2FKyITBGRXBHJ3b9/vw8fz5jAtM+5RLavtIKZE4cwJMPGYALdsanm64LokpkvSaahJxmpj3VOt7yu54DPVfUL5/VyoKuqXgA8A/ytoWBVdbqqZqtqdnJyckNVjAl4e0u8PZhjCcYG+YNDMM4w8yXJFADpdV6nAbsbqyMi4UA8cPAU+56yTRH5FZAM3HesTFVLVbXM2X4PiBCRJB/iNyao7Ck5yrjpC9l/uJLXJg0h2xJM0OjYNor2rSNDLsksBTJFJENEIvEO5OfUq5MDjHe2xwKfOGMqOcA4Z/ZZBpCJd5yl0TZFZDIwCrhJVY8/j1REUpxxHkRkiBN70Zl8aGMC1e5DRxk3fREHyqqYOXEIg7paggkmIkJWalvy9gTPNOYmH4mnqh4RmQosAMKAl1U1T0QeAXJVNQeYAbwuIvl4ezDjnH3zRGQesBbwAPeoag1AQ206h5wGbAcWOjnlLWcm2VjgbhHxAEeBcU4iMyYkHEswxUeqeG3SEAbajZZBKatTW179ahuemlrCwwL/VkYJ5u/p7Oxszc3NdTsMY87azoPl3PzSIg4dqea1SUPsTv4gNn9ZAf/vr6v4+Of/5tqq2SKyTFWzm6OtwE+TxgS5rQeOcOMLCykpr+b1yUMtwQS5Xh29iWXTvsMuR9I8LMkY48c27TvMj15YSIWnljlThjMgPcHtkEwL69nBm2Q27itzOZLm0eSYjDHGHXm7S7htxhLCWglzpwwj054HExJiI8NJbx/DRuvJGGNaysqdh7hp+iKiw1sx767hlmBCTK8ObdgUJD0ZSzLG+Jml2w5y60uLSYiNZO5dw8lIau12SOYcy+zYhi0HyqiuqW26sp+zJGOMH/kq/wC3z1hCh7ZRzLtrOOntY90OybigV8c4qmuUbQeOuB3KWbMkY4yf+HR9IRNeXUrXxFjmThlOSny02yEZl/RyLo8Gw+C/JRlj/MAHa/Yy5fVcenWM4407h5HcJsrtkIyLenaIo5UQFIP/NrvMGJf9bcUufv7XVVyQFs8rE4YQHxPhdkjGZdERYXRpH8umwsBPMtaTMcZFry3cxk/nrmRIt/a8NmmoJRhzXGbHNna5zBhzZlSVZz7exMN/z+OKrI68MmEwcVF2YcGc0LNDHNuLjuAJ8BlmlmSMOcdqa5XfvruO//twI9cP7MzztwwkOiLM7bCMn8lIak11jVJQfNTtUM6KJRljziFPTS33v7maGV9u5Y6LuvHE2AuCYqVd0/x6JHvvj9pyILAvmdm/bmPOkYrqGu6ZvZz5ywr42eW9+NW1WbRq1dBDYo2BjCTvGmZb9gf2vTJ2EdiYc6Cs0sNdr+fyVX4Rv7o2iwkjMtwOyfi59q0jSYiNYEuA35BpScaYFlZ8pIo7Xl3Kml0lPPmjC7h+YJrbIZkAkZHUmq0B3pOxy2XGtKB9pRXcOH0h6/aUMu3WQZZgzGnpnhRnYzLGmIblF5Zx/XNfs6v4KK9OGMwVWR3dDskEmO7JrdlXWklZpcftUM6YT0lGREaLyAYRyReRBxt4P0pE5jrvLxaRbnXee8gp3yAio5pqU0RmOeVrRORlEYlwykVEnnbqrxaRgWfzwY1pScu2FzN22tdUemqZe9dwLuqR5HZIJgB1d1bgDuSFMptMMiISBjwLXAlkATeJSFa9apOAYlXtCTwFPO7smwWMA/oCo4HnRCSsiTZnAX2AfkAMMNkpvxLIdH6mAM+fyQc2pqV9vG4ft7y0iISYCN66+yLO7xzvdkgmQGU405g37w/cS2a+9GSGAPmqukVVq4A5wJh6dcYAM53t+cBlIiJO+RxVrVTVrUC+016jbarqe+oAlgBpdY7xmvPWIiBBRFLP8HMb0yLmLt3BlNeX0btjG+bffRFdEm2pfnPmuiW2RgS2BnNPBugM7KzzusApa7COqnqAEiDxFPs22aZzmew24IPTiAMRmSIiuSKSu3//fh8+njFnT1V5+uNNPPDmN1zcM4nZdw4jKc5WUjZnJzoijE7xMQF9r4wvSaahu8XUxzqnW17Xc8DnqvrFacSBqk5X1WxVzU5OTm5gF2OaV02t8t9/W8OTzjIxL43PprWtQ2aaSffk1gE9w8yX34QCIL3O6zRgdyN1CkQkHIgHDjaxb6NtisivgGTgrtOMw5hzqqK6hnvfWME/1+7j7pE9uH9Ub7xXio1pHt0SW7Ny5yFUNSD/bfnSk1kKZIpIhohE4h3Iz6lXJwcY72yPBT5xxlRygHHO7LMMvIP2S07VpohMBkYBN6lqbb1j3O7MMhsGlKjqnjP4zMY0i5Lyam6bsZgP1+3jV9dm8cDoPgH5JWD8W9fEWA5XeDhUXu12KGekyZ6MqnpEZCqwAAgDXlbVPBF5BMhV1RxgBvC6iOTj7cGMc/bNE5F5wFrAA9yjqjUADbXpHHIasB1Y6PzCvqWqjwDvAVfhnTxQDkxojhNgzJnYUVTOHa8uoeDgUZ656UKu6d/J7ZBMkEpv7508suNgOe1aR7oczenz6cKxqr6H90u+btnDdbYrgBsa2fdR4FFf2nTKG4zJ6Rnd40u8xrSkFTuKmTwzF0+t8vqkIQztnuh2SCaIdXVmKG4/WM4F6QkuR3P6bHTSmNPw/jd7+OnclXRsG80rEwbTIznO7ZBMkOvi9GR2Hix3OZIzY0nGGB+oKi99sZXfvb+OAekJvHR7Nok2RdmcA7GR4STFRbG9KDCnMVuSMaYJnppafv2PPP6yaAdX9UvhyR8NsCdZmnOqa2Is24usJ2NM0DlS6WHq7OV8umE/d/1bdx4Y1cceNGbOua7tY1m0pcjtMM6IrcJsTCP2llRww7SFfL7pAL/7QT8euvI8SzDGFentY9lTWkGlp8btUE6b9WSMacDa3aVMmrmU0qPVzBifzcjeHdwOyYSwromxqMLOg0fp2SGwJptYT8aYehbk7WXstK8B+Ou/X2QJxrju2DTmQJxhZj0ZYxyqynOfbeYPCzYwID2B6bcNokPbaLfDMub4DZmBOMPMkowxeNcge/DN1fxt5W7GDOjE4z/sbzPIjN9IjosiNjKMHQePuh3KabMkY0Je4eEK7np9GSt2HOIXo3rz45E9bA0y41dEhC7tY9lx0HoyxgSUNbtKmPJaLsXl1Uy7dRCjz09xOyRjGpTePjYgH8NsA/8mZH2wZg83TFsIwPy7h1uCMX4tvV0suw4dxbuMY+CwnowJOarKs5/m88Q/N3oH+G8fRIc2NsBv/FvndjGUV9VQXF5N+wBajdmSjAkpRyo93D9/Ne9+s4fvD+jEYzbAbwJE54QYAHYVH7UkY4w/2l50hCmvLWNT4WEeurIPUy7pbgP8JmCktXOSzKFy+qXFuxyN7yzJmJDw2YZC7n1jBSLCzIlD+E5mstshGXNajvVkCooDaxqzJRkT1FSV5//lvcGyd8c2TL8tmy7O3dPGBJKE2AhiI8PYdciSjDF+oe74yzVaZifVAAAXXElEQVT9U/n92P7ERto/eROYRIS0djEB15PxaQqziIwWkQ0iki8iDzbwfpSIzHXeXywi3eq895BTvkFERjXVpohMdcpURJLqlI8UkRIRWen8HH/8szH1bS86wvXPfc37a/bw0JV9eOamCy3BmIDXOSGGXQGWZJr8rRORMOBZ4AqgAFgqIjmqurZOtUlAsar2FJFxwOPAjSKSBYwD+gKdgI9EpJezT2NtfgW8A3zWQDhfqOo1Z/A5TQix8RcTrDq3i2H5jkNuh3FafOnJDAHyVXWLqlYBc4Ax9eqMAWY62/OBy8Q7bWcMMEdVK1V1K5DvtNdom6q6QlW3neXnMiGotlZ5+uNNTHh1KZ0SYvjH1IstwZig0jkhlpKj1ZRVetwOxWe+JJnOwM46rwucsgbrqKoHKAEST7GvL202ZLiIrBKR90Wkb0MVRGSKiOSKSO7+/ft9aNIEg+IjVUx4dSlPfriRMRd04q0fX2QD/CbodG534l6ZQOHLReqGbiSov65BY3UaK28ouTW1VsJyoKuqlonIVcDfgMxvNaI6HZgOkJ2dHVjrL5gzsnLnIe6ZtZz9hyv57ffP55ahXez+FxOUTkxjLqd3ShuXo/GNLz2ZAiC9zus0YHdjdUQkHIgHDp5iX1/aPImqlqpqmbP9HhBRd2KACT2qyusLt3GD84Cx+XcP59ZhXS3BmKCVfvyGzMDpyfiSZJYCmSKSISKReAfyc+rVyQHGO9tjgU/Uu4pbDjDOmX2WgbfnscTHNk8iIinOOA8iMsSJvciXD2mCz5FKDz+du5Jf/j2Pi3sm8e69F9M/LcHtsIxpUUlxUUSGtQquy2Wq6hGRqcACIAx4WVXzROQRIFdVc4AZwOsiko+3BzPO2TdPROYBawEPcI+q1oB3qnL9Np3ye4H7gRRgtYi8p6qT8Savu0XEAxwFxmmgLUdqmkV+4WHu/styNu8v4/99rxc/HtmTVq2s92KCX6tWQqeEaAoCqCcjwfw9nZ2drbm5uW6HYZrRP1bt5oE3VxMTEcbTN13IiJ52xdSEllteWsSRyhr+ds+IFjuGiCxT1ezmaMvuTjMBoaK6hv95Zy2zFu9gUNd2/PnmC0mNj3E7LGPOuc4JMXyyPnBmzlqSMX4vv7CMqbOXs37vYaZc0p1fjOpNRJg9b8+Eps4JsRwoq6TSU0NUuP8/psKSjPFbqsr8ZQU8/Pc8YiLDeGXCYC7t3cHtsIxxVWq89wF7haWVpLf3/3vBLMkYv1RW6eGXf1vD2yt2Max7e/5444WkxNvTK41JTfD+Huw+dNSSjDFnIm93CVNnr2B70RF+enkm//HdTMJs9pgxwImezN7SCpcj8Y0lGeM3VJXXFm7n0XfX0a51BLPvHMaw7oluh2WMX0lxJrzsKbEkY4zPDpVX8cCbq1mQt49LeyfzxA0XkBgX5XZYxviduKhw2kSHsydA7pWxJGNc9/XmA9w3dxUHyir576vPY+KIDLu50phTSI2Ptp6MMU2p8tTyf//cwPQvtpCR2Jq3fnyRLQ1jjA9S42NsTMaYU8kvPMxP5qwkb3cpNw/twn9ffZ49udIYH6XGR5O3u9TtMHxiv9XmnFJV/rJ4B4++u5bYyHCm3zaI7/VNcTssYwJKSnw0B8oqqfLUEhnu3zcmW5Ix58yBskoemL+aj9cXckmvZJ4Y258Obe3eF2NOVydnhtm+0gq/v1fGkow5Jz7dUMgv/rqK0goPv7o2i/HDu9ngvjFn6NiNyXtKLMmYEFde5eHx99czc+F2+qS0YdbkYQHzRD9j/FWnhGNJxv+nMVuSMS1m2faD/HzeKrYVlTNxRAb3j+5NdIT/L+hnjL8LpBsyLcmYZlfpqeGpDzcx/fPNpMbH8Madwxjew+7cN6a5HLshc68lGRNq1uwq4efzVrFh32FuGpLOf12dRVyU/TMzprl5b8i0y2UmRFTX1PLcp5t55pNNtG8dySt3DObSPrYsvzEtJSU+JiAul/k0wVpERovIBhHJF5EHG3g/SkTmOu8vFpFudd57yCnfICKjmmpTRKY6ZSoiSXXKRUSedt5bLSIDz/RDm+a1ad9hrn/ua576aCNX90/lnz+7xBKMMS2sU4AsLdNkT0ZEwoBngSuAAmCpiOSo6to61SYBxaraU0TGAY8DN4pIFjAO6At0Aj4SkV7OPo21+RXwDvBZvVCuBDKdn6HA885/jUtqapUZX27hiX9uJC4qnOduGchV/VLdDsuYkBAoN2T6crlsCJCvqlsARGQOMAaom2TGAL92tucDfxYRccrnqGolsFVE8p32aKxNVV3hlNWPYwzwmqoqsEhEEkQkVVX3nM4HNs0jv/AwD7z5Dcu2F3NFVkd+94N+JLexVZONOVc6xceg6v83ZPqSZDoDO+u8LuDbPYjjdVTVIyIlQKJTvqjevp2d7aba9CWOzsBJSUZEpgBTALp06dJEk+Z0VdfUMv3zLfzpo03ERoXx5I8u4AcXdm7ojwJjTAtKqfPwskBPMg19e6iPdRorb6hvV7/NM4kDVZ0OTAfIzs5uqk1zGvJ2l3D//NXk7S7lqn4p/Oa68633YoxLjj0hc7efP1fGlyRTAKTXeZ0G7G6kToGIhAPxwMEm9m2qzTOJw7SASk8Nz3ycz7R/bSYhNpJptw5k9Pk29mKMm471ZPb5+ZL/vowWLQUyRSRDRCLxDuTn1KuTA4x3tscCnzhjJznAOGf2WQbeQfslPrZZXw5wuzPLbBhQYuMxLW/Z9mKufvpL/vxpPtcN6MRH911iCcYYP9AmOoK4qHC/n2HWZE/GGWOZCiwAwoCXVTVPRB4BclU1B5gBvO4M7B/EmzRw6s3DO0nAA9yjqjXgnapcv02n/F7gfiAFWC0i76nqZOA94CogHygHJjTXSTDfVl7l4YkFG3nl662kto3mlQmDubS3TUs2xp+kxEf7/V3/4u1wBKfs7GzNzc11O4yA8+WmA/zn29+w42A5tw7rwgOj+9AmOsLtsIwx9dw2YzFllR7e/vGIZm1XRJapanZztGV3/JvjDpRV8ui763h7xS66JcYyZ8owhnW3NceM8VcpbaP5Mv+A22GckiUZg6ry19wCfvf+Oo5Uerj3uz358aU9bcVkY/xcSnw0hYcr8dTUEh7mnzdkWpIJcfmFZfzn29+wZOtBBndrx+9+0I/Mjva8F2MCQUp8NDW1yoGyquOzzfyNJZkQVVFdw3Ofbeb5z/KJiQjjsev78aPsdHtapTEBJDX+xMPLLMkYv/H15gP899tr2HLgCN8f0In/ujrLbqo0JgCltPU+vMyfZ5hZkgkhRWWV/O699by5vICuibG8PmkI38lMdjssY8wZOtGTsSRjXFRTq8xesoMnFmzgSKWHey7twX98N9MG9o0JcAmxEUSFt/Lru/4tyQS5FTuK+eXf17BmVynDuyfyyJi+NrBvTJAQEVL8/LkylmSCVFFZJb//YANzc3fSsW0Uz9x0Idf0T7XVko0JMilt/fuuf0syQab+pbG7LunOf1yWSVyU/a82JhilxkezbEex22E0yr55gohdGjMm9KTEx7CvZC+1teqXtyBYkgkCRWWV/GHBBuYs9V4ae/qmC7nWLo0ZExJS46OpqqnlYHkVSXH+dyuCJZkAVuWp5bWF2/jTx5s4WlXDlEu6c69dGjMmpBx/QmZJhSUZ0zxUlU/WF/Lou+vYcuAI/9YrmV9ecx49O9ilMWNCTUrbE0nm/M7xLkfzbZZkAsymfYf5n3fX8fnG/XRPbs0rdwzm0j72nBdjQtXxGzL99F4ZSzIB4lB5FX/8aBOvL9pObGQYv7wmi9uHdyXCT1deNcacG4lxUYS3EvaWHHU7lAZZkvFznppaZi3ewVMfbaT0aDU3D+3CfVf0pn3rSLdDM8b4gbBWQse2/ntDpiUZP6WqfLqhkP99bz2bCsu4qEciD1+bRZ+Utm6HZozxM/78GGafrrWIyGgR2SAi+SLyYAPvR4nIXOf9xSLSrc57DznlG0RkVFNtikiG08Ymp81Ip/wOEdkvIiudn8ln88H92fIdxdw4fRETX82luqaWF24bxKzJQy3BGGMalBIfzd5AHZMRkTDgWeAKoABYKiI5qrq2TrVJQLGq9hSRccDjwI0ikgWMA/oCnYCPRKSXs09jbT4OPKWqc0RkmtP2884+c1V16ll+Zr+VX1jGHxasZ0HePpLiovif75/PuMHpNu5ijDmllLbRfLq+EFX1u/vjfLlcNgTIV9UtACIyBxgD1E0yY4BfO9vzgT+L95OOAeaoaiWwVUTynfZoqE0RWQd8F7jZqTPTafdYkglK+0or+ONHG5mXW0B0eCvuu6IXky7OoLXd72KM8UFqfDTlVTWUVniIj4lwO5yT+PIt1hnYWed1ATC0sTqq6hGREiDRKV9Ub9/OznZDbSYCh1TV00B9gB+KyCXARuBnqlq3jYBTcrSaF/61mZe/2kpNrXL78K5MvbQniX54Q5Uxxn/VvSEzEJNMQ30v9bFOY+UNXf85VX2AfwBvqGqliPw73l7Od78VrMgUYApAly5dGmjOfRXVNfxl0Xb+/Gk+h8qr+f6ATtx3RW+6JMa6HZoxJgDVfQxz7xT/uinblyRTAKTXeZ0G7G6kToGIhAPxwMEm9m2o/ACQICLhTm/meH1VLapT/0W8YzffoqrTgekA2dnZ9ZOhq2pqlb+t2MWTH25k16GjXNIrmftH9fbLu3SNMYEjJd5/H8PsS5JZCmSKSAawC+9A/s316uQA44GFwFjgE1VVEckBZovIk3gH/jOBJXh7LN9q09nnU6eNOU6bfwcQkVRV3eMc7zpg3Rl+5nNOVflsw34e/2A96/cepl/neH4/tj8jeia5HZoxJgh0aBOFCH45w6zJJOOMsUwFFgBhwMuqmicijwC5qpoDzABedwb2D+JNGjj15uGdJOAB7lHVGoCG2nQO+QAwR0R+C6xw2ga4V0Suc9o5CNxx1p/+HFixo5jH3l/P4q0H6ZoYy59vvpCrzk/1yyW5jTGBKSKsFUlxUX7ZkxFVv7qi1Kyys7M1NzfXlWNv3l/GEws28P6avSTFRfKTyzK5cXAXIsNtOrIxpvld9+cvaRcbycyJQ5qu3AQRWaaq2c0Qlt3x39wKSyv408ebmLN0J9HhrfjZ5b2Y/B2bjmyMaVkpbaPZXlTudhjfYt98zeRIpYcXPt/Ci59vwVNby23DujL1uz398vkOxpjgkxofzaItRU1XPMcsyZwlT00tc3N38tSHmzhQVsnV/VO5f1Rvuia2djs0Y0wISYmPobTCw5FKj19dOfGfSALQws1FPPz3NWwqLGNwt3a8ePsgLuzSzu2wjDEh6Ni9MntLK+iRHOdyNCdYkjkDB49U8ei763hzeQHp7WOYdusgRvXt6HdrBhljQkdH5wmZ+0osyQS0r/MP8NO5Kykur+KeS3sw9dJMYiLD3A7LGBPiTtz171/TmC3J+EhVefGLLfzv++vpntSaVycMIauTLb1vjPEPKXUul/kTSzI+UFUee389L3y+hav7p/KHsf2JjbRTZ4zxH9ERYbSLjWCPnz2G2b4pffDSF1t54fMt3DasK7+5rq/drW+M8Usp8TF+d9e/3X7ehOU7ivnf99dxVb8USzDGGL+WGh/td2MylmROoaZWeWD+alLaRvP4D/tbgjHG+LWObaPZ52djMpZkTuGd1bvZVFjGL6/Jok20fz0IyBhj6kuNj+ZAWRWVnhq3QznOkswpzPhyK706xjGqb4rboRhjTJOOzTArLK10OZITLMk0Ysv+MlYXlPCj7HS7TGaMCQj+eK+MJZlGfLyuEICr+6e6HIkxxvim7mOY/YUlmUYs215Ml/axpDqPNTXGGH/nj49htiTTiFUFhxjYJcHtMIwxxmdxUeHERYX71V3/lmQaUFFdw56SCrr70SJzxhjji5T4aHYV2+Uyv1bg/A9Kb2+XyowxgaVnchz5hWVuh3GcT0lGREaLyAYRyReRBxt4P0pE5jrvLxaRbnXee8gp3yAio5pqU0QynDY2OW1GNnWM5nbwSBUAyXHRLXUIY4xpEb1T2rCt6AhHq/zjXpkmk4yIhAHPAlcCWcBNIpJVr9okoFhVewJPAY87+2YB44C+wGjgOREJa6LNx4GnVDUTKHbabvQYLaG8ygNgS/gbYwLOealtqFXYVHjY7VAA33oyQ4B8Vd2iqlXAHGBMvTpjgJnO9nzgMvE+wWsMMEdVK1V1K5DvtNdgm84+33XawGnz+00co9lVVHv/Aoi1JGOMCTD907wTlr7YdMDlSLx8WYW5M7CzzusCYGhjdVTVIyIlQKJTvqjevp2d7YbaTAQOqaqngfqNHeOkMykiU4ApzssyESmqX8dXWS3WV3JNEmd4LoKQnQsvOw8nBNW5mPo4TD2zXZOArs0Vhy9JpqHegvpYp7HyhnpQp6rvaxyo6nRg+vHARHJVNbuBfUOOnYsT7Fx42Xk4wc6Fl3MeujVXe75cLisA0uu8TgN2N1ZHRMKBeODgKfZtrPwAkOC0Uf9YjR3DGGOMn/IlySwFMp1ZX5F4B/Jz6tXJAcY722OBT1RVnfJxzsywDCATWNJYm84+nzpt4LT59yaOYYwxxk81ebnMGf+YCiwAwoCXVTVPRB4BclU1B5gBvC4i+Xh7F+OcffNEZB6wFvAA96hqDUBDbTqHfACYIyK/BVY4bdPYMXwwvekqIcPOxQl2LrzsPJxg58KrWc+DWGfAGGNMS7E7/o0xxrQYSzLGGGNaTFAnmaaWwwl0IvKyiBSKyJo6Ze1F5ENnWZ4PRaSdUy4i8rRzLlaLyMA6+4x36m8SkfENHcvfiUi6iHwqIutEJE9EfuKUh9T5EJFoEVkiIquc8/Abp/y0l2tqbEmoQOOsMrJCRN5xXofkuRCRbSLyjYisFJFcp6zlfz9UNSh/8E4o2Ax0ByKBVUCW23E182e8BBgIrKlT9nvgQWf7QeBxZ/sq4H289xsNAxY75e2BLc5/2znb7dz+bGdwLlKBgc52G2Aj3iWLQup8OJ8nztmOABY7n28eMM4pnwbc7Wz/GJjmbI8D5jrbWc7vTBSQ4fwuhbn9+c7wnNwHzAbecV6H5LkAtgFJ9cpa/PcjmHsyviyHE9BU9XO+fa9Q3eV36i/L85p6LcJ7P1IqMAr4UFUPqmox8CHedeYCiqruUdXlzvZhYB3eVSJC6nw4n+fYErwRzo9y+ss1NbYkVEARkTTgauAl5/WZLF0VFOeiES3++xHMSaah5XA6N1I3mHRU1T3g/eIFOjjljZ2PoDtPzmWOC/H+FR9y58O5PLQSKMT7JbAZH5drAuouCRXQ58HxR+B+oNZ57fPSVQTfuVDgnyKyTLzLb8E5+P3wZVmZQOXTMjQh5HSX/glIIhIHvAn8VFVLpfE1VIP2fKj3XrQBIpIAvA2c11A1579Bex5E5BqgUFWXicjIY8UNVA36c+EYoaq7RaQD8KGIrD9F3WY7F8Hck/FlOZxgtM/p1uL8t9ApP90lfgKOiETgTTCzVPUtpzhkz4eqHgI+w3tN/XSXawqG8zACuE5EtuG9XP5dvD2bUDwXqOpu57+FeP/4GMI5+P0I5iTjy3I4waju8jv1l+W53Zk1MgwocbrHC4DviUg7Z2bJ95yygOJcO58BrFPVJ+u8FVLnQ0SSnR4MIhIDXI53fOp0l2tqbEmogKGqD6lqmnoXexyH97PdQgieCxFpLSJtjm3j/Xe9hnPx++H2jIeW/ME7Q2Ij3mvS/+V2PC3w+d4A9gDVeP/CmIT3GvLHwCbnv+2duoL3QXGbgW+A7DrtTMQ7mJkPTHD7c53hubgYb7d9NbDS+bkq1M4H0B/vckyrnS+Rh53y7ni/GPOBvwJRTnm08zrfeb97nbb+yzk/G4Ar3f5sZ3leRnJidlnInQvnM69yfvKOfR+ei98PW1bGGGNMiwnmy2XGGGNcZknGGGNMi7EkY4wxpsVYkjHGGNNiLMkYY4xpMZZkjDHGtBhLMsYYY1rM/we5gZknDF8RXwAAAABJRU5ErkJggg==\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": 11,
"metadata": {},
"outputs": [],
"source": [
"# total_f.update_integration_options(draws_per_dim=2000000, mc_sampler=None)\n",
"# inte = total_f.integrate(limits = (3090, 3102), norm_range=False)\n",
"# inte_fl = zfit.run(inte)\n",
"# print(inte_fl)\n",
"# print(pdg[\"jpsi_BR\"]/pdg[\"NR_BR\"], inte_fl/pdg[\"NR_auc\"])"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"# print(x_min)\n",
"# print(x_max)\n",
"# # total_f.update_integration_options(draws_per_dim=2000000, mc_sampler=None)\n",
"# total_f.update_integration_options(mc_sampler=lambda dim, num_results,\n",
"# dtype: tf.random_uniform(maxval=1., shape=(num_results, dim), dtype=dtype),\n",
"# draws_per_dim=1000000)\n",
"# # _ = []\n",
"\n",
"# # for i in range(10):\n",
"\n",
"# # inte = total_f.integrate(limits = (x_min, x_max))\n",
"# # inte_fl = zfit.run(inte)\n",
"# # print(inte_fl)\n",
"# # _.append(inte_fl)\n",
"\n",
"# # print(\"mean:\", np.mean(_))\n",
"\n",
"# _ = time.time()\n",
"\n",
"# inte = total_f.integrate(limits = (x_min, x_max))\n",
"# inte_fl = zfit.run(inte)\n",
"# print(inte_fl)\n",
"# print(\"Time taken: {}\".format(display_time(int(time.time() - _))))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Tensorflow scaling"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"# def scaling_func(x):\n",
"\n",
"# funcs = resonance(x, _mass = ztf.constant(jpsi_mass), scale = ztf.constant(jpsi_scale), phase = ztf.constant(jpsi_phase), width = ztf.constant(jpsi_width)) + resonance(x, _mass = ztf.constant(psi2s_mass), scale = ztf.constant(psi2s_scale), phase = ztf.constant(psi2s_phase), width = ztf.constant(psi2s_width))\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\n",
"\n",
"\n",
"# def s_func(x):\n",
" \n",
"# q = ztf.constant(x)\n",
" \n",
"# return zfit.run(scaling_func(q))\n",
" \n",
"\n",
"# print(integrate.quad(s_func, x_min, x_max, limit = 50))"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"# factor_jpsi = pdg[\"NR_auc\"]*pdg[\"jpsi_BR\"]/(pdg[\"NR_BR\"]*pdg[\"jpsi_auc\"])\n",
"# factor_jpsi = pdg[\"NR_auc\"]*pdg[\"jpsi_BR\"]/(pdg[\"NR_BR\"]*inte_fl)\n",
"# print(np.sqrt(factor_jpsi)*jpsi_scale)\n",
"# print(np.sqrt(factor_jpsi))\n",
"# # print(psi2s_scale)\n",
"# factor_psi2s = pdg[\"NR_auc\"]*pdg[\"psi2s_BR\"]/(pdg[\"NR_BR\"]*pdg[\"psi2s_auc\"])\n",
"# factor_psi2s = pdg[\"NR_auc\"]*pdg[\"psi2s_BR\"]/(pdg[\"NR_BR\"]*inte_fl)\n",
"# print(np.sqrt(factor_psi2s)*psi2s_scale)\n",
"# print(np.sqrt(factor_psi2s))"
]
},
{
"cell_type": "code",
"execution_count": 15,
"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",
"# funcs = psi2s_res(xq) + jpsi_res(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",
"# _ = np.array(x)\n",
"# probs = zfit.run(_t_f(_))\n",
"# return probs"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"# print(36000*(1+ pdg[\"jpsi_BR\"]/pdg[\"NR_BR\"] + pdg[\"psi2s_BR\"]/pdg[\"NR_BR\"]))"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"# start = time.time()\n",
"\n",
"# result, err = integrate.quad(lambda x: t_f(x), x_min, x_max, limit = 5)\n",
"# print(result, \"{0:.2f} %\".format(err/result))\n",
"# print(\"Time:\", time.time()-start)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Sampling\n",
"## One sample\n",
"! total_f.sample() always returns the same set !"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"# nevents = int(pdg[\"number_of_decays\"])\n",
"# event_stack = 5000\n",
"\n",
"# calls = int(nevents/event_stack + 1)\n",
"\n",
"# total_samp = []\n",
"\n",
"# start = time.time()\n",
"\n",
"# samp = total_f.sample(n=event_stack)\n",
"# s = samp.unstack_x()\n",
"\n",
"# for call in range(calls):\n",
"\n",
"# sam = zfit.run(s)\n",
"# clear_output(wait=True)\n",
" \n",
"# # if call != 0:\n",
"# # print(np.sum(_last_sam-sam))\n",
" \n",
"# # _last_sam = sam\n",
" \n",
"# c = call + 1 \n",
"# print(\"{0}/{1}\".format(c, calls))\n",
"# print(\"Time taken: {}\".format(display_time(int(time.time() - start))))\n",
"# print(\"Projected time left: {}\".format(display_time(int((time.time() - start)/c*(calls-c)))))\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": 19,
"metadata": {},
"outputs": [],
"source": [
"# print(\"Time to generate full toy: {} s\".format(int(time.time()-start)))\n",
"\n",
"# 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": 20,
"metadata": {},
"outputs": [],
"source": [
"# bins = int((x_max-x_min)/7)\n",
"\n",
"# # calcs = zfit.run(total_test_tf(samp))\n",
"\n",
"# plt.hist(total_samp[:event_stack], bins = bins, range = (x_min,x_max))\n",
"\n",
"# # plt.plot(sam, calcs, '.')\n",
"# # plt.plot(test_q, calcs_test)\n",
"# plt.ylim(0, 20)\n",
"# # plt.xlim(3000, 3750)\n",
"\n",
"# plt.savefig('test2.png')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Toys"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"nr_of_toys = 1\n",
"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",
"# sampler = total_f.create_sampler(n=event_stack)\n",
"\n",
"# for toy in range(nr_of_toys):\n",
" \n",
"# dirName = 'data/zfit_toys/toy_{0}'.format(toy)\n",
" \n",
"# if not os.path.exists(dirName):\n",
"# os.mkdir(dirName)\n",
"# print(\"Directory \" , dirName , \" Created \")\n",
"\n",
"# for call in range(calls):\n",
"\n",
"# sampler.resample(n=event_stack)\n",
"# s = sampler.unstack_x()\n",
"# sam = zfit.run(s)\n",
"# # clear_output(wait=True)\n",
"\n",
"# c = call + 1 \n",
"# print(\"{0}/{1}\".format(c, calls))\n",
"# print(\"Time taken: {}\".format(display_time(int(time.time() - start))))\n",
"# print(\"Projected time left: {}\".format(display_time(int((time.time() - start)/c*(calls-c)))))\n",
"\n",
"# with open(\"data/zfit_toys/toy_{0}/{1}.pkl\".format(toy, call), \"wb\") as f:\n",
"# pkl.dump(sam, f, pkl.HIGHEST_PROTOCOL)"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"# with open(r\"data/zfit_toys/toy_0/0.pkl\", \"rb\") as input_file:\n",
"# sam = pkl.load(input_file)\n",
"# print(sam[:10])\n",
"\n",
"# with open(r\"data/zfit_toys/toy_0/1.pkl\", \"rb\") as input_file:\n",
"# sam2 = pkl.load(input_file)\n",
"# print(sam2[:10])\n",
"\n",
"# print(np.sum(sam-sam2))"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Time to generate full toy: 0 s\n",
"(5404696,)\n"
]
}
],
"source": [
"print(\"Time to generate full toy: {} s\".format(int(time.time()-start)))\n",
"\n",
"total_samp = []\n",
"\n",
"for call in range(calls):\n",
" with open(r\"data/zfit_toys/toy_0/{}.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",
"data3 = zfit.data.Data.from_numpy(array=total_samp, obs=obs)\n",
"\n",
"print(total_samp[:nevents].shape)"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAD8CAYAAACPWyg8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XeYU2X2wPHvSTIZBEGKgAjooKKAii5iWVEXRRHLgl10dRFR1t/aXctQFBFRdC2ra0VRcUWKiICOjUVZVBABKYL0PoIOMlSBSbnn90fuQMAZmH5Tzud55snNmzfJufeZ5OQt972iqhhjjDGl4fM6AGOMMcnHkocxxphSs+RhjDGm1Cx5GGOMKTVLHsYYY0rNkocxxphS22/yEJE3RCRPRObFlf1TRBaKyFwR+UBEasc91ktElorIIhE5P668k1u2VESy48qbicg0EVkiIiNFJFiRO2iMMabilaTl8RbQaa+yCcBxqtoaWAz0AhCRVkBX4Fj3OS+JiF9E/MCLwAVAK+Aaty7AE8Czqtoc2Aj0KNceGWOMqXT7TR6qOhnI36vsc1WNuHe/BZq4212AEapaoKorgKXAKe7fUlVdrqohYATQRUQEOAcY7T5/KHBJOffJGGNMJQtUwGvcCIx0txsTSyaFct0ygDV7lZ8K1AM2xSWi+Pq/IyI9gZ4ANWrUOKlFixblDt6YtLJ2Vuz20D94G0dR8n6ESAE0aAWBTK+jSVkzZ878VVXrl/d1ypU8RKQPEAGGFRYVUU0puoWj+6hfJFUdDAwGaNu2rc6YMaNU8RqT9h4+yL1NwM/OcyfAxpVw2xg4uLnX0aQsEVlVEa9T5uQhIt2Ai4EOunuBrFygaVy1JsBad7uo8l+B2iIScFsf8fWNMWmlqN+SJlGVaaquiHQCHgA6q+r2uIfGA11FJFNEmgHNge+A6UBzd2ZVkNig+ng36XwJXOE+vxswrmy7YoxJCbZYa1IoyVTd4cBU4BgRyRWRHsALQE1ggojMFpFXAFR1PjAK+BH4FLhVVaNuq+I24DNgATDKrQuxJHSPiCwlNgYypEL30BhjTIXbb7eVql5TRHGxX/CqOhAYWET5x8DHRZQvJzYbyxhj2MewZ4UIh8Pk5uayc+fOSn0fr1WrVo0mTZqQkZFRKa9fEbOtjDEmaeTm5lKzZk2ysrKInS2QelSVDRs2kJubS7NmzSrlPWx5EmNMYqnkMY+dO3dSr169lE0cACJCvXr1KrV1ZcnDGJMYqvDLPJUTR6HK3kdLHsaYBGOzrZKBJQ9jjElAkyZN4uKLLwagoKCAc889lxNPPJGRI0fu55lVwwbMjTGJxc7z+J1Zs2YRDoeZPXu216HsYi0PY0yCSP1xiEIrV66kRYsWdOvWjdatW3PFFVewfft2Pv30U1q0aMEZZ5zBmDFjAMjLy+O6665j9uzZnHjiiSxbtszj6GOs5WGMSTBV2PL4JBt+/qFiX/OQ4+GCQfuttmjRIoYMGUK7du248cYbeeaZZ3j11Vf54osvOOqoo7j66qsBaNCgAa+//jpPPfUUH330UcXGWg7W8jDGGA80bdqUdu3aAXDdddcxY8YMmjVrRvPmzRERrrvuOo8j3DdreRhjEktVjnmUoIVQWfaeSrt58+akmkJsLQ9jTGJIoi/OirB69WqmTp0KwPDhwzn33HNZsWLFrjGN4cOHexneflnyMMYkmPSYbdWyZUuGDh1K69atyc/P5+6772bw4MFcdNFFnHHGGRx++OFeh7hP1m1ljDEe8Pl8vPLKK3uUderUiYULF/6ubvv27Wnfvn0VRVYy1vIwxiQWO88jKVjyMMaYKpaVlcW8efO8DqNcLHkYYxKMtTySgSUPY0yCSK/ZVsnOkocxJrHYmEdSsORhjDGm1GyqrjEmwVRtyyMrO6dCX2/loItK/ZyHH36YAw88kHvvvbfIx8eOHcvRRx9Nq1atyhtehbGWhzEmMaTZGealMXbsWH788Uevw9iDJQ9jTGJJkzGPgQMHcswxx3DuueeyaNEiAF577TVOPvlkTjjhBC6//HK2b9/OlClTGD9+PPfdd9+uJdmLqlfVLHkYY0wVmzlzJiNGjGDWrFmMGTOG6dOnA3DZZZcxffp05syZQ8uWLRkyZAinn346nTt35p///CezZ8/myCOPLLJeVbMxD2NMgkn9lsdXX33FpZdeSvXq1QHo3LkzAPPmzaNv375s2rSJbdu2cf755xf5/JLWq0yWPIwxCSK9xjyKWn79hhtuYOzYsZxwwgm89dZbTJo0qcjnlrReZbJuK2NMYkn9hgdnnXUWH3zwATt27GDr1q18+OGHAGzdupVGjRoRDocZNmzYrvo1a9Zk69atu+4XV68qWcvDGJPWyjK1trzatGnD1VdfzYknnsjhhx/OmWeeCcCAAQM49dRTOfzwwzn++ON3JYyuXbty88038/zzzzN69Ohi61Ul0f3MbBCRN4CLgTxVPc4tqwuMBLKAlcBVqrpRYu2w54ALge3ADar6vfucbkBf92UfVdWhbvlJwFvAAcDHwJ26v6CAtm3b6owZM0qzr8aYhw9ybzd7G0dRXjgFfl0EN02EJm0r7W0WLFhAy5YtK+31E0lR+yoiM1W13Ae4JN1WbwGd9irLBiaqanNgonsf4AKgufvXE3jZDbYu0A84FTgF6CciddznvOzWLXze3u9ljDEmwew3eajqZCB/r+IuwFB3eyhwSVz52xrzLVBbRBoB5wMTVDVfVTcCE4BO7mO1VHWq29p4O+61jDHpKE3O80h2ZR0wb6iq6wDc2wZueWNgTVy9XLdsX+W5RZQbY9JW5SePEvSMJ73K3seKnm1V1Fw7LUN50S8u0lNEZojIjPXr15cxRGNMYqqaL/Rq1aqxYcOGlE4gqsqGDRuoVq1apb1HWWdb/SIijVR1ndv1lOeW5wJN4+o1Ada65e33Kp/kljcpon6RVHUwMBhiA+ZljN0Yk4gKv8wr+Uu9SZMm5Obmkuo/QKtVq0aTJk32X7GMypo8xgPdgEHu7bi48ttEZASxwfHNboL5DHgsbpC8I9BLVfNFZKuInAZMA/4K/LuMMRljklrV/B7MyMigWbNmVfJeqWy/yUNEhhNrNRwsIrnEZk0NAkaJSA9gNXClW/1jYtN0lxKbqtsdwE0SA4Dpbr1HVLVwEP7/2D1V9xP3zxiTbna1OKxTIRnsN3mo6jXFPNShiLoK3FrM67wBvFFE+QzguP3FYYxJdZY0koktT2KMSSwpPJCdSix5GGMSgyWNpGLJwxiTIGzMI5lY8jDGJIYKankUXpO8oq9NbvZkycMYkyDKf57H3gnDEkjlseRhjEkM1luVVCx5GGMSRMnHPPbXorAWR+Wz5GGMSQwl6K4qa1LIys4p9rnx5ZZ0Ss6ShzEmsZRizKMsX/w2oF4xLHkYYxKDRktVvaSJw5JK5bDkYYxJDNGwu7Fny6O0X+YHsp2jZQ2n++bRyfcd5/u+42zfLI6VlZyQPXKP1yxJV5YpWllX1TXGmIrlhIt9qPhWhnJ2r9fp5p/LH30/0kpWcphv30utr3HqM0OPZkL0JCY6bSggWN7I05IlD2NMQti+YyfVhRKNeTSTdVzi/5ouvilk+X4BYLVTn7l6JCPCZ7NKD+FXPYjN1EBQMgnTUDZymPzCib6lnOH7gUv937BVD2B49BxOyd4I1NnjPbKyc1g56KJK2NPUYMnDGOO5rOwcFmdGALjm9WlMdbYB7PHl7SfK+b7p3Bj4lLa+xURV+MY5jtfDFzLZac1qbbjvNynMSVHw4XCqbwFX+7/kRv8ndPN/zsvRP/NypLO1RErIkocxxlOF3VBBiQ2Y+3H2eCyTEH/xT+TGwCc0kV9Z6TRkYPhaxkXbkbdXa6GkHHxMdY5lqnMsT8uV3BcYxV2BMXT2TeH28B3M16w9YrAWyO9Z8jDGJJTq7AQgQIQr/f/jjsAHNJJ8pjkteCR8Pf91TsKpwLk+a7Qhd4RvZ2S0PU9lvMqYYD8ejNzAqOjZFfYeqciShzEmAewe56hOAaf5fmRgYAhH+tbxvXMUd4f/zrdOq0qN4BvneC4qeIx/ZbzIkxmv0Yh80Asr9T2TmSUPY0yVKa4LqJn8vGv7X8GXAFjlNODm0D1McE4CpEriy6cWN4bv43F9nbsz3mfwgzuBa6vkvZONJQ9j0pEqSNV8IRcqbsXbA9hJ/8BbRFXwS6wF8kKkC/+OXOrJ4HWEAPdF/sZvVKNnIIdNWoOs7NixsrGP3ewkQWNMpSnqZLv4spayig+DfTnDN49ekZs4o+BfnFHwL56KXO3xrCehf+SvjImewf0Zo7jMNxmwkwfjWcvDGFPlsrI/4q/+z+kTeJdN1OC6cC+mOMd5HdYeFB/3h3tyCPk8ljGERaGmzNdmXoeVMCx5GJOOKrnbal+/0GuzlSczBtPRP5Mvoidyb/gW8qlVabGUR4QAt4Xv4MPMPrwafJYLCx73OqSEYd1WxphKV5hMTpUFfJzZi/a+2TwSvp4bw/clbOIolE8t/h66k0PIp1/G29Z15bLkYUxaqrrL9mVl5+Anyt2B0bwbfJSdGuSyUH/eiF5AVc2iKq85ehQvRi/hcv9XdPRN9zqchGDJwxhTqRqznhHBAdwZGMOY6JlcHHqMeXqE12GV2guRS5jnZDEw4w3YudnrcDxnycOYdBG/4GApLrhUEvFX6ovv1rnQ9y2fZPaihazhjtCt3Be5he1Uq9D3riphAmSHb6IeW2DSE16H4zlLHsakiwpOGPuSSYjHAq/xUvB5lmsjLgw9xninXZW9f2WZp0cwItqe8NRXIG+h1+F4ypKHMelCnfg7lfIWWdk5NJVfGBPsx7WBL3k58meuCPVjzf5WvE0iT0WuZjuZTPp3T69D8VS5koeI3C0i80VknogMF5FqItJMRKaJyBIRGSkiQbdupnt/qft4Vtzr9HLLF4nI+eXbJWNM0SouYcR3UcV3U3XwzSQn2IfG8ivdQ/fxROQaIil2RkA+tXgx0oX2/jlc0etpr8PxTJmTh4g0Bu4A2qrqcYAf6Ao8ATyrqs2BjUAP9yk9gI2qehTwrFsPEWnlPu9YoBPwkoj4yxqXMaYYlTDOUchPlPsCIxgSfJpV2oCLQwP50vlDhb5fInk72pH1ehD3ZrxHVvZHXofjifJ2WwWAA0QkAFQH1gHnAKPdx4cCl7jbXdz7uI93EBFxy0eoaoGqrgCWAqeUMy5jzN7iu60qMJHUYzNvZwzi1sB43o2czRWhh8nVBhX2+oloJ5m8ELmE03wLaOebB6Tf0iVlTh6q+hPwFLCaWNLYDMwENqlqxK2WCzR2txsDa9znRtz69eLLi3jOHkSkp4jMEJEZ69fv+zrFxpi9Vfw4RxtZzEeZfTjJt5h7w3+jd+TmtLkS3/DoOfysdfg///i0SxxQvm6rOsRaDc2AQ4EawAVFVC38jy3qbCDdR/nvC1UHq2pbVW1bv3790gdtTDqrgAHz+DGOa/0TGRkcQIFmcFmoP6Ojf6qAIJNHiAzejHTiDP98WslKr8OpcuXptjoXWKGq61U1DIwBTgdqu91YAE2Ate52LtAUwH38ICA/vryI5xhjKooT2X+dIuz9qzqDCI8GhvBYxhC+do6jc+hRfoy7bGs6GR49h21ajZsCH3sdSpUrT/JYDZwmItXdsYsOwI/Al8AVbp1uwDh3e7x7H/fxL1RV3fKu7mysZkBz4LtyxGWMKUo0LnmUcsyjsMVxMJsZFhzIdYGJvBTpTI/wfWyhRgUHmjy2UIOR0bP5s28qjdjgdThVqjxjHtOIDXx/D/zgvtZg4AHgHhFZSmxMY4j7lCFAPbf8HiDbfZ35wChiiedT4FZVjZY1LmNMMZxw3HbJWiHxrY5jZQXjMvtyvKzgjtBtPBnpWqHXEk9Wb0Y74ceha+ALr0OpUuWagK2q/YB+exUvp4jZUqq6E7iymNcZCAwsTyzGmP2IxiePcPH1itDZN4UnM15lA7W4ItTPrmsRJ1frM9lpzVX+/8Vad/7UOq+lOPazwZh0EZ8wonu2PIq74p8Ph+zAuzwffIE5eiSdCx61xFGE4dFzaCT5sHSC16FUmfRIkcYYcOJ6g4tpecQnkVps4/mMF2nvn8N/IufySOSvhO0ro0gTnT+Qp7WZ886TnDegqEmnqcf+E4xJF/HdVtF9d1sdLj/zRsY/aSp59Ar3YHi0QyUHl9wiBHgveha3+D+ELWuh1qFeh1TprNvKmHSxnwHzwlbHKbKAscGHqCNb+UuojyWOEhod/RN+UQY88ajXoVQJa3kYky5Cv+3eLma21eW+yTye8RprtAHdw/ezOoVWw61sK7QRc51mdPFP8TqUKmEtD2PSRcG23dtx3VZZ2TkIDvcGRvJ08BWmOy24NNTfEkcZjIueTmvfCs7u9brXoVQ6Sx7GpIvQ1t3bTnjXiX+ZhHgh43luC4zj3cjZdAs/wBYO9C7OJPZR9I84KnT2pX7rw5KHMeki/rrb7lTd+mxiZHAAF/im82j4L/SO3JRy19+oSr9Ql2+dlnT2T0n5pdoteRiTBrKyc2DzT7vu3/Dy5xwjq/kg8yGOllxuCd/F69GLKHqdUlMaHzl/5EjfOo6WXK9DqVSWPIxJF/nLcDSWHK71T2R0sD8BolwZeojPnZM9Di51/DfaBoDzfDNTeql2Sx7GpAE/UfLmfcknzsls0QPo6J/JCj2ELgUD7IzxCpZHHWY5R9HRP8PrUCqVJQ9jUljhL9+LfVNpIJsYH23H38L38Ej4eq4M9eMX6nocYWqaED2JE3zLOSSFV9q15GFMitrVZbJhGQ9nvM1s5wgmOCcx1TmWN6IXpM0V/7zwmdMWgHP933scSeWx5GFMCmvMeni7Cw7CXeFbbQn1KrJMD2W5cwgdfanbdWX/ScakiL0HZ5tLLiMzB7B50wauD/VipTbyKLJ0JEx02nCqbwGEtnsdTKWw5GFMisnKzuGa3k/yfvBhgkS4NtQ3bS8T66XJTmsyJcINDz/rdSiVws4GMibFXOL7miczXmWlHkL30P38RH2vQ0pL3zkt2KkZnOn7wetQKoW1PIxJAbEuK+WffW7mX8GXmOEcwxWhhy1xeKiAIN85LTjLNzclz/ew5GFMCggS5onAa9yXMYox0TPoFs5mCzW8DivtTXZa09z3E41ScMqudVsZk4SysnNYOegisrJzqMdm3gn+i1N8i3gucinPRq7AlhlJDJOd1sAwzvTP9TqUCmfJw5gklZWdQ0tZxWvBpzmYzdwWup2PnD96HZaJs1ib8LPW4awUHPew5GFMkjrfN51nMl5iCzW4ItSPeXqE1yGZ3xGmOMdypm8uqIKkTovQkocxSWL3oKtym38s92a8x2znSG4O3cN66ngamyned04LLvN/Db8ugfpHex1OhbHkYUwSqUYBT2YMprN/KmOiZ9ArfJMtM5LgpjktYxurvk6p5GGzrYxJYPFTPBuSz6jgI1zs+5ZB4a7cE/4/SxxJYIUeQp7WZuy497wOpUJZ8jAmAe19XsAlvZ7jw8y+HCHruDl8D69EO2MzqpKFMM1pwam+hbFxjxRhycOYBBKfNAq37+zdi5HBAezUDC4L9Weic5JX4Zkymua0pJHkw8YVXodSYcqVPESktoiMFpGFIrJARP4oInVFZIKILHFv67h1RUSeF5GlIjJXRNrEvU43t/4SEelW3p0yJpkVJg0/UXoHhvFc8CVm6VF0CQ1gsTb1ODpTFrvHPaZ4G0gFKm/L4zngU1VtAZwALACygYmq2hyY6N4HuABo7v71BF4GEJG6QD/gVOAUoF9hwjEmncS3OuqwhaEZg+gZyOGtSEeuD/ViI7U8jM6UxxJtzAatackDQERqAWcBQwBUNaSqm4AuwFC32lDgEne7C/C2xnwL1BaRRsD5wARVzVfVjcAEoFNZ4zIm2R0rK/kwsy8n+xZzb/hvPBy5gbBNjExywvdOc8id7nUgFaY8/5FHAOuBN0XkBGAmcCfQUFXXAajqOhFp4NZvDKyJe36uW1Zc+e+ISE9irRYOO+ywcoRuTOKIb3F09n3DExmvsZEDuTL0EHP1SA8jMxVplnMU5/06CnZshAOSv3OlPN1WAaAN8LKq/gH4jd1dVEUpamqI7qP894Wqg1W1raq2rV/fVgs1yS9+fKNv4D88H3yRuXoEfy4YaIkjxczS5rGNn2Z6G0gFKU/yyAVyVXWae380sWTyi9sdhXubF1c/frSvCbB2H+XGpKz41kZdtvCfjMe5KfAJb0bO5y+h3mzgIA+jM5VhrnMEjgrkpsalacucPFT1Z2CNiBzjFnUAfgTGA4UzproB49zt8cBf3VlXpwGb3e6tz4COIlLHHSjv6JYZk9KysnM4VlYwPrMvJ/mW8I/QLfSPdCNi4xsp6TcOYJE2SZlxj/L+l94ODBORILAc6E4sIY0SkR7AauBKt+7HwIXAUmC7WxdVzReRAUDhEX1EVfPLGZcxCe8S39cMyniNDdTiclvYMC3Mco6i0ZKp1HYc8CX3aXblSh6qOhtoW8RDHYqoq8CtxbzOG8Ab5YnFmERV2EW1ctBFsYJIiH6BoXQPfMa3TktuDd1h3VRpYpY251r5EvKXwcHNvQ6nXJI79RmTRLKyc2DzT3z/yB/pHviMIZELuC7UyxJHGpnlHBXbSIGuK0sexlSi+IHxdr4f2PDMqRwtufw9dAcDItfb+EaaWaaHslUPSIkZV5Y8jKlkgsOt/rH8J2MQG7QWnUOP8rFzmtdhGQ8oPn7Uw2Fd8l+W1n72GFOJarGNZzNepoN/FmOjp9M7fBPbqeZ1WMZD85xmnPrzJHCi4PN7HU6ZWfIwppJc3Ovf5ASfo6Hk0zfcnXei52LLqJt5ThZEdsSuLNighdfhlJl1WxlTAQrHNrKyc8jK/ghmvsX7wf74xOGqUD/eiZ6HJQ4DME+bxTbWzfE2kHKy5GFMBcnKzqEaBfwz8Cp8eCfTnBZcXDCQ2XqU16GZBLJcG7FDg5Y8jElX8a0NgCPlJ8YGH+Jy/1c8F7mMG8IP2DLq5nei+FmghyV98rAxD2PKYO/EcblvMgMy3mQ7mdwQvp/JzglehmcS3DynGUet/JpaSXymeXJGbUyCqM5Ons54maeDrzDHOZILCx63xGH2a55mUUt2JPVlaa3lYUwZtZDVvJDxPEfIOp4NX86/o5fi2O8xUwLznbhB83rJufS+/acbU0pZ2R9xjX8iY4MPUku285dwb56LXm6Jw5TYYm1CWP28OOIDr0MpM2t5GFMKx2W/x/MZQ+jsn8rk6PHcHf67rU1lSi1MgOXaiGNkzf4rJyj7qWTMfuxan2rtbD4M9uFC3zSeDF9Nt/ADljhMmS3SprTwJW/ysJaHMcUoTBqCw4A+f+f+wEgypRZdQ32Zocl7ZrBJDAudpnT2T4WdW6Ba8k3ptuRhTBEKE0d9NvJ0xiuc5f+Bz6MncX+4J5uo6XF0JhUsVvfq2+sXQtNTvA2mDKzbyhhX/PLpAB18M/k0M5uTfYvoFe5Bz/A9ljhMhVlYmDx+me9tIGVkLQ9j4hQuMTIgMIzrA/9lvnM4d4RvY5k29jo0k2J+0oPZptUYPfZjbmjb3etwSs2ShzHsbnW0lFU8n/ECzX0/MThyEU9FriJEhsfRmVSk+FisTThGcr0OpUwseZi0l5Wdg+Bwo/9T7g+MYBMHcl2oF187x3sdmklxi5ymnO+fDqogybXqso15mLSWlZ1DfTYyNOMJHsx4h8nOCXQqGGSJw1SJRdqUurINtv3idSilZsnDpJW9B8Uv9H3L55kPcLJvEX3CN3Jz+B5bCddUmUVJPGhuycOkvL0TBkDr7JGMffACXgo+zyptwMWhgQyzK/2ZKrbQcZNH3gJvAykDG/MwaaMwify192N8njmYemzh6fAVvBTtQpTkvZa0SV4bqcUGrUm9Xxd7HUqpWfIwaSErO4cD2EnvwLtcH/gvi53G3BT+B/P0CK9DM2luqTa25GFMomoji3km42UOkzwGRy7i6ciVFBD0OixjWOY04qhVM6jndSClZGMeJuXEj3EcnT2W+wMjeC/YHz8O14T68ljkL5Y4TMJYpodST7bCbxu8DqVUyt3yEBE/MAP4SVUvFpFmwAigLvA9cL2qhkQkE3gbOAnYAFytqivd1+gF9ACiwB2q+ll54zLpLSs7h2NlBeOCr9LSt5rhkbN5NHIdv3GA16EZs4dlemhsY8MSqJE87Y+KaHncCcRPFXgCeFZVmwMbiSUF3NuNqnoU8KxbDxFpBXQFjgU6AS+5CcmYEotvbQQJc29gJOOCD1JXttAj9A96RW62xGES0tLCpW+SbNyjXMlDRJoAFwGvu/cFOAcY7VYZClzibndx7+M+3sGt3wUYoaoFqroCWAok3xKTxnNZ2TmQO4OcYG9uC4zjg+gZnFfwJBOdk7wOzZhirdWD2akZsH6R16GUSnm7rf4F3A+7lhqtB2xS1Yh7PxcoXFGuMbAGQFUjIrLZrd8Y+DbuNeOfswcR6Qn0BDjssMPKGbpJJZmEuCfwHtHXPqaG1KFb6AH+55zgdVjG7JeDjxXaiJa/LvE6lFIpc/IQkYuBPFWdKSLtC4uLqKr7eWxfz9mzUHUwMBigbdu2RdYx6eeKXk/zSXAwR/h+ZlikA49HrmEb1b0Oy5gSW6qH0jLJuq3K0/JoB3QWkQuBakAtYi2R2iIScFsfTYC1bv1coCmQKyIB4CAgP668UPxzjClWy+z3uT8wklHBz/lJD+baUG+mOMd5HZYxpbZMD4VN30F4J2RU8zqcEinzmIeq9lLVJqqaRWzA+wtV/QvwJXCFW60bMM7dHu/ex338C1VVt7yriGS6M7WaA9+VNS6T2nYNjC+fxGfBB+ge+Iyh0Y6cH3rCEodJWsucQ0EdyF/mdSglVhknCT4AjBCRR4FZwBC3fAjwHxFZSqzF0RVAVeeLyCjgRyAC3Kqq0UqIy6SAOmyBD26BOcOJ0pArCx5iul1P3CS5PWZcNTzW22BKqEKSh6pOAia528spYraUqu4Erizm+QOBgRURi0lNWdkfsbLrViZm3kt49g5eiV7CC5FL7GQ/kxKWayMcFXzrk2fcw84wNwlrVxfVhmW8k/EYjL2FFdqIi0KP8XTkKkt4cTbCAAAQiElEQVQcJmUUEGQt9SB/udehlJitbWUSVgYRmPwUOycOorUvQJ/wjbwbPQe13zwmBa1yGtIkzcc8jCm3y3o9y4fBIfDFGiY6p9A/3I086ngdljGVZqUeQrv82V6HUWL2E854bo+LNe3YxH/6Xs7oYH9qynZ6hP7BreG7LHGYlLdSG8L2DbBjk9ehlIglD+OZPZKGKvf0foD1g47nWv9E3oqeT0dbWsSkkVXaMLaRJOMe1m1lPJWVncMxshrevJBnglOY5RzFDeEHmK/NvA7NmCq1QhvFNvKXQ+M23gZTApY8jCeysnOowQ7uCrxPd/+n5K+qzhORmxkV/ZMNiJu0tFobxDaSpOVhn1JTpbKyc0CVP/umMDHzXnr4P2FUtD3nFDzNyOjZljhM2iogyE+aPNN17ZNqqkTh+MaR8hPfPNSOfwdfIE9rc2moP70jN7Fp18LMxqSvVU5D2JAc03Wt28pUiers5KW+1/NJMIcdZNI33J13ox1w7PeLMbus1EM4PX+u12GUiH1yTYXbYxaV48CcEXyZeQ9/D4xnXLQd5xQ8zTvR8yxxGLOX2HTdX2HnZq9D2S9reZjKkzsDPnkAfprBOj2SW0J3M0ubex2VMQlrlR4S28hfDof+wdtg9sN++pkKU9jiaMBG3n/wYni9A7/kLuWe0C1cGupvicOY/VhRmDySYNzDWh6mwmQSgq+e5svMJwgQ5cVIZ16KdOE3DvA6NGOSwu7puiu8DaQELHmY8lPlb336MyH4Dkxcz1fOyQyMXMuawjNmjTElspNM1mpdDk2CBRIteZjyWTeXb166hVeD81noNLVLwRpTTqucQzg0Cc71sORhyuS07Le5N+M9LvN9RStfDfqGuzM8eg5R/F6HZkxSW6P1OWL1HBK93W7Jw5TOzi38+9HbmZSZg6AMjl7ES5EubKGG15EZkxJWawMayiYI74CMxB0vtORhSiYapu9D93FX4H1uD2xhbPR0nopcTa7W9zoyY1LKmsLP1KbVUP8Yb4PZB0seZt9UuanPAHoF3uXRjHVMc1rQPXQ/P+gRXkdmTEpaUzjjauMqSx4m+WRl59BaltEnYxivBxeyzGnETaF/8F+nDSBeh2dMytqdPFZ6Gsf+WPIwv/OnXkN4PuM9Ovun8qvWom+4OyOiZxOxfxdjKt16DmKHBjlg0yqvQ9kn+zYwu21Zx7Anb+W/wUmECfBCpAuvRP7MNqp7HZkxaURYo/U52loeJuFtz+eVx+6gm/8zrvJHeTd6Di9ELmU9tb2OzJi0tEYbcPRGa3mYBNUyO3YVv1sCH9HTv4OxTjuejVxuZ4Yb47HV2gA2TQVVkMQcY7TkkWaysnPIIMKSy/KYnDmA+rKZCdE2PBW5ikV6mNfhGWMgNgW+YAvs2AjV63odTpEseaQTJ0oX39f8I/AefLye5dqCW0J3MVMTdzqgMelodfyMqwRNHmVekl1EmorIlyKyQETmi8idbnldEZkgIkvc2zpuuYjI8yKyVETmikibuNfq5tZfIiLdyr9bJl6z7A9h3hgWP3w8zwVfYivVuSF0P1eHHrTEYUwCSobpuuVpeUSAf6jq9yJSE5gpIhOAG4CJqjpIRLKBbOAB4AKguft3KvAycKqI1AX6AW0BdV9nvKpuLEdshljSON83g0+C78PoNQiNuS10OznOqahdysWYhLX7LPPEHTQvc/JQ1XXAOnd7q4gsABoDXYD2brWhwCRiyaML8LaqKvCtiNQWkUZu3Qmqmg/gJqBOwPCyxpbusrI/oqNvBjnBMbTyrWKZ04g7QrfxkXOaXfrVmCTwGwewQWtSL0VbHruISBbwB2Aa0NBNLKjqOhFx2180BtbEPS3XLSuu3JSWKj36DODD4Psc71vJCqchd4X+znjndEsaxiSZNVqfegk8XbfcyUNEDgTeB+5S1S1S/LSyoh7QfZQX9V49gZ4Ahx1mM4N2UYUlE5jzzgMMCS5nldOAf4RuYazTzpZINyZJrdEGnJiqLQ8RySCWOIap6hi3+BcRaeS2OhoBeW55LtA07ulNgLVuefu9yicV9X6qOhgYDNC2bdsiE0xaUaV7n4HcEfiAP/iWUpf63BfuyQfRM2wpEWOS3BptAJtnghMFX+L9CCzPbCsBhgALVPWZuIfGA4UzproB4+LK/+rOujoN2Ox2b30GdBSROu7MrI5umSmO48CP45jX70TeDP6T+rKJ7PBNnBN6mvei7S1xGJMC1mh9cMKwZa3XoRSpPN8y7YDrgR9EZLZb1hsYBIwSkR7AauBK97GPgQuBpcB2oDuAquaLyABgulvvkcLBc7OXaIQ7H3yQWwPjONr3E9U5hHvDf2NstJ0lDGNSzB7netRuus+6XijPbKuvKX5t7g5F1Ffg1mJe6w3gjbLGkvIiIZgznFXjHuW5YB4LnabcHrqNHJs9ZUzK2nWux6ZVwJmexlIU+7mawI7J/oBFl+bBN8/Blp/YxBE8GrqO/zpt7DwNY1LcWq1HVAV//gqvQymSJY9EtHMzzHiTrzOfgU82M81pwQuRbL5yjscuxGRMeogQYK0eTNMEPVHQkkci2bIOpr3Mlq8HU0t2sMA5nlsjl/CdtvQ6MmOMB1ZrA5om6LkeljwSQd5CRj1/P1cFp4BG+Z9zCq9GLmaeXSfcmLS2WhvQbuN8r8MokiUPr6jC6m+Z8HpvzvN/z5/9Qd4Knc2Q6AV2PQ1jDOBO1/0tD0K/QbCG1+HswZJHFWuW/SErusH3wx+mjW8pJ/kO5Nnw5bwdPY+N1PI6PGNMAtm9uu4qaNjK22D2YsmjqoS2w9wR/Df4JIxcx8HU58HwDbwX/RM7yfQ6OmNMAtrjXA9LHunlj9lvM7XDMjZ9/Rq15Te2k8Vtodv5xDnF1p0yxuzT6j3O9Ugsljwqgyqs+Q6mvcxXmeNgCkxx2vJmpBPT9Rhsuq0xpiQ2UpOtegA1E3CBREseFSkSgh/HMnv045zoWw6ZBzEkegH/CXWMXZPYGGNKRcjV+rS05JGifvsVZrzJL1+8QEPZRE0a0TfcnTE7z2Q71byOzhiTxFZrA0seKeeX+Yz4d28u9X9DpoRZ6LTmgWhP/ue0tuVDjDEVYrU2gI3zY93hxV8vqcpZ8igtJwqLP+ObYY/Szj+fLv4g70XP4s1oJ5apXQDRGFOxVmsDiOyAbXlQM3HOAbPkUVI7t8DsYaz6+BkO9+XRzFeXx8PXMCJ6Nps50OvojDEpak38dF1LHklkwzL4bjBbvx1KTdlBHkfzRKgrnzkn21RbY0yl2+Ncj8NO9TSWeJY8irNuDkx+CufHD4ngY4LzR96MdOIHW2/KGFOFftKDYxsJNmhuyWNva2fBl4/Bks/ZogfwdrQzQyMdWU8dryMzxqShAoKs07o0SrATBS15FNr6C0zsD7OHka8HMiRyFf+JnscWEmsxMmNM+lmtDWhkLY8EtOBDNo64hers5M3on3kh0oVtVPc6KmOMAWKD5qduXOZ1GHtI7+ShCv97EiY9xhptxl3hW1muh3odlTHG7GG10wC2fAWRAggkxkKq6X0m2/+egEmPMTp6FpeH+lviMMYkpNiMK4VNa7wOZZf0TR6LPoFJjzM6ehb3hv9GOM0bYcaYxLWmcG28BBr3SM/kEfoNcu5lgXMYvcM9sFVujTGJbPe5Hiu8DSROeiaPWcNgSy4PhW8gRIbX0RhjzD6tpzbbNRPyl3sdyi7pmTxmDGG2cyTTtYXXkRhjTAkIS7Qx5C3wOpBd0i955K+A9QsZG23ndSTGGFNiS7SJJQ9PrZgMwGSntceBGGNMyS1ymsC2n2F7vtehAOmYPNYvZLtmskIP8ToSY4wpsSXaJLaxfqG3gbgSJnmISCcRWSQiS0Uku9LeaMMyVuohdrEmY0xSWeQ0jW38Mt/bQFwJ8Q0qIn7gReACoBVwjYi0qpQ327mJjWrX3zDGJJd11OVnrQOrvvE6FCBBkgdwCrBUVZeraggYAXSplHcq2MZvdl1xY0zSEb6KHg9LJybEuEeinFbdGIg/7z4X+N1VT0SkJ9DTvbtNRDYAv5btLb8s29MS08GU+TikHDsWu9mxiEmZ43BV4UbvemV9iYOBwysilkRJHkWd4q2/K1AdDAze9SSRGaratjIDSwZ2HHazY7GbHYsYOw67ucciqyJeK1G6rXKBpnH3mwBrPYrFGGPMfiRK8pgONBeRZiISBLoC4z2OyRhjTDESottKVSMichvwGeAH3lDVksxHG7z/KmnBjsNudix2s2MRY8dhtwo7FqL6u6EFY4wxZp8SpdvKGGNMErHkYYwxptSSMnlU2VImHhKRN0QkT0TmxZXVFZEJIrLEva3jlouIPO8ej7ki0ibuOd3c+ktEpJsX+1IeItJURL4UkQUiMl9E7nTL0/FYVBOR70Rkjnss+rvlzURkmrtfI91JJ4hIpnt/qft4Vtxr9XLLF4nI+d7sUfmIiF9EZonIR+79dD0OK0XkBxGZLSIz3LLK/3yoalL9ERtQXwYcAQSBOUArr+OqhP08C2gDzIsrexLIdrezgSfc7QuBT4idL3MaMM0trwssd2/ruNt1vN63Uh6HRkAbd7smsJjYEjbpeCwEONDdzgCmufs4Cujqlr8C/J+7/XfgFXe7KzDS3W7lfm4ygWbu58nv9f6V4XjcA7wLfOTeT9fjsBI4eK+ySv98JGPLo+qWMvGQqk4G9l6DoAsw1N0eClwSV/62xnwL1BaRRsD5wARVzVfVjcAEoFPlR19xVHWdqn7vbm8FFhBbkSAdj4Wq6jb3bob7p8A5wGi3fO9jUXiMRgMdRETc8hGqWqCqK4ClxD5XSUNEmgAXAa+794U0PA77UOmfj2RMHkUtZdLYo1iqWkNVXQexL1XAvbBxscckpY6V293wB2K/uNPyWLhdNbOBPGIf8GXAJlWNuFXi92vXPruPbwbqkRrH4l/A/YDj3q9Heh4HiP2A+FxEZkpsCSeogs9HQpznUUolWsokzRR3TFLmWInIgcD7wF2quiX2w7HoqkWUpcyxUNUocKKI1AY+AFoWVc29TcljISIXA3mqOlNE2hcWF1E1pY9DnHaqulZEGgATRGRfF/yosGORjC2PdF7K5Be3iYl7m+eWF3dMUuJYiUgGscQxTFXHuMVpeSwKqeomYBKxfuvaIlL4QzB+v3bts/v4QcS6QpP9WLQDOovISmLd1ucQa4mk23EAQFXXurd5xH5QnEIVfD6SMXmk81Im44HCWRDdgHFx5X91Z1KcBmx2m6qfAR1FpI4726KjW5Y03L7pIcACVX0m7qF0PBb13RYHInIAcC6xMaAvgSvcansfi8JjdAXwhcZGR8cDXd1ZSM2A5sB3VbMX5aeqvVS1icYW+OtKbL/+QpodBwARqSEiNQu3if1fz6MqPh9ezxQo4+yCC4nNulkG9PE6nkrax+HAOiBM7FdBD2L9tBOBJe5tXbeuELuY1jLgB6Bt3OvcSGwgcCnQ3ev9KsNxOINY83kuMNv9uzBNj0VrYJZ7LOYBD7nlRxD70lsKvAdkuuXV3PtL3cePiHutPu4xWgRc4PW+leOYtGf3bKu0Ow7uPs9x/+YXfh9WxefDlicxxhhTasnYbWWMMcZjljyMMcaUmiUPY4wxpWbJwxhjTKlZ8jDGGFNqljyMMcaUmiUPY4wxpfb/PF3sT8+FSGsAAAAASUVORK5CYII=\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[:nevents], bins = bins, range = (x_min,x_max), label = 'data')\n",
"plt.plot(test_q, calcs_test*nevents*4.5 , label = 'pdf')\n",
"\n",
"# plt.plot(sam, calcs, '.')\n",
"# plt.plot(test_q, calcs_test)\n",
"# plt.yscale('log')\n",
"plt.ylim(0, 12000)\n",
"# plt.xlim(3080, 3110)\n",
"\n",
"plt.legend()\n",
"\n",
"plt.savefig('test2.png')"
]
},
{
"cell_type": "code",
"execution_count": 25,
"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": 26,
"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 = 18537900.293873098</td>\n",
" <td title=\"Total number of call to FCN so far\">TOTAL NCALL = 194</td>\n",
" <td title=\"Number of call in last migrad\">NCALLS = 194</td>\n",
" </tr>\n",
" <tr>\n",
" <td title=\"Estimated distance to minimum\">EDM = 0.03059806673262826</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:#FF7878\">False</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:#FF7878\">True</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=\"$('#kJHicCUbKd').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>jpsi_s</td>\n",
" <td>464.524</td>\n",
" <td>0.229484</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>76.5024</td>\n",
" <td>0.0505517</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>psi2s_p</td>\n",
" <td>0.30538</td>\n",
" <td>0.0245293</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_p</td>\n",
" <td>-9.42323</td>\n",
" <td>0.0114089</td>\n",
" <td></td>\n",
" <td></td>\n",
" <td></td>\n",
" <td></td>\n",
" <td>No</td>\n",
" </tr>\n",
"</table>\n",
"<pre id=\"kJHicCUbKd\" 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 & $jpsi_{s}$ & 464.524 & 0.229484 & & & & & No\\\\\n",
"\\hline\n",
"1 & $psi2s_{s}$ & 76.5024 & 0.0505517 & & & & & No\\\\\n",
"\\hline\n",
"2 & $psi2s_{p}$ & 0.30538 & 0.0245293 & & & & & No\\\\\n",
"\\hline\n",
"3 & $jpsi_{p}$ & -9.42323 & 0.0114089 & & & & & No\\\\\n",
"\\hline\n",
"\\end{tabular}\n",
"</textarea>\n",
"</pre>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
"<hr>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"ename": "RuntimeError",
"evalue": "Function mimimum is not valid. Make sure migrad converge first",
"output_type": "error",
"traceback": [
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[1;31mRuntimeError\u001b[0m Traceback (most recent call last)",
"\u001b[1;32m<ipython-input-26-17821450db47>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m 5\u001b[0m \u001b[0mresult\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mminimizer\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mminimize\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnll\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 6\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 7\u001b[1;33m \u001b[0mparam_errors\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mresult\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0merror\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 8\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 9\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mvar\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0merrors\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mparam_errors\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mitems\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32mc:\\users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\minimizers\\fitresult.py\u001b[0m in \u001b[0;36merror\u001b[1;34m(self, params, method, error_name, sigma)\u001b[0m\n\u001b[0;32m 227\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 228\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0muncached_params\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 229\u001b[1;33m \u001b[0merror_dict\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_error\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mparams\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0muncached_params\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mmethod\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mmethod\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0msigma\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0msigma\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 230\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_cache_errors\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0merror_name\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0merror_name\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0merrors\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0merror_dict\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 231\u001b[0m \u001b[0mall_errors\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mOrderedDict\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mp\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mparams\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mp\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0merror_name\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mp\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mparams\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32mc:\\users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\minimizers\\fitresult.py\u001b[0m in \u001b[0;36m_error\u001b[1;34m(self, params, method, sigma)\u001b[0m\n\u001b[0;32m 238\u001b[0m \u001b[1;32mexcept\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 239\u001b[0m \u001b[1;32mraise\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"The following method is not a valid, implemented method: {}\"\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmethod\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 240\u001b[1;33m \u001b[1;32mreturn\u001b[0m \u001b[0mmethod\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mresult\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mparams\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mparams\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0msigma\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0msigma\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 241\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 242\u001b[0m \u001b[1;31m# def set_error_method(self, method):\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32mc:\\users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\minimizers\\fitresult.py\u001b[0m in \u001b[0;36m_minos_minuit\u001b[1;34m(result, params, sigma)\u001b[0m\n\u001b[0;32m 46\u001b[0m \"`MinuitMinimizer`.\")\n\u001b[0;32m 47\u001b[0m result = [minimizer._minuit_minimizer.minos(var=p.name, sigma=sigma)\n\u001b[1;32m---> 48\u001b[1;33m for p in params][-1] # returns every var\n\u001b[0m\u001b[0;32m 49\u001b[0m \u001b[0mresult\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mOrderedDict\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mp\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mresult\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mname\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mp\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mparams\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 50\u001b[0m \u001b[1;32mreturn\u001b[0m \u001b[0mresult\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32mc:\\users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\minimizers\\fitresult.py\u001b[0m in \u001b[0;36m<listcomp>\u001b[1;34m(.0)\u001b[0m\n\u001b[0;32m 46\u001b[0m \"`MinuitMinimizer`.\")\n\u001b[0;32m 47\u001b[0m result = [minimizer._minuit_minimizer.minos(var=p.name, sigma=sigma)\n\u001b[1;32m---> 48\u001b[1;33m for p in params][-1] # returns every var\n\u001b[0m\u001b[0;32m 49\u001b[0m \u001b[0mresult\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mOrderedDict\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mp\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mresult\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mname\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mp\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mparams\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 50\u001b[0m \u001b[1;32mreturn\u001b[0m \u001b[0mresult\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32miminuit\\_libiminuit.pyx\u001b[0m in \u001b[0;36miminuit._libiminuit.Minuit.minos\u001b[1;34m()\u001b[0m\n",
"\u001b[1;31mRuntimeError\u001b[0m: Function mimimum is not valid. Make sure migrad converge first"
]
}
],
"source": [
"nll = zfit.loss.UnbinnedNLL(model=total_f, data=data3, 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": null,
"metadata": {},
"outputs": [],
"source": [
"(-7.95933+2*np.pi)/np.pi+np.pi"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"display_time(int(395*pdg[\"number_of_decays\"]/100000))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(display_time(22376))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# 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": null,
"metadata": {},
"outputs": [],
"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, 5e-4)\n",
"# plt.yscale('log')\n",
"# plt.xlim(3080, 3110)\n",
"plt.savefig('test3.png')\n",
"print(jpsi_width)"
]
},
{
"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.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
