{
"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+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8VfWd//HXJzcbS9gDCGGJEoWguEXq1pZWW9GqdGEqdmyx1Trj6Pj7tdPp1Fm6ONNO/U0fdaatSxmxP/TXCtZx2rSltbVqx7YoRkEqCBpki6AEEsISstzk8/vjfsFrvDf3BnNyk/B+Ph73cc/5nu/5nu+Xc8k7Z7kn5u6IiIhEIS/XHRARkcFLISMiIpFRyIiISGQUMiIiEhmFjIiIREYhIyIikckqZMxsvpltMrNaM/tSiuVFZrYiLH/GzKYnLbs1lG8ys0t60OZ3zexgNtsQEZH+KWPImFkMuBO4FKgErjazyi7VrgMa3X0GcAdwe1i3ElgEzAbmA3eZWSxTm2ZWBYzKZhsiItJ/ZXMkMxeodfdX3b0NWA4s6FJnAbAsTD8MXGRmFsqXu3uru28BakN7adsMAfRvwBez3IaIiPRT+VnUmQzsSJqvA96Vro67x82sCRgbyp/usu7kMJ2uzZuBanff1SVD0m1jT3IlM7sBuAFg2LBhZ8+cOTOLIYpkZ3N94gzuSaXDc9yT3tPY3EZd42FOnlBCUb4u0wo899xze9y9tDfayiZkUh0tdH0WTbo66cpTfZLdzCYBfwbMO8Z+4O5LgCUAVVVVXlNTk2I1kZ5rbosz56u/5rPvOZG/mz94fnl58bUmLv/u7/nG1WdyxemTct0d6QfMbFtvtZXNry11wJSk+TJgZ7o6ZpYPjAQaulk3XfmZwAyg1sy2AkPNrDbDNkT6xDNbGoh3OuedODbXXelVJ08ooSBmrN+5P9ddkUEom5B5Fqgws3IzKyRxIb+6S51qYHGYXgg87oknb1YDi8KdYeVABbA6XZvu/gt3n+ju0919OtAcLvR3tw2RPvHkxt0UF+Qxt3xMrrvSqwrz86gYX8L6nU257ooMQhlPl4XrHzcDjwIx4D53X29mtwE17l4NLAUeCEcdDSRCg1DvIWADEAducvcOgFRtZuhKym2I9AV354lN9Vxw0jiKC2K57k6vO3XyCH6z4Q3cHd1PI70pm2syuPtKYGWXsi8nTbeQuJaSat2vA1/Pps0UdYYnTafdhkjUtuw5xPaGZj777vJcdyUSZ04dzUM1dWzd20z5uGG57k6/1t7eTl1dHS0tLbnuyjtWXFxMWVkZBQUFkW0jq5AROd49sakegHmnjM9xT6Jx1tTRADy/rVEhk0FdXR0lJSVMnz59QB/1uTt79+6lrq6O8vLofnnS/YoiWXh84xucVDqMKWOG5rorkagYP5ySonye296Y6670ey0tLYwdO3ZABwyAmTF27NjIj8gUMiIZ7D3YyqrNe7n01BNy3ZXI5OUZZ04bzfPbFDLZGOgBc0RfjEMhI5LBr9a/TqfDZacN3pABOGvqKDa9cYADLe257ooMIgoZkQxW/mkX5eOGMeuEklx3JVJnTxuNO6zZvi/XXZFe8OSTT3L55ZcD0NraysUXX8wZZ5zBihUr+rQfuvAv0o0jp8punHfSoDlFks5ZU0dTEDP+uHkv7zm5V54oIv3EmjVraG9vZ+3atX2+bR3JiHTjZy/spNPh8jmD/3Erw4ryOXPqaP5QuydzZcmprVu3MnPmTBYvXsycOXNYuHAhzc3N/OpXv2LmzJlceOGFPPLIIwDs3r2ba665hrVr13LGGWewefPmPu2rjmRE0nB3VtTUcdrkkcw6YUSuu9MnLpwxjjsee5nGQ22MHlaY6+70e1/72Xo29PLjeConjeArV8zOWG/Tpk0sXbqUCy64gM985jN8+9vf5vvf/z6PP/44M2bM4KqrrgJg/Pjx3HvvvXzrW9/i5z//ea/2NRs6khFJY/3O/by0az8fryrLdVf6zAUzxuEOq17dm+uuSAZTpkzhggsuAOCaa66hpqaG8vJyKioqMDOuueaaHPcwQUcyImk8VLODovw8rjxjcubKg8TpZSMZXpTP72v3DPq76XpDNkccUel6jbCpqalfXjfUkYxICoda4/z3mteYf+pERg6J7pEb/U1+LI/zTxrLkxt3o+fP9m/bt29n1apVADz44INcfPHFbNmy5eg1lwcffDCX3TtKISOSwo9rdnCgJc6150/PdVf63AcqJ7CzqUWP/u/nZs2axbJly5gzZw4NDQ187nOfY8mSJXzoQx/iwgsvZNq0abnuIqDTZSJv09Hp3PeHrZw1dRRnhmd6HU8umjWBPINfr3+dUyePzHV3JI28vDzuueeet5TNnz+fjRs3vq3uvHnzmDdvXh/17K10JCPSxWMvvcH2hmauu/DEXHclJ8YMK+Sc6WP49YY3ct0VGQQUMiJJ3J07n6ilbPQQLpk9IdfdyZkPzp7IxtcPsG3voVx3RVKYPn06L774Yq67kRWFjEiS3760m3V1Tdzy/gryY8fvf4/5p04EEl9GlbcbLDdF9MU4jt//RSJddHY63/7Ny0wbO5SPnHX83LacyuRRQzj3xDE88vxrg+YHam8pLi5m7969A/7f5cjfkykuLo50O7rwLxKsfHEXG3bt59sfP52C4/go5oiPnlnGF/9rHWt37Dsub4BIp6ysjLq6Ourr63PdlXfsyF/GjFJWIWNm84H/AGLAve7+zS7Li4D7gbOBvcBV7r41LLsVuA7oAG5x90e7a9PMlgJVgAEvA9e6+0Ezuxb4N+C1sNnvufu9xzZskbc63NbBN37xErNOGMGC4+jLl9259LSJ/NNPX+SR519TyCQpKCiI9C9JDjYZf10zsxhwJ3ApUAlcbWaVXapdBzS6+wzgDuD2sG4lsAiYDcwH7jKzWIY2P+fup7v7HGA7cHPSdla4+xnhpYCRXnP37zazs6mFr105m1he//vWdC6UFBcw/9SJ/GTtaxxqjee6OzJAZXNOYC5Q6+6vunsbsBxY0KXOAmBZmH4YuMgSzzdYACx391Z33wLUhvbStunu+wHC+kOAgX3iU/q9rXsOcc/vNnPl6ZOYWz4m193pVz513nQOtMR5ZM1rmSuLpJBNyEwGdiTN14WylHXcPQ40AWO7WbfbNs3sB8DrwEzgu0n1PmZm68zsYTObkkXfRbrV2en87cMvUJSfx99fNivX3el3zpo6ijllI1n2x60D/kK35EY2IZPq3EHXT1u6Oj0tT0y4fxqYBLwEXBWKfwZMD6fRHuPNI6e3dsTsBjOrMbOawXBhTqL1gz9u5dmtjXz1itlMHBntXTYDkZlx7fnTqd19kKde0d+ZkZ7LJmTqgOSjhjKg683zR+uYWT4wEmjoZt2Mbbp7B7AC+FiY3+vurWHxf5K4yeBt3H2Ju1e5e1Vpqf66n6S38fX9/J9fbeTiWeP56HF+y3J3PjTnBCaMKOJ7T9TqaEZ6LJuQeRaoMLNyMyskcSG/ukudamBxmF4IPO6JT2M1sMjMisysHKgAVqdr0xJmwNFrMlcAG8N88nPHryRxlCNyTPa3tHPj/3uekUMK+NePzumXj0jvL4ryY/zVvBms3tKgvzMjPZbxFmZ3j5vZzcCjJG43vs/d15vZbUCNu1cDS4EHzKyWxBHMorDuejN7CNgAxIGbwhEKadrMA5aZ2QgSp9ReAG4MXbnFzK4M7TQA1/bKv4Acdzo7nS/+eB3bG5pZfsO5lJYU5bpL/d5V50zhridr+fffvMJ5J45VKEvWbDAf/lZVVXlNTU2uuyH9zDd/uZF7freZf/zQLK5/9/H5EMxjcf+qrXz5p+u591NVXFx5/D7X7XhgZs+5e1VvtKWvNctx5f5VW7nnd5v583dN5boL9YW6nrh67lRmjB/OP/9iA63xjlx3RwYIhYwcN/57TR1fqV7PxbMmcNuCU3XKp4cKYnl85YpKtu1t5t6ntuS6OzJAKGTkuPDwc3V8/qEXOO/EsXz36jP1rf5j9O6KUj5YOYHv/PYVNtcfzHV3ZABQyMigd/+qrfztwy9w4YxxLF18DkMKY7nu0oD2Lx8+leKCGF/48QvEOzpz3R3p5xQyMmh1dDq3/WwDX/7pei6aOZ7//FSVAqYXjB9RzD9/+FTWbN/H3U9uznV3pJ/To/5lUNrX3MbnH3qBxzfu5tMXTOcfP1SpU2S96Io5J/Dbl97gjsde5vQpo3jPyfris6SmIxkZdJ7b1sBl//EUT71Szz8vmM1XrtCTlXubmfGvHz2NkyeU8NcPrmH73uZcd0n6KYWMDBqt8Q6+/ZuX+fj3nyY/lsd/3Xg+nzxveq67NWgNLczn+588G3fn2h+sZs/B1swryXFHISODwrNbE0cv3/ntK1wx5wR+fsuFzCkbletuDXrTxg5j6bXnsLPpMNf+YDUHWtpz3SXpZxQyMqC9tu8wn1uxlj+7ZxUt7Z384NPn8O+LzmREcUGuu3bcOGf6GO7+87PZuOsA1yxdTeOhtlx3SfoRhYwMSPua2/jXX77E+771JL/40y7+8r0n8evPvYf3nTI+1107Lr1v5njuvuZsXtq1n49/fxWvN7XkukvST+jZZTKg7Go6zNKntvCj1ds53N7BR86czN988BQmjxqS664JsGrzXj57fw1DC2Pcfc1ZnD1Nf2l0IOrNZ5cpZKTfc3ee376PHz2zneoXXqPT4crTJ/GX7z2JUyaW5Lp70sXG1/fzFw88x859h/mnyyv55LnT9AifAUYhkyWFzMC2r7mNn67dyYOrt7Px9QMMK4zxsbPL+Oy7T2TKmKG57p50o6m5nf+9Yg1PbKrnPSeXcvvHTuOEkTraHCgUMllSyAw8TYfb+c2GN/j5up38/pU9xDud0yaP5BPvmsoVp09ieJG+PzxQdHY6P3xmG99YuZH8mPE3HziZPz93GgUxXQru7xQyWVLIDAxb9hziyU27eXJTPas276Wto5PJo4Zw+ZwTuOL0SZw6eWSuuyjvwNY9h/jHn7zI72v3MGP8cP7+spm875TxOoXWjylksqSQ6Z/qD7RSs7WBp1/dy5Mv17MtfFu8fNww3j9zPJfPOYEzpozSD6FBxN157KXdfP0XG9i6t5nTJo/k5vfP4AOzJpCnpzH0OwqZLClkcq+j09my5xDPb2+kZmsDz25tZMueQwAUF+Rx/knjmHdKKe89uZRpY4fluLcStbZ4J/+9po67ntzMtr3NnFg6jE/MncrCs8sYNbQw192TQCGTJYVM32qLd1K7+yAv7mxi/WtNvLhzPxt27udwe+KvKI4aWkDVtDGcM30055SP4dRJIynM1/n541G8o5Ofr9vFslVbWbN9H4X5eVx66kQunzOJd1eMo7hAT8vOpT4PGTObD/wHEAPudfdvdlleBNwPnA3sBa5y961h2a3AdUAHcIu7P9pdm2a2FKgCDHgZuNbdD3a3jXQUMtHY39LO5t0H2Vx/iNrdB9lcf5DNuw+yraGZjs7E52lYYYzZk0Yye/IIZk8ayellIzmpdLhOjcjbvLRrf7g9fSdNh9spKcrn4soJXDRrPBecNI7Rw3SE09f6NGTMLEbih/0HgDrgWeBqd9+QVOevgDnu/pdmtgj4iLtfZWaVwIPAXGAS8BhwclgtZZtmNsLd94d2vw3sdvdvpttGd31XyBybtngnO/cdZkdjMzsaDlPX2MyOxsPsaGimrvHwWx6EWBAzpo8dxkmlw5kxfjgnTyzh1EkjmD52mAJFeqS9o5M/1O5h5Z928ej6N2g63I4ZzJk8kndXlHJO+RjOmDKKkUP0yKCo9WbIZHM/6Fyg1t1fDRtfDiwANiTVWQB8NUw/DHzPEldtFwDL3b0V2GJmtaE90rWZFDAGDAG8u234YD7f18s6O53G5jbe2N/KGwda2L2/hd1h+o39reze33J0WfK/an6eMWnUEMpGD+GimeOZNm4oM0qHc9L44UwdM1S3pEqvKIjlMe+U8cw7ZTzf+Egn615r4qmX9/DUK/Xc/bvNfO+JWgBmjB/OmVNGMWfKKGZOLOGUiSV6Vl0/lk3ITAZ2JM3XAe9KV8fd42bWBIwN5U93WXdymE7bppn9ALiMRJD9TYZt7EnuiJndANwAMHXq1CyGN3C1xTvZ19xGQ3MbDQcT742H2mg41E5jcxt7Dx2Zb6OxuY09B1tp73h7Jo8eWsCEEcWUlhRRMaGEySFQpowZypQxQ5lQUkS+gkT6UH4sj7OmjuasqaP5XxdXcLA1zrod+3h+eyPPb9/HYy+9wY+fqztaf9LIYk6ZWMLJE0qYOnYo08YMY+qYoUwaVazPbo5lEzKpznl0/UmVrk668lR7/Wib7v7pcJruu8BVwA+y7AfuvgRYAonTZSnW6RfiHZ00t3fQ3NrBgZZ29re0s78lzv7D7RxoibO/Jbx3M9/c1pG2/RHF+YwZVsjoYYWcMLKYykkjKC0pYnxJERNGFDNhRBHjSxLBoous0t8NL8rn/BnjOH/GOCBxS/Rr+w7z8hsH2Pj6ATaF1x9qE9+zOiI/z5g8eghTRg89+rlPvL85XVpSpKPxCGUTMnXAlKT5MmBnmjp1ZpYPjAQaMqzbbZvu3mFmK4C/JREy6bYRiXhHJy3xTlraO2gN78nTre2dNLd1cKgtzuHw3tzaQXNbB81tcQ61dXC4Lc6h1o4QJvG3LGuLd2bsQ0HMGFFcQElxPiOGJN7HlwynpDifkuICRg4pYMywwkSYDA3vwwoYPbRQ/2lkUDMzykYPpWz0UN4/c8LR8o5O5439LWzb28z2hkNs29vMtoZmXms8zKub97D7QCvxzrf/7llSlM+o8H8n8SpgVPg/NWpoAcOL8t98FeczrCifkqLE+9DCmL7T1Y1sQuZZoMLMyoHXgEXAJ7rUqQYWA6uAhcDj7u5mVg38KFzAnwRUAKtJHJW8rc1wHeYkd68N01cAG7vbRncdf2N/C99Y+dLRUGiJd3QJjU5a4520dgmQlnjn0bukeiLPEn8tcGhh7OiHb2hhjFFDCpg8qpghBfkMK4odrZN45b8lREYUFzAizBfl5+nDK9IDsXD9cNKoIZx30ti3Le/sdBqa23i9qYXd4Vpk/YFWGpvb2NfcTsOhNvY1t7FlzyEaD7VxoDWecZt5BsMKE4EzpDBGUX4eRQUxivPzKC6IUVyQeC86Op9YVhTKCmJ55MeMglgeBTEjPy/xnijPoyDPEu9HyxJ1CsN0LM/IMyPPEuM3O1JGKH9zPhc/TzKGTLj+cTPwKInbje9z9/VmdhtQ4+7VwFLggXBhv4FEaBDqPUTi2kocuMndOwDStJkHLDOzESSC6AXgxtCVlNvozu4DrSz749ajO7oo/607vKQ4n3H5MYoK8ijOT/VhePs6iQ9PYp1hSWExrChfoSDSz+XlGeOGFzFueBGJkyHda+/oZF9zO4da4xw88mqJc6gtzoGW+NvKW5LPerQnrpm2JP2Cm/jFNvGeCxaCJ2aGhVA6ElB5eUfKe/dnmL6MKSLSx9w9cRYl3km8o5N4p9MWT7zHOzpp73DaOzqJdyam42G+PdRt7+g8WtbhTqcnjtI63enodNwJ5R7KCeWeoj6Jev7m/Dc/NqdPb2EWEZFeZGZHT531R9/MXCVrujosIiKRUciIiEhkFDIiIhIZhYyIiERGISMiIpFRyIiISGQUMiIiEhmFjIiIREYhIyIikVHIiIhIZBQyIiISGYWMiIhERiEjIiKRUciIiEhkFDIiIhIZhYyIiERGISMiIpHJKmTMbL6ZbTKzWjP7UorlRWa2Iix/xsymJy27NZRvMrNLMrVpZj8M5S+a2X1mVhDK55lZk5mtDa8vv5OBi4hI9DKGjJnFgDuBS4FK4Gozq+xS7Tqg0d1nAHcAt4d1K4FFwGxgPnCXmcUytPlDYCZwGjAEuD5pO0+5+xnhdduxDFhERPpONkcyc4Fad3/V3duA5cCCLnUWAMvC9MPARWZmoXy5u7e6+xagNrSXtk13X+kBsBooe2dDFBGRXMkmZCYDO5Lm60JZyjruHgeagLHdrJuxzXCa7JPAr5KKzzOzF8zsl2Y2O1VnzewGM6sxs5r6+voshiciIlHJJmQsRZlnWaen5cnuAv7H3Z8K888D09z9dOC7wE9Sddbdl7h7lbtXlZaWpqoiIiJ9JJuQqQOmJM2XATvT1TGzfGAk0NDNut22aWZfAUqBzx8pc/f97n4wTK8ECsxsXBb9FxGRHMkmZJ4FKsys3MwKSVzIr+5SpxpYHKYXAo+HayrVwKJw91k5UEHiOkvaNs3seuAS4Gp37zyyATObGK7zYGZzQ9/3HsugRUSkb+RnquDucTO7GXgUiAH3uft6M7sNqHH3amAp8ICZ1ZI4glkU1l1vZg8BG4A4cJO7dwCkajNs8h5gG7AqZMoj4U6yhcCNZhYHDgOLQpCJiEg/ZYP553RVVZXX1NTkuhsiIgOKmT3n7lW90Za+8S8iIpFRyIiISGQUMiIiEhmFjIiIREYhIyIikVHIiIhIZBQyIiISGYWMiIhERiEjIiKRUciIiEhkFDIiIhIZhYyIiERGISMiIpFRyIiISGQUMiIiEhmFjIiIREYhIyIikVHIiIhIZLIKGTObb2abzKzWzL6UYnmRma0Iy58xs+lJy24N5ZvM7JJMbZrZD0P5i2Z2n5kVhHIzs++E+uvM7Kx3MnAREYlexpAxsxhwJ3ApUAlcbWaVXapdBzS6+wzgDuD2sG4lsAiYDcwH7jKzWIY2fwjMBE4DhgDXh/JLgYrwugG4+1gGLCIifSebI5m5QK27v+rubcByYEGXOguAZWH6YeAiM7NQvtzdW919C1Ab2kvbpruv9ABYDZQlbeP+sOhpYJSZnXCM4xYRkT6QTchMBnYkzdeFspR13D0ONAFju1k3Y5vhNNkngV/1oB+Y2Q1mVmNmNfX19VkMT0REopJNyFiKMs+yTk/Lk90F/I+7P9WDfuDuS9y9yt2rSktLU6wiIiJ9JT+LOnXAlKT5MmBnmjp1ZpYPjAQaMqybtk0z+wpQCvxFD/shIiL9SDZHMs8CFWZWbmaFJC7kV3epUw0sDtMLgcfDNZVqYFG4+6ycxEX71d21aWbXA5cAV7t7Z5dtfCrcZXYu0OTuu45hzCIi0kcyHsm4e9zMbgYeBWLAfe6+3sxuA2rcvRpYCjxgZrUkjmAWhXXXm9lDwAYgDtzk7h0AqdoMm7wH2AasStw7wCPufhuwEriMxM0DzcCne+MfQEREomOJA47BqaqqymtqanLdDRGRAcXMnnP3qt5oS9/4FxGRyChkREQkMgoZERGJjEJGREQio5AREZHIKGRERCQyChkREYmMQkZERCKjkBERkcgoZEREJDIKGRERiYxCRkREIqOQERGRyChkREQkMgoZERGJjEJGREQio5AREZHIKGRERCQyWYWMmc03s01mVmtmX0qxvMjMVoTlz5jZ9KRlt4byTWZ2SaY2zezmUOZmNi6pfJ6ZNZnZ2vD68rEOWkRE+kZ+pgpmFgPuBD4A1AHPmlm1u29IqnYd0OjuM8xsEXA7cJWZVQKLgNnAJOAxMzs5rJOuzT8APweeTNGdp9z98mMYp4iI5EA2RzJzgVp3f9Xd24DlwIIudRYAy8L0w8BFZmahfLm7t7r7FqA2tJe2TXdf4+5b3+G4RESkH8gmZCYDO5Lm60JZyjruHgeagLHdrJtNm6mcZ2YvmNkvzWx2qgpmdoOZ1ZhZTX19fRZNiohIVLIJGUtR5lnW6Wl5d54Hprn76cB3gZ+kquTuS9y9yt2rSktLMzQpIiJRyiZk6oApSfNlwM50dcwsHxgJNHSzbjZtvoW773f3g2F6JVCQfGOAiIj0P9mEzLNAhZmVm1khiQv51V3qVAOLw/RC4HF391C+KNx9Vg5UAKuzbPMtzGxiuM6Dmc0Nfd+bzSBFRCQ3Mt5d5u5xM7sZeBSIAfe5+3ozuw2ocfdqYCnwgJnVkjiCWRTWXW9mDwEbgDhwk7t3QOJW5a5thvJbgC8CE4F1ZrbS3a8nEV43mlkcOAwsCkEmIiL9lA3mn9NVVVVeU1OT626IiAwoZvacu1f1Rlv6xr+IiERGISMiIpFRyIiISGQUMiIiEhmFjIiIREYhIyIikVHIiIhIZBQyIiISGYWMiIhERiEjIiKRUciIiEhkFDIiIhIZhYyIiERGISMiIpFRyIiISGQUMiIiEhmFjIiIREYhIyIikckqZMxsvpltMrNaM/tSiuVFZrYiLH/GzKYnLbs1lG8ys0sytWlmN4cyN7NxSeVmZt8Jy9aZ2VnHOmgREekbGUPGzGLAncClQCVwtZlVdql2HdDo7jOAO4Dbw7qVwCJgNjAfuMvMYhna/ANwMbCtyzYuBSrC6wbg7p4NVURE+lo2RzJzgVp3f9Xd24DlwIIudRYAy8L0w8BFZmahfLm7t7r7FqA2tJe2TXdf4+5bU/RjAXC/JzwNjDKzE3oyWBER6VvZhMxkYEfSfF0oS1nH3eNAEzC2m3WzafNY+oGZ3WBmNWZWU19fn6FJERGJUjYhYynKPMs6PS1/p/3A3Ze4e5W7V5WWlmZoUkREopRNyNQBU5Lmy4Cd6eqYWT4wEmjoZt1s2jyWfoiISD+STcg8C1SYWbmZFZK4kF/dpU41sDhMLwQed3cP5YvC3WflJC7ar86yza6qgU+Fu8zOBZrcfVcW/RcRkRzJz1TB3eNmdjPwKBAD7nP39WZ2G1Dj7tXAUuABM6slcQSzKKy73sweAjYAceAmd++AxK3KXdsM5bcAXwQmAuvMbKW7Xw+sBC4jcfNAM/Dp3vpHEBGRaFjigGNwqqqq8pqamlx3Q0RkQDGz59y9qjfa0jf+RUQkMgoZERGJjEJGREQio5AREZHIKGRERCQyChkREYmMQkZERCKjkBERkcgoZEREJDIKGRERiYxCRkREIqOQERGRyChkREQkMgoZERGJjEJGREQio5AREZHIKGRERCQyChkREYlMViFjZvPNbJOZ1ZrZl1IsLzKzFWH5M2Y2PWnZraF8k5ldkqlNMysPbbwS2iwM5deaWb2ZrQ2v69/JwEVEJHoZQ8bMYsCdwKVAJXC1mVV2qXYd0OjuM4A7gNvDupX6duj7AAAIIElEQVTAImA2MB+4y8xiGdq8HbjD3SuAxtD2ESvc/YzwuveYRiwiIn0mmyOZuUCtu7/q7m3AcmBBlzoLgGVh+mHgIjOzUL7c3VvdfQtQG9pL2WZY5/2hDUKbHz724YmISC5lEzKTgR1J83WhLGUdd48DTcDYbtZNVz4W2BfaSLWtj5nZOjN72MymZNF3ERHJoWxCxlKUeZZ1eqsc4GfAdHefAzzGm0dOb+2I2Q1mVmNmNfX19amqiIhIH8kmZOqA5KOGMmBnujpmlg+MBBq6WTdd+R5gVGjjLdty973u3hrK/xM4O1Vn3X2Ju1e5e1VpaWkWwxMRkahkEzLPAhXhrq9CEhfyq7vUqQYWh+mFwOPu7qF8Ubj7rByoAFanazOs80Rog9DmTwHM7ISk7V0JvNSzoYqISF/Lz1TB3eNmdjPwKBAD7nP39WZ2G1Dj7tXAUuABM6slcQSzKKy73sweAjYAceAmd+8ASNVm2OTfAcvN7F+ANaFtgFvM7MrQTgNw7TsevYiIRMoSBw+DU1VVldfU1OS6GyIiA4qZPefuVb3Rlr7xLyIikVHIiIhIZBQyIiISGYWMiIhERiEjIiKRUciIiEhkFDIiIhIZhYyIiERGISMiIpFRyIiISGQUMiIiEhmFjIiIREYhIyIikVHIiIhIZBQyIiISGYWMiIhERiEjIiKRUciIiEhkFDIiIhKZrELGzOab2SYzqzWzL6VYXmRmK8LyZ8xsetKyW0P5JjO7JFObZlYe2ngltFmYaRsiItI/ZQwZM4sBdwKXApXA1WZW2aXadUCju88A7gBuD+tWAouA2cB84C4zi2Vo83bgDnevABpD22m3ISIi/Vc2RzJzgVp3f9Xd24DlwIIudRYAy8L0w8BFZmahfLm7t7r7FqA2tJeyzbDO+0MbhDY/nGEbIiLST+VnUWcysCNpvg54V7o67h43syZgbCh/usu6k8N0qjbHAvvcPZ6ifrpt7EnuiJndANwQZlvN7MUsxjhQjaPL+AcZjW/gGsxjg8E/vlN6q6FsQibV0YJnWSddeaojqO7qZ9sP3H0JsATAzGrcvSrFeoOCxjewDebxDeaxwfExvt5qK5vTZXXAlKT5MmBnujpmlg+MBBq6WTdd+R5gVGij67bSbUNERPqpbELmWaAi3PVVSOJCfnWXOtXA4jC9EHjc3T2ULwp3hpUDFcDqdG2GdZ4IbRDa/GmGbYiISD+V8XRZuP5xM/AoEAPuc/f1ZnYbUOPu1cBS4AEzqyVxdLEorLvezB4CNgBx4CZ37wBI1WbY5N8By83sX4A1oW3SbSODJVnUGcg0voFtMI9vMI8NNL6smQ4GREQkKvrGv4iIREYhIyIikRlQIWNmxWa22sxeMLP1Zva1UJ7uUTRTzewJM1tjZuvM7LKktlI+7iaXemt8ZjbdzA6b2drwuieX4zriGMY3zcx+G8b2pJmVJbW1ONR/xcwWp9tmX+rl8XUk7b+uN9rkRDfjuzn8X3IzG5dU38zsO2HZOjM7K2lZv9p/vTy2wbDvZprZKjNrNbMvdGmr28eMvY27D5gXie/KDA/TBcAzwLnAQ8CiUH4PcGOYXpI0XQlsTZp+ASgCyoHNQGwQjW868GKux9ML4/sxsDhMvx94IEyPAV4N76PD9OjBMr4wfzDX4+nB+M4Mn7mtwLik+pcBvwzrnQs801/3X2+NbRDtu/HAOcDXgS8klcdI/Lw8ESgk8XO0srttD6gjGU84GGYLwstJ/ygaB0aE6ZG8+Z2bdI+7yaleHF+/dAzjqwR+G6af4M3HGV0C/MbdG9y9EfgNiWfj5VQvjq9fSjc+d1/j7ltTrLIAuD+s9zSJ78CdQD/cf704tn6pp+Nz993u/izQ3mVRNo8Ze4sBFTKQeGCnma0FdpP4cG4m/aNovgpcY2Z1wErgr0N5qkflTKYf6KXxAZSH02i/M7N3903vM+vh+F4APhamPwKUmNmRxxUNhv2XbnwAxWZWY2ZPm9mRUMq5ruNz92e6qZ5uP/XL/ddLY4PBse/S6fG+G3Ah4+4d7n4GiacBzAVmpaoW3q8G/q+7l5E4vH3AzPLI8hE1udBL49sFTHX3M4HPAz8ysxEp2ulzPRzfF4D3mtka4L3AayS+bzVY9l+68UFi/1UBnwD+3cxOirbn2ek6PjM7tZvqPX3cVE710thgcOy7dHq87wZcyBzh7vuAJ0mcV0z3KJrrSJwPx91XAcUkHmyXzaNycuqdjC+cBtwbyp8j8dv0yX3X+8yyGZ+773T3j4aw/IdQ1sQg2X/djA93P1Ln1dDOmX3Z/0ySxtfdaa6ePlaqX3iHYxss+y6dHu+7ARUyZlZqZqPC9BDgYuAl0j+KZjtwUag/i8QP4XrSP+4mp3prfKGdWCg/kcT4Xu2rcaTT0/GZ2bhwZAZwK3BfmH4U+KCZjTaz0cAHQ1lO9db4wriKjtQBLiDx1IycSjO+jd2sUg18KtyJdS7Q5O676If7r7fGNoj2XTrZPGbsrbq7K6C/vYA5JB41sw54EfhyKD+RREjUkrhjpyiUVwJ/IHHuey3wwaS2/oHEb/ibgEtzPbbeHB+J8/zrQ/nzwBW5Htsxjm8h8ArwMnDvkfKw7DOhfi3w6VyPrTfHB5wP/Cnsvz8B1+V6bBnGdwuJ33DjJH6rvTeUG4k/Trg5jKOqv+6/3hrbINp3E0P5fmBfmB4Rll0WPrObgX/ItG09VkZERCIzoE6XiYjIwKKQERGRyChkREQkMgoZERGJjEJGREQio5AREZHIKGRERCQy/x8EJ308JJET5QAAAABJRU5ErkJggg==\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": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"211.4\n",
"4781.599999999999\n",
"807.8926831881124\n",
"Time taken: 2 s\n"
]
}
],
"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": 33,
"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": 34,
"metadata": {},
"outputs": [],
"source": [
"# print(36000*(1+ pdg[\"jpsi_BR\"]/pdg[\"NR_BR\"] + pdg[\"psi2s_BR\"]/pdg[\"NR_BR\"]))"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"nan nan %\n",
"Time: 1000.0239598751068\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"c:\\users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\ipykernel_launcher.py:3: IntegrationWarning: The maximum number of subdivisions (5) has been achieved.\n",
" If increasing the limit yields no improvement it is advised to analyze \n",
" the integrand in order to determine the difficulties. If the position of a \n",
" local difficulty can be determined (singularity, discontinuity) one will \n",
" probably gain from splitting up the interval and calling the integrator \n",
" on the subranges. Perhaps a special-purpose integrator should be used.\n",
" This is separate from the ipykernel package so we can avoid doing imports until\n"
]
}
],
"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": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAD8CAYAAABZ/vJZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XmUVPWd9/H3t6o3tmYVRRqhxVbAJS7tEsWEqCNuEWNMJNGEqNHn5DEzSWaSGR0zj5lJnDEnmZgxJjoojsQxoONKEhNCiE40KgqCyiI7YrMINtACDb1UfZ8/6jYW2Out6r5V1Z/XOX3q3t/93dvf37nd/em71C1zd0RERLoqFnUBIiKSnxQgIiISigJERERCUYCIiEgoChAREQlFASIiIqF0GCBm9qCZbTOzpWltPzKzt83sTTN7yswGpS271czWmNlKM5uc1n5R0LbGzG5Ja680swVmttrMHjWzkmwOUEREukdnjkAeAi46pG0ecIK7nwSsAm4FMLMJwFTg+GCdX5hZ3MziwM+Bi4EJwBeCvgA/BO5y9ypgJ3BDRiMSEZEe0WGAuPufgR2HtP3B3ZuD2VeAimB6CjDb3RvcfT2wBjgj+Frj7uvcvRGYDUwxMwPOAx4P1p8JXJHhmEREpAcUZWEb1wOPBtMjSQVKi5qgDeDdQ9rPBIYCu9LCKL3/R5jZTcBNAP369Ttt3LhxGRcvheutTXVZ29YA6hljW1nrR1JPWda2my0n2Hre94FsZUjG2zpx5MAsVCS5aNiwYcydO3euux96VimUjALEzG4DmoFHWppa6ea0fqTj7fRvlbtPB6YDVFdX+8KFC7tUr/QuY275bda2dW18Hj8o/i/O2P9vbGNw1rabLb8p/TovJE7k75v/T8bbWnjnpVmoSHKVmQ3L1rZCB4iZTQMuA873Dx+oVQOMSutWAWwOpltrfx8YZGZFwVFIen+RnHGk1dLocbaTm/+dv+eDONx2Rl2G9DKhbuM1s4uAfwAud/f6tEVzgKlmVmpmlUAV8CrwGlAV3HFVQupC+5wgeJ4DrgrWnwY8E24oIt1npL3PVh+C5+id71t9KEfYjo47imRRZ27jnQW8DBxnZjVmdgNwDzAAmGdmS8zsPgB3XwY8BiwHfg/c7O6J4Oji68BcYAXwWNAXUkH0t2a2htQ1kRlZHaFIFhxl29jow6Muo01bfAgjFCDSwzo8heXuX2iluc0/8u5+B3BHK+3PAs+20r6O1F1aIjlrtG3ld8kzoy6jTe/5YMqtnn7sYy99oi4npzU1NVFTU8P+/fujLqVblZWVUVFRQXFxcbd9j2zchSVS0MrZwxDbwwY/POpS2rTFU3dfHWE7WOtt3sgoQE1NDQMGDGDMmDGk3klQeNyd2tpaampqqKys7Lbvk5sndEVyyGjbBsA7ORwgW30ogC6kd8L+/fsZOnRowYYHgJkxdOjQbj/KUoCIdGCMbQVyO0C2BO//GIGug3RGIYdHi54YowJEpANHBUcguXwR/T1PvTdFd2JJT1KAiHRgjG1lqw9mXw6+A71FAyXU+gBGWG3UpUiWPP/881x22WUANDQ0cMEFF3DyySfz6KOPdrBmz9FFdJEOjI69l9Onr1ps9SE6AilQixcvpqmpiSVLlkRdykF0BCLSgUrbyjvJ3A8QvRckf2zYsIFx48Yxbdo0TjrpJK666irq6+v5/e9/z7hx45g4cSJPPvkkANu2bePaa69lyZIlnHzyyaxduzbi6j+kIxCRdpSzh8OsjjV+ZNSldGizD+P02Mqoy8gvv7sFtr6V3W0ecSJcfGeH3VauXMmMGTM455xzuP766/nJT37Cf/7nf/KnP/2JY445hquvvhqA4cOH88ADD/DjH/+Y3/zmN9mtNUM6AhFpR5VtAmC1V3TQM3rv+mEMtHrK2Rt1KdIJo0aN4pxzzgHg2muvZeHChVRWVlJVVYWZce2110ZcYcd0BCLSjmNjNUB+BEjLXWKjbDvLvF/E1eSJThwpdJdDb7Otq6vLu9uLdQQi0o4q20S9l7IpeKNeLnv3QIBsi7gS6YyNGzfy8ssvAzBr1iwuuOAC1q9ff+Aax6xZs6Isr1MUICLtqLIaVvvInH0Kb7qWADnK3ou4EumM8ePHM3PmTE466SR27NjBt771LaZPn86ll17KxIkTGT16dNQldkinsETaURXbxF+SJ0RdRqfspi87vf+BNz5KbovFYtx3330HtV100UW8/fbbH+k7adIkJk2a1EOVdV7u/1slEpFy9nKE7WRVMvevf7TY6MMVINJjFCAibaiylgvo+fN023d9OBW2PeoypANjxoxh6dKlUZeRMQWISBuOj20AYFlyTKR1dMXGIEBiJKMuRXoBBYhIG06wDWz3ct5jcNSldNpGH06JJRiBnokl3U8BItKGE2IbWJasBPLn3vx1yREAjI1tjrgS6Q0UICKtKKWRKqthqY+JupQuWRs8cuUYU4BI99NtvCKtONZqKLYES5Pd93Gg3aGWcnZ6f8YqQDptzC2/zer2Ntx5aZfX+d73vkf//v359re/3eryp59+mmOPPZYJEyZkWl5W6QhEpBUnxNYD5N0RCBhr/EiOiW2KuhDJoqeffprly5dHXcZHKEBEWnGiraPO+1Ljh0VdSpetTR6pI5A8cMcdd3DcccdxwQUXsHJl6inK999/P6effjof+9jH+OxnP0t9fT0vvfQSc+bM4Tvf+c6Bx7m31i8KChCRVpwWW83rySry6QJ6izU+kmH2AYPYHXUp0oZFixYxe/ZsFi9ezJNPPslrr70GwJVXXslrr73GG2+8wfjx45kxYwZnn302l19+OT/60Y9YsmQJY8eObbVfFHQNROQQ5ezhuFgNv276eNSlhNLy2SVjbTOL/LiIq5HWvPDCC3zmM5+hb9++AFx++eUALF26lO9+97vs2rWLPXv2MHny5FbX72y/7qYjEJFDnBpbDcAiPzbiSsJZE7xzvkrXQXJaa49u/8pXvsI999zDW2+9xe23387+/ftbXbez/bqbAkTkEKfFVtPsMZYkx0ZdSiibfBi7vQ/j7Z2oS5E2fOITn+Cpp55i37597N69m1//+tcA7N69mxEjRtDU1MQjjzxyoP+AAQPYvfvDU5Jt9etpOoUlcohqW8VyH80+yqIuJRQnxgo/igkxBUhnhLntNlOnnnoqV199NSeffDKjR4/m3HPPBeD73/8+Z555JqNHj+bEE088EBpTp07lxhtv5O677+bxxx9vs19PM3dvv4PZg8BlwDZ3PyFoGwI8CowBNgCfd/edljom+w/gEqAe+Iq7vx6sMw34brDZH7j7zKD9NOAhoA/wLPAN76gooLq62hcuXNiVsUovE+b+/hKaeKP0RmYnPsU/N0/rhqp6xu1FM/lc/H85seGBLn+WSRR/UHvSihUrGD9+fNRl9IjWxmpmi9y9Ohvb78xP1kPARYe03QLMd/cqYH4wD3AxUBV83QTcGxQ8BLgdOBM4A7jdzFoeMHRv0LdlvUO/l0iPOTW2mj7WmDefAdKW5T6a/rafMfpwKelGHQaIu/8Z2HFI8xRgZjA9E7girf2XnvIKMMjMRgCTgXnuvsPddwLzgIuCZeXu/nJw1PHLtG2J9LhzYktp9hgLkvn9H+ry4AnCx9uGSOuQwhb2Ivrh7r4FIHgdHrSPBN5N61cTtLXXXtNKu0gkzo29xRI/ht30jbqUjKzyCho9fuCR9HKwTpwlz3s9McZs34XV2ruuPER76xs3u8nMFprZwu3b9aE5kl3l7OFEW5f3p68AmihitVfoCKQVZWVl1NbWFnSIuDu1tbWUlXXvjSBh78J6z8xGuPuW4DRUy2do1gCj0vpVAJuD9kmHtD8ftFe00r9V7j4dmA6pi+ghaxdp1cdjy4mb80Ii/wMEYEnyGD4dfxkj2eUL6YWsoqKCmpoaCv2f0LKyMioquvfjmMMGyBxgGnBn8PpMWvvXzWw2qQvmdUHIzAX+Ne3C+YXAre6+w8x2m9lZwALgy8DPQtYkkpG/ir9OnfdliR8TdSlZsShZxTVF86myTazyUR2v0EsUFxdTWZlfT1nOVR0GiJnNInX0MMzMakjdTXUn8JiZ3QBsBD4XdH+W1C28a0jdxnsdQBAU3wdeC/r9i7u3XJj/Gh/exvu74EukR8VJcF7sdeYnT6W5QN4e1fJO+tNiq1iVUIBI9nX4m+LuX2hj0fmt9HXg5ja28yDwYCvtC4HCOGcgeavaVjHE9vCHRFZuj88J7/jhvO/lnBZbzazER35dRTKmE6MiwIXxhTR4MX9OnhR1KVlkvJ6s4lRbFXUhUqAUINLrGUkmx1/jxeQJ1Ofp40vasih5LEfHtjKUuqhLkQKkAJFer9pWUWHv8+tEfj6+vT0tb4g8O7Ys4kqkEClApNf7TPwF9nopc5OFc/2jxZt+NHXel4mxpVGXIgVIASK9WimNXBZfwNzk6Xn79N32JInxUvJ4Jsbfop336IqEogCRXu282GLKrZ6nEhOjLqXbvJg8kZFWS6VtjboUKTAKEOnVro3/kRofVhCPL2nLC8kTATg39mbElUihUYBIr3WM1XBOfBm/aj6fZAH/Kmz0w1mfPJzzY4ujLkUKTOH+1oh04EvxeTR4EbMTn4q6lG43N3k6H48to5y9UZciBUQBIr1SOXu5Mv4iv0mexQ7Koy6n281NnE6JJfiUjkIkixQg0itNi89lgO3jgebC/vjWFkt8LFt9MJPj+hhoyR4FiPQ6/djH9UW/Z17iVFb46KjL6RFOjD8kqpkUe4O+7I+6HCkQChDpdb4Un8dg28M9zb3r05OfSZxNX2vg4tirUZciBUIBIr3KIHbztaI5PJ/4GG8UyOd+dNYiP5b1ycO5Kv7nqEuRAqEAkV7lb4qeoj/7+NfmL0ZdSgSMxxOf5OPx5VTYto67i3RAASK9RqVt4UvxeTyamNRrP6HvqcREkm58Pv581KVIAVCASC/h3FE0g32UcFfz5zruXqA2M4z5yVO4Jj6fUhqjLkfynAJEeoWr489zdnw5/9b8RbYzKOpyIvVg4mKG2m6uiP8l6lIkzylApOAdQS23FT3CK8nxveJd5x15OTmB5cnRXB//HXpCr2RCASIFLU6Cu0vuIU6CW5q+iutHHjDub76E42I1TI7pjYUSnn6bpKB9s+gJzoit5B+bbmCDj4i6nJwxJ3k2a5Mj+Nui/yFGMupyJE8pQKRgTY69xs3xZ5jdPIlnkoX7eR9hJIhzV/NVHBer4dOxl6IuR/KUAkQKU80iflr8c5b4WG5v/krU1eSk3ybPZHlyNN8pfowyGqIuR/KQAkQKT+1amDWV7T6QGxv/jgZKoq4oJzkxvtf0ZSrsfb5e9HTU5UgeUoBIYaldCw9dCp7kuqa/p5aBUVeU01718TyRmMhN8d8w1jZFXY7kGQWIFI7tq+ChyyDRCNPmsNZHRl1RXvi3pmvYSx9+UnwvRTRHXY7kEQWIFIYNf4EZF0CyCb48Bw4/PuqK8sb7DOTWpq/ysdg6vln0RNTlSB5RgEj+W/IrePgK6DccvvpHOOKEqCvKO79PnsFjzZ/k/8bnwOo/Rl2O5ImMAsTMvmVmy8xsqZnNMrMyM6s0swVmttrMHjWzkqBvaTC/Jlg+Jm07twbtK81scmZDkl6jsR6euRme/hqMOhNu+AMMHhN1VXnr9uZpvO1HwePXpU4HinQgdICY2Ujgb4Bqdz8BiANTgR8Cd7l7FbATuCFY5QZgp7sfA9wV9MPMJgTrHQ9cBPzCzOJh65JeYtPrcP95sPgR+MR34EtPQ98hUVeV1/ZRxlcb/w6KSuFXn4fdW6MuSXJcpqewioA+ZlYE9AW2AOcBjwfLZwItH/s2JZgnWH6+mVnQPtvdG9x9PbAGOCPDuqRQNdbD3NvggfNh/y649gk477sQL4q6soKwmWEwdRbs2Qa/nAJ734+6JMlhoQPE3TcBPwY2kgqOOmARsMvdW27lqAFaboUZCbwbrNsc9B+a3t7KOgcxs5vMbKGZLdy+fXvY0iUfJZPwxqPw8zPg5Xvg1Glw8wI45vyoKys8o06Hax6Dne/AzMvhgy1RVyQ5KpNTWINJHT1UAkcC/YCLW+na8rhPa2NZW+0fbXSf7u7V7l592GGHdb1oyT/usGY+3D8JnropdZrqK8/Cp38KZXqPR7cZMxG+MAt2vQMPXADbVkRdkeSgTE5hXQCsd/ft7t4EPAmcDQwKTmkBVACbg+kaYBRAsHwgsCO9vZV1pLdKJmDZ0zB9Evz3lVC/A668H258HsacE3V1vcPYT8F1z6ZujZ5xISx/JuqKJMdkEiAbgbPMrG9wLeN8YDnwHHBV0Gca0PJTNyeYJ1j+J3f3oH1qcJdWJVAFvJpBXZLP9myDF++Cn50G/zMNGvfA5T+Dv14EJ30eYrrzvEeN+Bh8dT4MPQYe+zI8+x1o2hd1VZIjQl95dPcFZvY48DrQDCwGpgO/BWab2Q+CthnBKjOAh81sDakjj6nBdpaZ2WOkwqcZuNndE2HrkjzU3ABr/5R6P8fKZyHZDKMnwgXfg/GfhphuyovUoFFw/Vz44+3wyi9gzR/hsp/C0Z+MujKJmKUOAvJPdXW1L1yoD8PJW037Yd1zsOwpWPk7aPgA+gyBk7+YukB+2LEZf4sxt/w2C4X2PhvuvLTtheueh19/E3auhxM/B5+6DYZU9lhtkjkzW+Tu1VnZlgIk9xXGH0JnrG1mUuwNPhF7kzNjKyizJnZ5P/6QqOa3ybN4KXk8TeEPiqWHlNLIzUVPc2P8WeIkmJU4j/uaL2cLQ6MuLXLthm+OyGaA6LdVukWcBOPsXU6LraQ6torq2EqOtB0ArE2OYFbiPJ5PnqzQyEMNlPCT5s/z381/xTeKnuSa+Hyuic/nt8mzmNF8MW/62KhLlB6i31zJWJwER9sWxts7TIi9wwm2npNja+lv+wHY4kNYlDyWnyWP58+Jk9iEbsEuBNsYzG3NN3Bv4nKmxedydfw5ppS+xPLkaJ5ITGRO4hy2MyjqMqUbKUCk04po5ijbRqVtodK2coxtYnxsI+PsXUqtCYAGL2KVV/BE4lwWJY9jUbKKTQyj9bf7SCGo8cO4o/la/qP5Sj4Tf5HPxl/gn4of4daiWbyaHMf85CnMT56qz6QvQAoQSeMMZC8j7X1G2vscabWMsm1U2lYqbQtH2TaKLHmgd60PYEXyKGYmL2RF8iiW+2jW+pE068eqV9pDXx5OXMjDiQsZa5u4Iv4X/iq2iH8qfoR/4hHWJkfwSnICC5LjeDU5jq26ZpL39JveS5TSyDDqOMzqGGZ1HGa7OIw6jrAdHGm1HBkERstppxb7vZj1PoIVfhTPJs9kXXIE630E63wEdfSPaDSS69b6SP69+fP8O5+nwrZxfmwxk2JL+HT8Ja4pmg/AxuRhvOlHsyxZyTIfzbLkGH2CZJ5RgOShOAkGspdBtodB7El73cvAYHqo7U4FBbs4zOoot/pWt1XrA9jkw1jnR/Ji8kQ2+TA2+TA2+1A2+TBqKUennyQTNT6cmYnJzExMJkaS8fYOZ8bepjq2kpNsHZcVLzjQ9z0fxLrkkaz3I1jrH/6zstmH0UhxhKPonMK4Y7Lz8vY23tIRVT5i2k+jLqOLnL400J999Ld9B7+yj362nwHso99By/bTn32U214GkQqIcmv7ncBJN+rox07vz3YGsd0H8r4PZLsPYjtp0z6QWgbqDiiJXDl7OD72DhNsA+PsXSpjWznaNjPE9hzUb7sPZLMPZYsPZXPw9Z4PppZyar2cHV7ODgaQQG88bc87P7xMt/FmS4wkJTRRQhOlNFNqTcF8MyU00YdG+lgDZTTShwb6WOq17KD2RsqsIdU36FPGofON9GM/Mes4sJs8zh76sMf7pF4pY7sPYhUV1CX7s8v7s4t+7PL+1NGfXd6PnQxgl/djN31xfdCk5JEP6M/LyeN5mYM/hngQuzk6uGFjpL3PCKvlSKtlrG1mYuytj5xubbHD+7PDy6mlnLrgd+ID78sH9GW39z0wv5u+7A5+x/Z5KfsooZ4yGihGR92dk7cBcjg7uaXoV5Sm/bEvseYDf/xLaTwwX8qH7SXWdNA66ReFuyrhRj1l7KeEfV7CPlI/hA2UsMsHsIVU2/5kCfspYQ9lQSj0ZY+XsYc+7A1+gHfTh71Bm36ARWAXA3jdB/C6t/ZUAqeceobbToaym6FWx5DgtO2QYH6YfcAo284Aq6ecevqzr1P/wAHUB4Gyj1L2ewn1lB6YbqSYBopooogmT702Bl8tbY0U00g8NR98NXoxzcRIECNJjGbiJIP5BDES3s6ytK+kpy83HCOJ4cRwIEksmE9f9uFrNuVtgAy3nVwXn0sDLTsrtYMaWqYpptGLqaM02OGpnd6YTFtOEQ1ecvA8xTR6EQ2k2lt+aFrCYZ+XpgKDUpqIoz/0IlEwPqAfH3g/1kAbHwBx6BpJ+pE6TVxuexlAPQNsHwPYR2lwZqHvIWcQ+lrDQWcTBtpeSmimOPgqiTWlXmmmmETG/5T2hGz+xcrbAHnLj+a4hny7BiIiUXFiqaN/+rLFh7Y0Zl2M5IGAORAu1hz8y5kMjh2SFJE4MB3HiVviI8tj+IG2OAni1tL/w6+WY43YgWOO1PTB80kMiFkSeChrY83bABERyUVJYjQEp7IP6ExQ9dj9TA9lbUu62ioiIqEoQEREJBQFiIiIhKIAERGRUBQgIiISigJERERCUYCIiEgoChAREQlFASIiIqEoQEREJBQFiIiIhKIAERGRUBQgIiISigJERERCyShAzGyQmT1uZm+b2Qoz+7iZDTGzeWa2OngdHPQ1M7vbzNaY2ZtmdmradqYF/Veb2bRMByUiIt0v0yOQ/wB+7+7jgI8BK4BbgPnuXgXMD+YBLgaqgq+bgHsBzGwIcDtwJnAGcHtL6IiISO4KHSBmVg58ApgB4O6N7r4LmALMDLrNBK4IpqcAv/SUV4BBZjYCmAzMc/cd7r4TmAdcFLYuERHpGZkcgRwNbAf+y8wWm9kDZtYPONzdtwAEr8OD/iOBd9PWrwna2mr/CDO7ycwWmtnCRH1dBqWLiEimMgmQIuBU4F53PwXYy4enq1rT2me5ezvtH210n+7u1e5eHe87sKv1iohIFmUSIDVAjbsvCOYfJxUo7wWnpghet6X1H5W2fgWwuZ12ERHJYaEDxN23Au+a2XFB0/nAcmAO0HIn1TTgmWB6DvDl4G6ss4C64BTXXOBCMxscXDy/MGgTEZEcVpTh+n8NPGJmJcA64DpSofSYmd0AbAQ+F/R9FrgEWAPUB31x9x1m9n3gtaDfv7j7jgzrEhGRbpZRgLj7EqC6lUXnt9LXgZvb2M6DwIOZ1CIiIj1L70QXEZFQFCAiIhKKAkREREJRgIiISCgKEBERCUUBIiIioShAREQkFAWIiIiEogAREZFQFCAiIhKKAkREREJRgIiISCgKEBERCUUBIiIioShAREQkFAWIiIiEogAREZFQFCAiIhKKAkREREJRgIiISCgKEBERCUUBIiIioShAREQkFAWIiIiEogAREZFQFCAiIhKKAkRERELJOEDMLG5mi83sN8F8pZktMLPVZvaomZUE7aXB/Jpg+Zi0bdwatK80s8mZ1iQiIt0vG0cg3wBWpM3/ELjL3auAncANQfsNwE53Pwa4K+iHmU0ApgLHAxcBvzCzeBbqEhGRbpRRgJhZBXAp8EAwb8B5wONBl5nAFcH0lGCeYPn5Qf8pwGx3b3D39cAa4IxM6hIRke6X6RHIT4G/B5LB/FBgl7s3B/M1wMhgeiTwLkCwvC7of6C9lXUOYmY3mdlCM1uYqK/LsHQREclE6AAxs8uAbe6+KL25la7ewbL21jm40X26u1e7e3W878Au1SsiItlVlMG65wCXm9klQBlQTuqIZJCZFQVHGRXA5qB/DTAKqDGzImAgsCOtvUX6OiIikqNCH4G4+63uXuHuY0hdBP+Tu18DPAdcFXSbBjwTTM8J5gmW/8ndPWifGtylVQlUAa+GrUtERHpGJkcgbfkHYLaZ/QBYDMwI2mcAD5vZGlJHHlMB3H2ZmT0GLAeagZvdPdENdYmISBZlJUDc/Xng+WB6Ha3cReXu+4HPtbH+HcAd2ahFRER6ht6JLiIioShAREQkFAWIiIiEogAREZFQFCAiIhKKAkREREJRgIiISCgKEBERCUUBIiIioShAREQkFAWIiIiEogAREZFQFCAiIhKKAkREREJRgIiISCgKEBERCUUBIiIioShAREQkFAWIiIiEogAREZFQFCAiIhKKAkREREJRgIiISCgKEBERCUUBIiIioShAREQkFAWIiIiEEjpAzGyUmT1nZivMbJmZfSNoH2Jm88xsdfA6OGg3M7vbzNaY2ZtmdmratqYF/Veb2bTMhyUiIt0tkyOQZuDv3H08cBZws5lNAG4B5rt7FTA/mAe4GKgKvm4C7oVU4AC3A2cCZwC3t4SOiIjkrtAB4u5b3P31YHo3sAIYCUwBZgbdZgJXBNNTgF96yivAIDMbAUwG5rn7DnffCcwDLgpbl4iI9IysXAMxszHAKcAC4HB33wKpkAGGB91GAu+mrVYTtLXVLiIiOSzjADGz/sATwDfd/YP2urbS5u20t/a9bjKzhWa2MFFf1/ViRUQkazIKEDMrJhUej7j7k0Hze8GpKYLXbUF7DTAqbfUKYHM77R/h7tPdvdrdq+N9B2ZSuoiIZCiTu7AMmAGscPefpC2aA7TcSTUNeCat/cvB3VhnAXXBKa65wIVmNji4eH5h0CYiIjmsKIN1zwG+BLxlZkuCtn8E7gQeM7MbgI3A54JlzwKXAGuAeuA6AHffYWbfB14L+v2Lu+/IoC4REekBoQPE3V+k9esXAOe30t+Bm9vY1oPAg2FrERGRnqd3oouISCgKEBERCUUBIiIioShAREQkFAWIiIiEogAREZFQFCAiIhKKAkREREJRgIiISCgKEBERCUUBIiIioShAREQkFAWIiIiEogAREZFQFCAiIhKKAkREREJRgIiISCgKEBERCUUBIiIioShAREQkFAWIiIiEogAREZFQFCAiIhKKAkREREJRgIiISCgKEBERCUUBIiIioShAREQklJwJEDPeQ/IjAAAFW0lEQVS7yMxWmtkaM7sl6npERKR9OREgZhYHfg5cDEwAvmBmE6KtSkRE2pMTAQKcAaxx93Xu3gjMBqZEXJOIiLSjKOoCAiOBd9Pma4AzD+1kZjcBNwWzDe/88LKlPVBbFIYB70ddRDfS+PKbxpe/ysji2HIlQKyVNv9Ig/t0YDqAmS109+ruLiwKhTw20PjyncaXv4KxXZSt7eXKKawaYFTafAWwOaJaRESkE3IlQF4Dqsys0sxKgKnAnIhrEhGRduTEKSx3bzazrwNzgTjwoLsv62C16d1fWWQKeWyg8eU7jS9/ZXVs5v6RSw0iIiIdypVTWCIikmcUICIiEkrOBIiZlZnZq2b2hpktM7N/DtorzWyBma02s0eDi+yY2VFm9pyZLTazN83skrRt3Ro8EmWlmU2OakzpsjU+MxtjZvvMbEnwdV+U4wpq6urYRpvZ/GBcz5tZRdq2pgX9V5vZtKjGlC7L40uk7bucuFGknfF9Pfg9cjMbltbfzOzuYNmbZnZq2rJ82n9hxlcI+2+cmb1sZg1m9u1DttW1R0q5e058kXovSP9guhhYAJwFPAZMDdrvA74WTE9Pm54AbEibfgMoBSqBtUC8gMY3Blga9XgyHNv/ANOC6fOAh4PpIcC64HVwMD24UMYXzO+JejxdGN8pwc/bBmBYWv9LgN8F650FLMjT/del8RXQ/hsOnA7cAXw7rT1O6u/l0UAJqb+jE9r73jlzBOIpe4LZ4uDLSf0CPh60zwSuaFkFKA+mB/Lh+0amALPdvcHd1wNrSD0qJVJZHF/OCTG2CcD8YPo5PnxszWRgnrvvcPedwDwga296CiuL48tJbY3P3Re7+4ZWVpkC/DJY7xVgkJmNIM/2X4jx5aSujs/dt7n7a0DTIYu6/EipnAkQSD1U0cyWANtI/fCtBXa5e3PQpYbUY08Avgdca2Y1wLPAXwftrT0WZSQ5IEvjA6gMTm39r5md2zPVt6+LY3sD+Gww/RlggJkNpXD2XVvjAygzs4Vm9oqZtQRO5A4dn7svaKd7W/spb/ZfyPFBYey/tnR5/+VUgLh7wt1PJvVO9DOA8a11C16/ADzk7hWkDjkfNrMYnXwsShSyNL4twFHufgrwt8CvzKy8le30qC6O7dvAJ81sMfBJYBPQTOHsu7bGB6l9Vw18EfipmY3t3so759DxmdkJ7XRvaz/lzf4LOT4ojP3Xli7vv5wKkBbuvgt4ntR5vEFm1vKGx/RHnNxA6hw07v4yqYeEDSMPHouSyfiCU3O1QfsiUv8JH9tz1bevM2Nz983ufmUQgrcFbXUUyL5rZ3y4e0ufdcF2TunJ+juSNr72Tj21tZ/yaf+FGV+h7L+2dHn/5UyAmNlhZjYomO4DXACsIHUO+aqg2zTgmWB6I3B+0H88qT+w20k9AmWqmZWaWSVQBbzaU+NoS7bGF2wnHrQfTWp863pqHK3p6tjMbFhwNAVwK/BgMD0XuNDMBpvZYODCoC1S2RpfMK7Slj7AOcDynhpHW9oY39vtrDIH+HJwt9JZQJ27byG/9l+Xx1dA+68tXX+kVHtX2HvyCzgJWAy8CSwF/l/QfjSpAFhD6u6W0qB9AvAXUueblwAXpm3rNlL/ma8ELo56bNkcH6lz68uC9teBT+fh2K4CVgOrgAda2oNl1wf91wDXRT22bI4POBt4K9h3bwE3RD22Dsb3N6T+K20m9Z/oA0G7kfoAuLXBOKrzdP91aXwFtP+OCNo/AHYF0+XBskuCn9u1wG0dfW89ykRERELJmVNYIiKSXxQgIiISigJERERCUYCIiEgoChAREQlFASIiIqEoQEREJJT/D3KEMKnbI7XxAAAAAElFTkSuQmCC\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 = 23448660.25141733</td>\n",
" <td title=\"Total number of call to FCN so far\">TOTAL NCALL = 73</td>\n",
" <td title=\"Number of call in last migrad\">NCALLS = 73</td>\n",
" </tr>\n",
" <tr>\n",
" <td title=\"Estimated distance to minimum\">EDM = nan</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:#FF7878\">False</td>\n",
" <td align=\"center\" style=\"background-color:#FF7878\">False</td>\n",
" <td align=\"center\" style=\"background-color:#FF7878\">True</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=\"$('#mfLVoBAVNv').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_p</td>\n",
" <td>-1.5</td>\n",
" <td>nan</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_p</td>\n",
" <td>-1.5</td>\n",
" <td>nan</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_s</td>\n",
" <td>23.07</td>\n",
" <td>nan</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>184.39</td>\n",
" <td>nan</td>\n",
" <td></td>\n",
" <td></td>\n",
" <td></td>\n",
" <td></td>\n",
" <td>No</td>\n",
" </tr>\n",
"</table>\n",
"<pre id=\"mfLVoBAVNv\" 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_{p}$ & -1.5 & nan & & & & & No\\\\\n",
"\\hline\n",
"1 & $psi2s_{p}$ & -1.5 & nan & & & & & No\\\\\n",
"\\hline\n",
"2 & $psi2s_{s}$ & 23.07 & nan & & & & & No\\\\\n",
"\\hline\n",
"3 & $jpsi_{s}$ & 184.39 & nan & & & & & 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": "exception was raised in user function\nUser function arguments:\n jpsi_p = -nan(ind)\n psi2s_p = +nan\n psi2s_s = +nan\n jpsi_s = +nan\nOriginal python exception in user function:\nKeyboardInterrupt: \n File \"c:\\users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\minimizers\\minimizer_minuit.py\", line 48, in func\n loss_evaluated = self.sess.run(loss_val)\n File \"c:\\users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\client\\session.py\", line 929, in run\n run_metadata_ptr)\n File \"c:\\users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\client\\session.py\", line 1152, in _run\n feed_dict_tensor, options, run_metadata)\n File \"c:\\users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\client\\session.py\", line 1328, in _do_run\n run_metadata)\n File \"c:\\users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\client\\session.py\", line 1334, in _do_call\n return fn(*args)\n File \"c:\\users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\client\\session.py\", line 1319, in _run_fn\n options, feed_dict, fetch_list, target_list, run_metadata)\n File \"c:\\users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\client\\session.py\", line 1407, in _call_tf_sessionrun\n run_metadata)\n",
"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: exception was raised in user function\nUser function arguments:\n jpsi_p = -nan(ind)\n psi2s_p = +nan\n psi2s_s = +nan\n jpsi_s = +nan\nOriginal python exception in user function:\nKeyboardInterrupt: \n File \"c:\\users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\minimizers\\minimizer_minuit.py\", line 48, in func\n loss_evaluated = self.sess.run(loss_val)\n File \"c:\\users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\client\\session.py\", line 929, in run\n run_metadata_ptr)\n File \"c:\\users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\client\\session.py\", line 1152, in _run\n feed_dict_tensor, options, run_metadata)\n File \"c:\\users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\client\\session.py\", line 1328, in _do_run\n run_metadata)\n File \"c:\\users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\client\\session.py\", line 1334, in _do_call\n return fn(*args)\n File \"c:\\users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\client\\session.py\", line 1319, in _run_fn\n options, feed_dict, fetch_list, target_list, run_metadata)\n File \"c:\\users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\client\\session.py\", line 1407, in _call_tf_sessionrun\n run_metadata)\n"
]
}
],
"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
}
