{
"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:57: UserWarning: Not running on Linux. Determining available cpus for thread can failand be overestimated. Workaround (only if too many cpus are used):`zfit.run.set_n_cpu(your_cpu_number)`\n",
" warnings.warn(\"Not running on Linux. Determining available cpus for thread can fail\"\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"WARNING: The TensorFlow contrib module will not be included in TensorFlow 2.0.\n",
"For more information, please see:\n",
" * https://github.com/tensorflow/community/blob/master/rfcs/20180907-contrib-sunset.md\n",
" * https://github.com/tensorflow/addons\n",
"If you depend on functionality not listed there, please file an issue.\n",
"\n"
]
}
],
"source": [
"import os\n",
"\n",
"# os.environ[\"CUDA_VISIBLE_DEVICES\"] = \"-1\"\n",
"\n",
"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\n",
"import tensorflow_probability as tfp\n",
"tfd = tfp.distributions"
]
},
{
"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, q) * _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, 200000)\n",
"\n",
"probs = total_f.pdf(test_q)\n",
"\n",
"calcs_test = zfit.run(probs)\n",
"res_y = zfit.run(jpsi_res(test_q))\n",
"f0_y = zfit.run(formfactor(test_q,\"0\"))\n",
"fplus_y = zfit.run(formfactor(test_q,\"+\"))\n",
"fT_y = zfit.run(formfactor(test_q,\"T\"))"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ8AAAD8CAYAAACo9anUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXwc1ZXo8d/pbu2SJWuzjSRbtiUvMtgGBBhswDEEDCExSWAwCQEmJAwTePMmeW8CJJONgbyQmUBeJiSBB2QIyWAIIcRhc1jMDrZlvGPLllfJm3ZZttbuvu+PKglZ6k2y1NVSn+/no49b1VX33ipLfXTvPXVLjDEopZRS0eRyugFKKaXijwYfpZRSUafBRymlVNRp8FFKKRV1GnyUUkpFnQYfpZRSURdR8BGRpSJSKSJVInJXgPeTRORp+/01IlLc57277e2VInJ5uDJFZKpdxi67zMQI6pgrIh+IyDYR2SIiyUO5GEoppaIjbPARETfwEHAFUAZcLyJl/Xa7BWgyxpQADwL328eWAcuBOcBS4Fci4g5T5v3Ag8aYUqDJLjtUHR7g98Btxpg5wGKge5DXQSmlVBRF0vM5F6gyxuwxxnQBK4Bl/fZZBjxhv34WuERExN6+whjTaYzZC1TZ5QUs0z5miV0GdplXh6njMmCzMWYTgDGmwRjji/wSKKWUijZPBPsUANV9vq8Bzgu2jzHGKyItQI69/cN+xxbYrwOVmQM0G2O8AfYPVscMwIjIKiAPK9j9tP9JiMitwK0AaWlpZ8+aNSuCU1cqfhkDWw+1MHFcMnkZSQPe7/L5qTzSSmFWCuPTEh1ooaW2tZOjxzo4/bRMRKxt1Y1ttHX5mDkxw7F2jUXr16+vN8bkDUdZkQQfCbCt/5o8wfYJtj1QjyvU/qHq8ACLgHOANuB1EVlvjHn9pB2NeQR4BKC8vNxUVFQEKE4p1eNEp5c5P1jFd66cxa0XTR/w/pGWDhb8n9e59wtncP25kx1ooeWXb+ziP/62kw/uvYJEj/XR8q1nNrJmTyPv3bXEsXaNRSKyf7jKimTYrQYo6vN9IXAo2D72HEwm0Bji2GDb64Esu4z+dYWq4y1jTL0xpg14CTgrgvNSSoXg9Vt/97kk0N994HbJSfs5pWd5SlefZia4XHj9fmcapCISSfBZB5TaWWiJWAkEK/vtsxK4yX59DfCGsVYsXQkstzPVpgKlwNpgZdrHrLbLwC7zL2HqWAXMFZFUOyhdDHwc+SVQSgXit4OKxxU4+PRs9/mc/ZDviX3SJ0h63ILXp4smx7Kww272/ModWB/ybuBxY8w2EbkHqDDGrAQeA54UkSqs3shy+9htIvIMVjDwArf3JAMEKtOu8k5ghYjcC2ywyyZEHU0i8gBWQDPAS8aYF0/pqiilens0bnfgv1Hd7hjp+dDTQ/tkW4Lb5Xi7VGiRzPlgjHkJazir77bv93ndAVwb5Nj7gPsiKdPevgcrG67/9lB1/B4r3VopNUx8PcEnyLBbb8/H4Q/5QD0ft0vwRrlH1t3dTU1NDR0dHVGtdyQkJydTWFhIQkLCiNURUfBRSsUfnwk97BYrcz4EeCaZxy10R7ldNTU1ZGRkUFxcfFIgHG2MMTQ0NFBTU8PUqVNHrB5dXkcpFZDPnjNxB53zsT4+nO75GE4ecgM74SDKPZ+Ojg5ycnJGdeABqweZk5Mz4j04DT5KqYB6ssWCBZ+ezU73fPzGDPjA97gFv/kkaSJaRnvg6RGN89Dgo5QKyG9C93xEBI9L8Dmc0mzMwJ6PJ1aGBFVQGnyUUgH1ZrsFCT497zn9Ae83IPTv+VgfbXqvz0BvvvkmV111FQCdnZ1ceumlzJ8/n6effjqq7dCEA6VUQL4Igo/HJb1zQ04xmAHrn/T0fLr1Xp+QNmzYQHd3Nxs3box63drzUUoF5AtzkynERs+HAMNuCT09H4dvgI22ffv2MWvWLG666Sbmzp3LNddcQ1tbG6+88gqzZs1i0aJFPPfccwDU1tZyww03sHHjRubPn8/u3buj2lbt+SilAupdXidUz8ftcjzbzW/MgGE3t8P3IP3or9v4+NCxYS2z7LRx/OCzc8LuV1lZyWOPPcbChQv56le/ygMPPMDDDz/MG2+8QUlJCddddx0A+fn5PProo/zHf/wHL7zwwrC2NRLa81FKBRRueR2IjZ5PoISDBHv1hWjf6xMLioqKWLhwIQA33HADFRUVTJ06ldLSUkSEG264weEWWrTno5QKKJKEg1jIdvObganBPfcgOTXsFkkPZaT0vxYtLS0xmQKuPR+lVEDhlteBGOn5YAY8b8Xjjt+EgwMHDvDBBx8A8NRTT3HppZeyd+/e3jmdp556ysnm9dLgo5QKqDfhwB2u5+P8sFv/+JgQx6nWs2fP5oknnmDu3Lk0NjbyzW9+k0ceeYTPfOYzLFq0iClTpjjdRECH3ZRSQfQElWDP84EY6fkEWuGg5ybTOOz5uFwufvOb35y0benSpezYsWPAvosXL2bx4sVRatnJtOejlArI25twEPxjwuNyxcB9PgMTDpIS3AB0euOv5zNaaPBRSgXU2/MJ8SkRCz2fQGu7JdrDbp1enxNNckxxcTFbt251uhkR0eCjlAqoZ74kIcjD5MCaD3I6280MXOCARI/V5q4o93xMgMc7jEbROA8NPkqpgLp94YNPLPR8DAPTi5McCD7Jyck0NDSM+gDU8zyf5OTkEa1HEw6UUgH1pCmHusk0NrLdzIBst57gE805n8LCQmpqaqirq4tanSOl50mmI0mDj1IqoJ5MsZjv+QRY4SDJYyUcRLPnk5CQMKJP/hxrdNhNKRVQz5xP6Pt8YnNtt0QHej5qcDT4KKUC6hl2SwiR7uZ2ieMrRwe6yfSThIP4ynYbTTT4KKUC6k048ATv+SS4XY4vYWPd5xMk4SDOHqkwmmjwUUoF1NOjCXWTaaJHeoOUU/wBsst6h926NfjEKg0+SqmAeofdQsz5WD0fhz/gzcAbYT0uwSXa84llGnyUUgF5/X7cLgm5HH8sDLsFSjgQERI9rqjfZKoip8FHKRWQ12dC3uMDVvBxuncRaG03sNKtNdstdmnwUUoF1OXz966RFkyiOxbmfAaucADWvI8Gn9gVUfARkaUiUikiVSJyV4D3k0Tkafv9NSJS3Oe9u+3tlSJyebgyRWSqXcYuu8zEUHWISLGItIvIRvvr5LXElVJD4vWZkPf4gD3s5vAHvDEDHyYH1uKi8baw6GgSNviIiBt4CLgCKAOuF5GyfrvdAjQZY0qAB4H77WPLgOXAHGAp8CsRcYcp837gQWNMKdBklx20DttuY8x8++u2QV0BpVRAXr8fT5ieT4LH+Tkfa223gduTEnTOJ5ZF0vM5F6gyxuwxxnQBK4Bl/fZZBjxhv34WuESsfvAyYIUxptMYsxeosssLWKZ9zBK7DOwyrw5Th1JqBHT7DAkRzvk4uZhmoIfJgdXz0eATuyIJPgVAdZ/va+xtAfcxxniBFiAnxLHBtucAzXYZ/esKVgfAVBHZICJviciFEZyTUioMry98zyfRHpZzcn23QI9UAOuBcjrnE7siWVg00P9r/5+0YPsE2x7oJzrU/qHqOAxMNsY0iMjZwPMiMscYc+ykBorcCtwKMHny5ABFKaX66vZHNucD1moIoRYgHUk+v8EdoIeWpD2fmBbJT0sNUNTn+0LgULB9RMQDZAKNIY4Ntr0eyLLL6F9XwDrsIb0GAGPMemA3MKP/SRhjHjHGlBtjyvPy8iI4baXiW7c3fLZbb/DxOtfz8ZuBy+uANeejCQexK5Lgsw4otbPQErESCFb222clcJP9+hrgDWMNAq8EltuZalOBUmBtsDLtY1bbZWCX+ZdQdYhInp3AgIhMs+vYE/klUEoF4o2k5xMDa6j5jQn4qO8kj5t2XV4nZoUddjPGeEXkDmAV4AYeN8ZsE5F7gApjzErgMeBJEanC6vEst4/dJiLPAB8DXuB2Y4wPIFCZdpV3AitE5F5gg102weoALgLuEREv4ANuM8Y0Dv2SKKXAGkoLta4bfDLn4+S9Pn5jcAfo+aQmumnv8gY4QsWCiB4mZ4x5CXip37bv93ndAVwb5Nj7gPsiKdPevgcrG67/9oB1GGP+BPwp7EkopQbF6zMh13WDk+d8nBLsJtPURDdtXTrsFqt0hQOlVEBef/ieT0wEH78JuLxOSqKbdg0+MUuDj1IqoC6f6Z3TCaYn+HQ5mnAQONstNdFNW7fP0XuQVHAafJRSAXl9/rA3mSZ6nJ/z8fkD32SamujB5zeOL3yqAtPgo5QKKNK13cDZ4GMMARMOUhLcADr0FqM0+CilAuqOZG03d+ymWqcmWsHnhAafmKTBRykVkDfCtd0ARxcX9RkT8CbTlMSeno+mW8ciDT5KqYC6vH4SwyQcJPaucOBsqnWg4JOaaN1JounWsUmDj1IqoC6fnySPO+Q+CTGQcBAs1bpn2E2DT2zS4KOUCqiz2xe25xMrcz6BUq0/GXbT4BOLNPgopQLq8g1i2M3JOZ+gqdba84llGnyUUgP4/YZunyEp4ptMYy/VOjWhZ85HEw5ikQYfpdQAPcNo4Xo+PcHJyUcXBE21TtKeTyzT4KOUGqDTfhRBuISDZPtGzg4HH13gC/IY7fQkq+fT2tEd7SapCGjwUUoN0OmzeguR9nw6up3rXQQbdktOcJPkcXGsQ4fdYpEGH6XUAD1zOOHmfFwuIdHtotPBOR9fkFRrgHEpCRxr155PLNLgo5QaoDPC4APW46qd7PlYcz6Bo8+4ZA/HdNgtJmnwUUoN0NPzSQyzthtYw1uOJhz4Ay+vAz09Hx12i0UafJRSA/T2fBIiCT4uRxMOrOV1Ar83LjlBEw5ilAYfpdQAn/R8Qme7ASR73I4PuwVa4QDsno8mHMQkDT5KqQF6htEi6/k4H3wCpVqDPeejCQcxSYOPUmqAwc35OD/sFijVGnp6Pt36KO0YpMFHKTVAb/CJINvN6YSDkKnWyQl0+4yjwVEFpsFHKTXAoFKtPW6Hez7BU60zUxIAaG7vimaTVAQ0+CilBhhMzycpwUWHgz0fE+RhcgDZaYkANBzX4BNrNPgopQboTTgIs7YbWNlunU6u7RZi2C0vwwo+9cc7o9giFQkNPkqpAToHNecTuysc5KQlAVCvPZ+Yo8FHKTXAYOZ8YiHVOtiwW26GFXwatOcTczT4KKUGaO/y4ZJIg4+LDgcXFg2Vap2W6CY5waXDbjEoouAjIktFpFJEqkTkrgDvJ4nI0/b7a0SkuM97d9vbK0Xk8nBlishUu4xddpmJ4eqw358sIsdF5H8P9iIopU7W3u0jJcEd9ObNvpI9bnx+g9fnTAAKNecjIuSkJWnCQQwKG3xExA08BFwBlAHXi0hZv91uAZqMMSXAg8D99rFlwHJgDrAU+JWIuMOUeT/woDGmFGiyyw5aRx8PAi9HeuJKqeDaunykJHoi2rfngXLtDgy99dw8GipI5mYkUac9n5gTSc/nXKDKGLPHGNMFrACW9dtnGfCE/fpZ4BKxfhqWASuMMZ3GmL1AlV1ewDLtY5bYZWCXeXWYOhCRq4E9wLbIT10pFUxHt4+UxMhG5XseV93uwOOqfX4r+HiCdX2A3LRETTiIQZH8dBUA1X2+r7G3BdzHGOMFWoCcEMcG254DNNtl9K8rYB0ikgbcCfwo1EmIyK0iUiEiFXV1dWFOWan41t5lDbtFIs3uIR3vjP4Cnl47+LjdwYPPhMxkjrS0R6tJKkKRBJ9A/6v9F0oKts9wbQ9Vx4+whumOB3j/kx2NecQYU26MKc/Lywu1q1Jxr6078mG3tCRrvxOd0e/5eCPo+RSOT6GprZu2Ll3dOpZE8tNVAxT1+b4QOBRknxoR8QCZQGOYYwNtrweyRMRj92767h+sjvOAa0Tkp0AW4BeRDmPMLyM4N6VUAB1dPlIiWNEaIM0ednOi5+Pz9QSf4G0tyEoB4GBTO6UTMqLSLhVeJD9d64BSOwstESuBYGW/fVYCN9mvrwHeMNZM4EpguZ2pNhUoBdYGK9M+ZrVdBnaZfwlVhzHmQmNMsTGmGPg58GMNPEqdmp5st0ik9/Z8nBh2szLsPCGG3QrHW8GnplmH3mJJ2J6PMcYrIncAqwA38LgxZpuI3ANUGGNWAo8BT4pIFVZvZLl97DYReQb4GPACtxtjfACByrSrvBNYISL3AhvssglWh1Jq+LV1eUlNTI1o31R7eO6EA8NaPQkHwR4mB1CQZZ3HwSYNPrEkokFdY8xLwEv9tn2/z+sO4Nogx94H3BdJmfb2PVjZcP23B62jzz4/DPW+UioyHd3+3hTqcNJjfM4nPyOJBLdwUHs+MUVXOFBKDdA+iFTrnjkfJ4bdPun5BG+ryyUUZKVwoKEtWs1SEdDgo5QawBp2iyzbLTUGUq1D9XwASvLTqaoNmRCrokyDj1LqJH6/GdSwm9slpCS4nUk4sJf0CTXnA1CSn8Ge+uOOLQGkBtLgo5Q6Sc+K1pFmu4F1r88JB1Y4iLTnU5qfTrfPsL9Rh95ihQYfpdRJem7GjPQ+H7DmfZyc8/G4Q7e1dEI6ALuO6tBbrNDgo5Q6SU/WWnpyQsTHpCV6HLrPJ7Kez/S8dERgx5Fj0WiWioAGH6XUSVo7u4FPUqgjkZ7kodWRnk9kcz5pSR5m5Gewsbo5Gs1SEdDgo5Q6yfEOK4hkJEcefMalJHCsvXukmhSU1xdZzwfgzMlZbDjQ3PsYBuUsDT5KqZO0DiH4ZKUm0OJA8IlkhYMe84uyaGnvZm/9iZFuloqABh+l1El67tcZzLDb+NQEmtsc6Pn0JhyEDz5nTRkPQMW+phFtk4qMBh+l1El65m7SB9XzSaS920dHlJ9m6u2d8wn/UVaan86EcUm8tVOf5xULNPgopU7SO+eTFHm2W2aKtW+0530GM+cjIlw8I493dtXpzaYxQIOPUuokxzu7cbuE5EHc55OVagWf5igHH98ght0AFs/M51iHl/X7dejNaRp8lFInOd7hJSPZg0hkH+gAWSmJAFGf94n0Pp8eF8/IIyXBzfMb+z8PU0WbBh+l1ElaO72DSjaAPj2ftq6RaFJQPT0fV4SBMi3Jw9LTJ/Li5kNRn59SJ9Pgo5Q6SWvH4INPz5xPtIfduuy5m4Qwy+v09YWzCjjW4eXFzYdHqlkqAhp8lFIn6Rl2G4yenk9LlIfduuxFUJM8kX+ULZyey8wJGTz89m78fr3h1CkafJRSJ2lq6yLTnsOJVHqSh0SPi/oTnSPUqsB6gk/iIIKPyyXctngaO48e55VtR0aqaSoMDT5KqZM0t3UzPjXyNGuw0pjz0pOoOxbl4OMbfPAB+Ozc05g1MYP7Xtyucz8O0eCjlDpJU1sX2WmD6/kA5GUkUdvqUM9nEHM+YD2C4Yefm8PB5nYeeHXnSDRNhaHBRynVq73LR6fXT1bq4INPfkYSta0dI9Cq4Lq8flwS/nk+gSyYlsOXz5vMI2/vYfWO2hFonQpFg49SqleTnSo92GE3gPxxDvR8fP5BD7n19b2rypg9aRz/46kNbD3YMowtU+Fo8FFK9Wo8YQWfofV8kmlu66bTG705lC6vf9BDbn0lJ7h5/OZyMlMSuPHxtWyp0QAULRp8lFK9elYoGFLPJyMJgLoo9n46vX4SPe5TKmNSZgp/+Np5pCS4+buHP+BvmgEXFRp8lFK9eobdhpJwkD/OCj5Ho5jx1uX1D+oen2CKc9P48+0XUDohnVufXM/3/7JVs+BGmAYfpVSvnuVxhjLsdlpWCgCHmtuHtU2hnOqcT1/5Gcn88bbzuWXRVH73wX4u//nbrK7URISRosFHKdWr/ngXIkMbdisanwrAgca24W5WUF1e3ynN+fSX5HHzvavK+MPXzsPtEv7+t+v4+u8q2F13fNjqUJaI/tdEZKmIVIpIlYjcFeD9JBF52n5/jYgU93nvbnt7pYhcHq5MEZlql7HLLjMxVB0icq6IbLS/NonI54d6MZSKd7WtHeSkJQ0pdTktyUNuehIHGqIZfIav59PXwpJcXvmfF3HXFbN4v6qeyx58m7uf28yRluimko9lYf/XRMQNPARcAZQB14tIWb/dbgGajDElwIPA/faxZcByYA6wFPiViLjDlHk/8KAxphRosssOWgewFSg3xsy363hYRAa3MJVSCrDmaybYczdDMTk7Jbo9n2Ecdusv0ePitoun89a3P8VXFkzh2fU1XPzvq/nJyzuivobdWBTJ/9q5QJUxZo8xpgtYASzrt88y4An79bPAJWI9DGQZsMIY02mM2QtU2eUFLNM+ZoldBnaZV4eqwxjTZozx2tuTAV0pUKkhqm3t6M1aG4rJ2alRDT7tXT5SEk4t2y2c3PQkfvi5Obz+rcVcecYkHn57Nxf+9A3+39t7eldYUIMXSfApAKr7fF9jbwu4jx0IWoCcEMcG254DNPcJJn3rClYHInKeiGwDtgC39Tm+l4jcKiIVIlJRV6fPcFcqEKvnkzzk4ydnp3K4pT1qH8ptXT5SEkc2+PSYnJPKg9fN58X/cSFnTh7PfS9tZ+nP39bVEYYokuAT6ClN/XsXwfYZru0h22GMWWOMmQOcA9wtIgN+e4wxjxhjyo0x5Xl5eQGKUiq+eX1+Go53nlLPZ3p+On4De+qjM0Hf3u0jNUrBp0fZaeN44qvn8tubzwHg7/9rHTf/di3VUezxjQWRBJ8aoKjP94VA/2fQ9u5jz7dkAo0hjg22vR7I6jNn07euYHX0MsZsB04Ap0dwXkqpPhpOdOE3kH8KPZ9ZE8cBsONw63A1K6QTnT5SE52Z4v3UrHxe+eeL+NfPzGbd3kaW/vxt/rBmP8boyH8kIgk+64BSOwstESuBYGW/fVYCN9mvrwHeMNb/wEpguZ2pNhUoBdYGK9M+ZrVdBnaZfwlVh12GB0BEpgAzgX0RXwGlFABHj1mZXKfS85mWl0ai28X2I8eGq1khtXd5o97z6SvR4+JrF05j1TcvYv7kLL77563c+Pha6o9Hd4270Shs8LHnT+4AVgHbgWeMMdtE5B4R+Zy922NAjohUAd8C7rKP3QY8A3wMvALcbozxBSvTLutO4Ft2WTl22UHrABYBm0RkI/Bn4BvGmPqhXQ6l4ld1o3VzaFF26pDLSHC7KMlPj0rPxxhDmwPDboEUjk/l97ecx79dfTpr9zZy1S/e5aMDTU43K6ZF1F81xrwEvNRv2/f7vO4Arg1y7H3AfZGUaW/fg5UN1397wDqMMU8CT4Y9CaVUSNVN1pzFqQQfgNmTxvHWzlqMMVgJrCOjo9uPMTg27NafiPCVBVM4a3IW//j7j7ju4Q+49+rTue6cyU43LSbpCgdKKcBamWB8agLpSaf2YX7WlCzqj3exb4RvNm3rspJaY6Hn09ec0zL56x2LOH96Lnf+aQv/+founQcKQIOPUgqA6sa2U+71AJxbnA3Aur2NYfY8NW1d1sKf0Uq1HozM1AQeu6mcL5xVwM9e3ckPV27TANSPBh+lFAA1Te3DEnxK8tPJTktk7b6RDT6tHVbPZ1zy4Nehi4YEt4ufXTuPr184lSc+2M+/vbBdA1AfsTFYqpRylNfnp6apjcvnTDzlskSEBdOyeXtnHX6/weUamXmf5vaeFbhjM/iAdS2+c+VsvH7D4+/tJTnBxbeXznK6WTFBez5KKfY3ttHtM5Tkpw9LeZeVTaS2tZMN1c3DUl4gPQ++i+XgA1YA+v5VZXzpvMn86s3d/P7D/U43KSZo8FFKsfOIlRo9Y8LwBJ9PzcrH4xJWjeBTQXuDT8rgnz0UbSLCvy07nSWz8vnBym28u0vvBtHgo5Ri59HjiDBsPZ/MlAQWz8zjuY8Ojtg6b6Nh2K0vt0v4v8vnU5KXzjf+sD6qj56IRRp8lFLsrG2laHzqsN4zc8OCKdQf7+TlrYeHrcy+Wtq6SfK4SB7hVa2HU0ZyAo/eVA7AHU99FNerYmvwUUpReaR12IbcelxUmse03DR+tXo3Pv/wZ3kdPdZB/ik8e8gpRdmp/PSauWyuaeGnr+xwujmO0eCjVJxr7ehmd91x5hZmDWu5Lpfwvy6bSeXRVv5YUR3+gEE61NLBpHEpw15uNCw9fRI3nj+FR9/dG7fzPxp8lIpzm2taMAbmFw1v8AG48oyJnFM8nvte3D7sjxw40tLBpKyhr8DttO9cOZtpuWnc9dzm3tUa4okGH6Xi3AZ7Acx5IxB8RIQH/m4+AP/w5Hpa2ofn8dM+v+FISwcTM0dv8ElOcHP/NXOpaWrn31dVOt2cqNPgo1Sc23CgmWl5aWSmjEzWWFF2Kr/88lnsqm3lK4+t4UhLxymXub/hBF0+P9PzhneeKtrOKc7mxvOn8F/v72PTCN4TFYs0+CgVx7p9ftbsbWTBtJwRrefiGXn8+stnU1V7nKv+8x3+vKHmlJaaqbTvS5o1MWO4muiYf7l8JjlpSfzor/G1/psGH6Xi2MbqZo53ermwJHfE67q0bALP376QgvGpfPPpTVz+87dZsfbAkIbiNlQ3k+h2MWPC6A8+GckJfPvymXx0oJmVm/o/JHrs0uCjVBx7Z2cdLoELpo988AGYMSGDP//jBTx43TxcItz13BbK732Vmx5fy+Pv7mXHkWP4w6RlG2N4s7KW8uLxo+oen1CuObuQ0wvG8ZOXd8RN8oEuLKpUHHttey3zi7LIjOIqAS6X8PkzC7l6fgEbq5t5ZesR/vbxUe554WMActISWTA9hwVTszlnajYz8jNOWpz0gz0N7Dx6nBvPL45am0eayyV8/6o5/N3DH/C7D/Zz28XTnW7SiNPgo1Sc2lN3nI8PH+NfPzPbkfpFhDMnj+fMyeO5+8rZHGxu54PdDbxfVc/7uxt4cbO1MkJGsofyKeMpL84mwS38v3f2Ujg+hS+cVeBIu0fKuVOzuXhGHg+/tZsbFkw55Yf6xbqxfXZKqaBesD/cPzN3ksMtsRRkpXDN2YVcc3YhxhhqmtpZt6+RdfuaqNjXyOpKKx151sQMfnH9mTHz+Ozh9M1Pz+Dqh97jiff3cfunSpxuzogae/97Sqmw/H7Dcx/VcG5xNpMyY2+VABGhKDuVouxUvnBWIQAt7d34/IbxqQmIjMwzgpw2vyiLS7PH1U8AABjoSURBVGbl88jbe7jx/ClkxOiD8oaDJhwoFYfe3lXHvoY2vrxgstNNiVhmSgLZaYljNvD0+OdLZ9DS3s1Taw843ZQRpcFHqTj0xPv7yMtI4orTY2PITX3ijMJMzp+Ww2/f20e3b+yueq3BR6k4s7mmmdWVddy4YAqJHv0IiEVfv2gqh1s6eGnLyDyOIhboT55SceZnf9vJ+NQEbl5Y7HRTVBCLZ+QzPS+NR97eM2ZXPdDgo1QceWtnHW/trOO2i6eP6cns0c7lEr524TS2HTrGmr2NTjdnRGjwUSpOtHf5+NfntzA9L017PaPA1fMLyEj2jNnEAw0+SsWJn7y8nerGdn78+TNI8oyNZWnGspREN184s4CXtxyh6USX080ZdhEFHxFZKiKVIlIlIncFeD9JRJ62318jIsV93rvb3l4pIpeHK1NEptpl7LLLTAxVh4h8WkTWi8gW+98lQ70YSo1VL285zBMf7OerC6dy3givYK2Gz/XnTabL5+dPH9U43ZRhFzb4iIgbeAi4AigDrheRsn673QI0GWNKgAeB++1jy4DlwBxgKfArEXGHKfN+4EFjTCnQZJcdtA6gHvisMeYM4CbgycFdAqXGth1HjvHtZzczryiLu66Y5XRz1CDMmjiOsyZn8d9rD4y5xINIej7nAlXGmD3GmC5gBbCs3z7LgCfs188Cl4h1J9gyYIUxptMYsxeosssLWKZ9zBK7DOwyrw5VhzFmgzGmZx3ybUCyiCRFegGUGssONrdz0+NrSU1y8+svn6Wp1aPQl86bwp66E1Tsb3K6KcMqkp/EAqC6z/c19raA+xhjvEALkBPi2GDbc4Bmu4z+dQWro68vAhuMMZ39T0JEbhWRChGpqKurC3PKSo1+h5rb+cqja2jr8vHEV8/ltKzYW0ZHhXfF6RNJSXDz5w0HnW7KsIok+ARay6J//y/YPsO1PWw7RGQO1lDcPwTYD2PMI8aYcmNMeV5eXqBdlBoz9jec4NrffEBdaye/vfkcZk0c53ST1BClJXm4bM4EXtpymC7v2FnxIJLgUwMU9fm+EOj/uL3efUTEA2QCjSGODba9Hsiyy+hfV7A6EJFC4M/AjcaY3RGck1Jj1vr9jXzx1x/Q1uXlqVsXUF6c7XST1Cm6en4BzW3dvLVz7IzaRBJ81gGldhZaIlYCwcp++6zEmuwHuAZ4w1izYyuB5Xam2lSgFFgbrEz7mNV2Gdhl/iVUHSKSBbwI3G2MeW8wJ6/UWLNi7QGWP/Ih6Ulu/njbBZxekOl0k9QwWFSaS3ZaIs+PoaG3sI9UMMZ4ReQOYBXgBh43xmwTkXuACmPMSuAx4EkRqcLqjSy3j90mIs8AHwNe4HZjjA8gUJl2lXcCK0TkXmCDXTbB6gDuAEqA74nI9+xtlxljaod2SZQafU50ernnrx/zdEU1F5bm8svrz4rq00nVyEpwu/js3EmsWFfNsY5uxo2B1SlkrKXvRaK8vNxUVFQ43QylhsXG6mb+ecUG9je28Y3F0/nmpTPwuDWrbazpGU79+XXzufpMZ57iKiLrjTHlw1GW/oQqNUp1dPt44G+VfPHX79PtM6z4+gL+5fJZGnjGqDOLxpOfkcSqbUecbsqw0CeZKjUKvV9Vz3ef38re+hN8/swCfvi5OWSmjP6hGBWcyyVcNmcCf1p/kI5uH8kJo3uJJP0TSalRpP54J996ZiNfenQNfmN48pZzefC6+Rp44sTlcybS3u3j7TGQ9aY9H6VGgY5uH799bx8Pra6i0+vjjk+VcMeSklH/168anAXTchiX7GHVtqNcNmei0805JRp8lIphxhhe2HyYn7y8g4PN7Vw6ewJ3XzmL6XnpTjdNOSDB7eLS2RN4bftRun1+Ekbx/J4GH6Vi1Lp9jfz4pe1sONDM7Enj+Pdr5nJBSa7TzVIOu2zORJ7bcJB1extH9c+DBh+lYszG6mZ+9rdK3tlVT35GEj/94ly+eHYhblegFaZUvFlUmkuCW3hzZ50GH6XUqdt2qIUHX93Ja9tryU5L5LtXzuaGBVNISdR5HfWJ9CQP5xRn82ZlLd+5crbTzRkyDT5KOWzn0VZ+/tpOXtpyhHHJHv73ZTO4eeFU0pP011MFtnhmHj9+aQeHmttH7Wrl+tOtlEM21zTz0OoqVm07Slqim39aUsItF07TtGkV1uKZ+fz4pR28tbOO68+d7HRzhkSDj1JRZIxhzd5GHlpdxTu76hmX7OGflpRw88KpZKclOt08NUqU5qdzWmYyb1bWavBRSgVnjGF1ZS0Prd7N+v1N5KYncdcVs/jyeZPJGAOLRKroEhEunpnPXzcdosvrH5VPqNXgo9QI8vr8vLT1CL9+czfbDx+jICuFf1s2h2vLi/QGUXVKFs/M46m1B/joQBMLpvV/qHPs0+Cj1Ag41tHN02ur+a/393GwuZ3peWn87Np5fG7+aaP6xkAVOxZMy8El8P7uBg0+SsW7mqY2fvvePp5eV83xTi8LpmXzo8/NYcmsfFx6n44aRpkpCZxekMmHuxvg0063ZvA0+Cg1DDYcaOLRd/fy8pbDuES4au4kblk0jTMK9UmiauScPz2Hx9/dS3uXb9TdD6bBR6kh8vr8vLb9KI++s5eK/U1kJHv4+kXTuPmCYiZljs57L9Tocv60HB5+aw8V+xu5sDTP6eYMigYfpQap8UQXK9Yd4A8fHuBgcztF2Sn84LNlXFtepDeGqqg6pzgbj0t4f3eDBh+lxqrNNc088f5+/rrZSm+9YHoO37uqjE+XTdB115Qj0pI8zCvK4oPdDU43ZdA0+CgVQqfXx0tbDvPE+/vZWN1MaqKb68qLuPH8KZROyHC6eUpxwfQcfvXmblo7ukfVPWMafJQK4HBLO3/48ABPrT1Aw4kupuWm8cPPlvGFswsZN4p+wdXYd/60HP7zjSrW7WtkyawJTjcnYhp8lLL5/YZ3qur57zX7eW17LX5juGTWBG66YAoLp+dqqrSKSWdOHo/HJazb16TBR6nR5OixDv5YUc1Ta6s52NxOdloiX1s0lRsWTKEoO9Xp5ikVUkqimzkFmazf1+R0UwZFg4+KSz6/4e1ddTy15gCv76jF5zcsLMnh7itn8emyCSR5Rtc9Eyq+lU8Zz+8/3D+q1nnT4KPiytFjHTyzrpoV66xeTk5aIl+7cCrXnzOZ4tw0p5un1JCUTxnPY+/uZeuhFs6aPN7p5kREg48a87w+v9XLWVvNG3YvZ1FJLt+5cjafLpswav5SVCqYs4utgLN+X5MGH6WcVlXbyh/X1/Dnjw5S29pJbnoiX79wGsvPKdJejhpT8jOSmZKTSsX+Rr7ONKebE5GIgo+ILAX+L+AGHjXG/KTf+0nA74CzgQbgOmPMPvu9u4FbAB/wT8aYVaHKFJGpwAogG/gI+IoxpitYHSKSAzwLnAP8lzHmjiFeCzUGHOvo5q+bDvHHiho2VjfjdgmfmpnPteWFfGpmvvZy1Jh19pTxvFVZhzEGkdjPzAwbfETEDTyEtW5qDbBORFYaYz7us9stQJMxpkRElgP3A9eJSBmwHJgDnAa8JiIz7GOClXk/8KAxZoWI/MYu+9fB6gA6gO8Bp9tfKs74/Yb3dtfz7PoaXtl6hE6vnxkT0vnulbO5+swC8jKSnG6iUiOufEo2z310kH0NbUwdBT37SHo+5wJVxpg9ACKyAlgG9A0+y4Af2q+fBX4pVuhdBqwwxnQCe0Wkyi6PQGWKyHZgCfAle58n7HJ/HawOY8wJ4F0RKRnEeasxYH/DCZ5dX8Of1tdwqKWDccke/q68iGvLCzmjIHNU/PWn1HApt+d9KvY1jpngUwBU9/m+Bjgv2D7GGK+ItAA59vYP+x1bYL8OVGYO0GyM8QbYP1gd9RGcgxojmk508eKWwzy/4SAV+5sQgQtL87jbTh7Qp4OqeFWSl864ZA8fHWji2vIip5sTViTBJ9CfjybCfYJtDzTwHmr/SNsRlIjcCtwKMHny5EgPUzGgo9vH69treX7jQd6srKXbZyjNT+dfLp/JF84q0McXKAW4XMK8oiw2Vbc43ZSIRBJ8aoC+YbQQOBRknxoR8QCZQGOYYwNtrweyRMRj93767h+sjogYYx4BHgEoLy+POGgpZ/j8hjV7Gnh+40Fe3nKE1k4vE8YlcfMFxVx9ZgFlk8bpsJpS/cwrzOLXb+0eFQ+XiyT4rANK7Sy0g1gJBF/qt89K4CbgA+Aa4A1jjBGRlcB/i8gDWAkHpcBarF7MgDLtY1bbZaywy/xLqDqGdtoqVm0/fIznNxzkLxsPceRYB+lJHpaePpHPn1nAgmk5+ugCpUKYV5SFz2/4+HALZ0/Jdro5IYUNPvb8yh3AKqy06MeNMdtE5B6gwhizEngMeNJOKGjECibY+z2DlZzgBW43xvgAApVpV3knsEJE7gU22GUTrA67rH3AOCBRRK4GLuuXjadi2IGGNl7YcoiVGw+x40grHpeweGYe/3rVbC6drfM4SkVqnv3Y9o3VsR98JB47D+Xl5aaiosLpZsS1g83tvLj5EC9sPszmGmuM+uwp47n6zAI+c8YkstMSHW6hUqPTBf/ndcqLs/nF9WcOe9kist4YUz4cZekKBypqjrR08OKWw7yw+RAbDjQDMLcwk+9cOYsrz5hE4XhdQVqpUzWvKItNNc1ONyMsDT5qRNW2dvDK1iO8sOkw6/Y3YgyUTRrHt5fO5KozTmNyjgYcpYbTvKIsXt56hKYTXYyP4REEDT5q2NUf72TVNivgrNnbgN/AzAkZfPPSGVw1dxLT8tKdbqJSY9a8wiwANtU0s3hmvsOtCU6DjxoWh5rbWbXtCK9sPcK6fY34DUzLS+OOJaVcNXcSMyZkON1EpeLCGYWZiMCm6hYNPmps2lt/gle2HuGVrYfZZCcNzJyQwR1LSlk6ZyKzJ2XovThKRVl6kofS/PSYn/fR4KMiZoxh++FWXtl2hFVbj1B5tBWw0jvvXDqLy+dM0CE1pWLAvMIs3thRG9MrXGvwUSH5/YaNNc12D+cIBxrbcAmcU5zNDz5bxmVzJlKQpcvbKBVL5hZl8cf1NdQ0tVOUHZtJPRp81AAd3T7eq6rnte1HeW17LXWtnSS4hQum5/KNxdO5tGwCuen6mAKlYtV8O+lgc02LBh8V2+paO3ljx1Fe/biWd6vq6Oj2k57k4eIZeVxals+SWRPITElwuplKqQjMnJhBotvF5ppmPjN3ktPNCUiDT5wyxrDz6HFe236UVz8+yqaaZoyBgqwUrisv4tKyCZw3NUef/KnUKJTocTH7tHFsrI7dpAMNPnGk2+dn7d5GXv34KK/vOEp1YztgJQx869IZXDJ7gmaoKTVGzC/M5Nn1Nfj8JiYX5NXgM8YdPdbBW5V1vLmzlnd21tPa6SXJ42JRSS7fWFzCkln5TBiX7HQzlVLDbG5hFk98sJ/ddcdj8j47DT5jjNfn56MDzbxZWcublXV8fPgYABPGJXHlGZO4ZHY+i0pzSU3U/3qlxrJ5RfZKB9XNGnzUyKg91sGbO+t4q7KOd3bVcazDi9slnD1lPHcuncXimXnMmqjDaUrFk2m5aWQkedhc0xKTj9XW4DMKeX1+NlY3s9ru3Ww7ZPVu8jOSuOL0SSyemcfC0lzGJWt2mlLxyuUSTi/IjNmVDjT4jALGGPbUn+DdXfW8W1XPh7sbaO20ezeTx/PtpTNZPCNfkwWUUieZV5TFY+/uodPrI8kTWw9l1OATo+qPd/JeVT3v7qrnvap6DrV0AFCUncJV805jUUkui0pz9d4bpVRQ8woz6fZZy2LNt+eAYoUGnxjR3uVjzd4G3quq551d9ew4Yq2blpmSwMKSHG4vyeXCkjx9/o1SKmI9SQeba5o1+CiLz2/YcrDFDjZ1fLS/mS6fn0S3i/Li8fzL5TO5sDSXOadlxmSOvlIq9k3KTCY3PYmN1c3ceL7TrTmZBp8oMcawv6GNd6usYbT3dzfQ0t4NwOxJ47h5YTGLSnI5pziblMTYGptVSo1OIsK8wkw22488iSUafEZQw/FO3t/d0JsocLDZWlHgtMxkLiubwKLSXBaW5OoinUqpETOvKIs3Kmtp7egmI4YyYDX4DCNjDJVHW3ntY2s16J51lTKSPVwwPYfbFk9nUUkuxTmpmpWmlIqKuYWZGANbDrZwwfRcp5vTS4PPMKht7eBP6w/yx4pq9tSfAKy/Nr716RlcNCOPMwp03kYp5Yy5fR6voMFnjGg43skvXt/FU2ur6fL5Oad4PLdcOJVPz55Avq6XppSKAdlpiUzOTmVTjK1wrcFniN7eWce3ntlIc1s315YX8vULp+kjpJVSMWluYSYbDmjwGfVW76jl1icrmJ6Xzh++toCZE2Nv0T6llOoxvyiLFzYfpq61k7yM2Ehw0ieFDVJdayfffGYjpfkZPP0P52vgUUrFvE/mfWKn96PBZ5AefXcPx9q7+cX1Z+rSNkqpUeH0gnF4XMLavY1ON6VXRMFHRJaKSKWIVInIXQHeTxKRp+3314hIcZ/37ra3V4rI5eHKFJGpdhm77DITh1rHcDPG8MKmwyyZNYGSfJ3fUUqNDqmJHs6dms3qylqnm9IrbPARETfwEHAFUAZcLyJl/Xa7BWgyxpQADwL328eWAcuBOcBS4Fci4g5T5v3Ag8aYUqDJLnvQdQz2QkSi7ngnB5vbWViSMxLFK6XUiFkyK5+dR49zoKHN6aYAkfV8zgWqjDF7jDFdwApgWb99lgFP2K+fBS4R6y7KZcAKY0ynMWYvUGWXF7BM+5gldhnYZV49xDqGXV1rJ2Ctl6SUUqPJFWdMIsEt/Pqt3U43BYgs260AqO7zfQ1wXrB9jDFeEWkBcuztH/Y7tsB+HajMHKDZGOMNsP9Q6uglIrcCt9rfHheRBqA+6FmHcMX9QzkqpuUyxGsxBum1sOh1+MSYuhY/sb+GIBeYMlztiCT4BLo130S4T7DtgXpcofYfSh0nbzDmEeCRnu9FpMIYUx7g2Lij1+ITei0seh0+odfCYl+H4uEqL5Jhtxqg7wPAC4FDwfYREQ+QCTSGODbY9nogyy6jf12DrUMppVSMiiT4rANK7Sy0RKzJ/ZX99lkJ3GS/vgZ4wxhj7O3L7Uy1qUApsDZYmfYxq+0ysMv8yxDrUEopFaPCDrvZ8yt3AKsAN/C4MWabiNwDVBhjVgKPAU+KSBVWb2S5few2EXkG+BjwArcbY3wAgcq0q7wTWCEi9wIb7LIZSh1hPBJ+l7ih1+ITei0seh0+odfCMqzXQazOg1JKKRU9usKBUkqpqNPgo5RSKuriMviEWy5oLBCRx0WkVkS29tmWLSKv2ksXvSoi4+3tIiK/sK/HZhE5q88xN9n77xKRmwLVFctEpEhEVovIdhHZJiL/094eV9dCRJJFZK2IbLKvw4/s7TG7nNVIs1db2SAiL9jfx+W1EJF9IrJFRDaKSIW9beR/P4wxcfWFleCwG5gGJAKbgDKn2zUC53kRcBawtc+2nwJ32a/vAu63X18JvIx1z9QCYI29PRvYY/873n493ulzG+R1mAScZb/OAHZiLekUV9fCPp90+3UCsMY+v2eA5fb23wD/aL/+BvAb+/Vy4Gn7dZn9O5METLV/l9xOn98Qr8m3gP8GXrC/j8trAewDcvttG/Hfj3js+USyXNCoZ4x5GysrsK++SxT1X7rod8byIda9VpOAy4FXjTGNxpgm4FWs9fNGDWPMYWPMR/brVmA71goYcXUt7PM5bn+bYH8ZYng5q5EkIoXAZ4BH7e9jemkvB4z470c8Bp9AywUNWI5njJpgjDkM1ocykG9vD3ZNxtS1sodLzsT6qz/uroU9zLQRqMX6cNhNhMtZAX2XsxrV18H2c+DbgN/+PuKlvRh718IAfxOR9WItQwZR+P2IxyeZRrQcT5w5paWLRgMRSQf+BPyzMeaY9Ydr4F0DbBsT18JY97/NF5Es4M/A7EC72f+O2esgIlcBtcaY9SKyuGdzgF3H/LWwLTTGHBKRfOBVEdkRYt9huxbx2POJ5+V4jtpdZOx/ex7uMdhlkEYVEUnACjx/MMY8Z2+Oy2sBYIxpBt7EGrOPx+WsFgKfE5F9WMPuS7B6QvF4LTDGHLL/rcX6o+RcovD7EY/BJ5LlgsaqvksU9V+66EY7k2UB0GJ3tVcBl4nIeDvb5TJ726hhj80/Bmw3xjzQ5624uhYikmf3eBCRFOBSrPmvuFvOyhhztzGm0FiLZC7HOrcvE4fXQkTSRCSj5zXWz/VWovH74XSmhRNfWBkbO7HGvL/rdHtG6ByfAg4D3Vh/ldyCNU79OrDL/jfb3lewHu63G9gClPcp56tYE6lVwN87fV5DuA6LsLr/m4GN9teV8XYtgLlYy1Vttj9cvm9vn4b1gVkF/BFIsrcn299X2e9P61PWd+3rUwlc4fS5neJ1Wcwn2W5xdy3sc95kf23r+TyMxu+HLq+jlFIq6uJx2E0ppZTDNPgopZSKOg0+Simlok6Dj1JKqajT4KOUUirqNPgopZSKOg0+Simlou7/A2bqLdwLyQsTAAAAAElFTkSuQmCC\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, f0_y, label = '0')\n",
"# plt.plot(test_q, fT_y, label = 'T')\n",
"# plt.plot(test_q, fplus_y, label = '+')\n",
"# plt.plot(test_q, res_y, label = 'res')\n",
"plt.legend()\n",
"plt.ylim(0.0, 6e-6)\n",
"# plt.yscale('log')\n",
"# plt.xlim(3080, 3110)\n",
"plt.savefig('test.png')\n",
"# print(jpsi_width)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"dtype = tf.float64\n",
"mixture = tfd.MixtureSameFamily(mixture_distribution=tfd.Categorical(probs=[tf.constant(0.007, dtype=dtype),\n",
" tf.constant(0.965, dtype=dtype),\n",
" tf.constant(0.04, dtype=dtype),\n",
" tf.constant(0.03, dtype=dtype),\n",
" tf.constant(0.006, dtype=dtype)]),\n",
" components_distribution=tfd.Uniform(low=[tf.constant(x_min, dtype=dtype), \n",
" tf.constant(3092, dtype=dtype),\n",
" tf.constant(3682, dtype=dtype), \n",
" tf.constant(3070, dtype=dtype),\n",
" tf.constant(3660, dtype=dtype)], \n",
" high=[tf.constant(x_max, dtype=dtype),\n",
" tf.constant(3100, dtype=dtype), \n",
" tf.constant(3690, dtype=dtype),\n",
" tf.constant(3110, dtype=dtype), \n",
" tf.constant(3710, dtype=dtype)]))\n",
"\n",
"# probs = mixture.prob(test_q)\n",
"# probs_np = zfit.run(probs)\n",
"# probs_np *= np.max(calcs_test) / np.max(probs_np)\n",
"# plt.figure()\n",
"# plt.semilogy(test_q, probs_np,label=\"importance sampling\")\n",
"# plt.semilogy(test_q, calcs_test, label = 'pdf')\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Adjust scaling of different parts"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"# total_f.update_integration_options(draws_per_dim=20000000, mc_sampler=None)\n",
"# inte = total_f.integrate(limits = (3080, 3112), norm_range=False)\n",
"# inte_fl = zfit.run(inte)\n",
"# print(inte_fl)\n",
"# print(pdg[\"jpsi_BR\"]/pdg[\"NR_BR\"], inte_fl*pdg[\"psi2s_auc\"]/pdg[\"NR_auc\"])"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"# print(\"jpsi:\", inte_fl)\n",
"# print(\"Increase am by factor:\", np.sqrt(pdg[\"jpsi_BR\"]/pdg[\"NR_BR\"]*pdg[\"NR_auc\"]/inte_fl))\n",
"# print(\"New amp:\", pdg[\"jpsi\"][3]*np.sqrt(pdg[\"jpsi_BR\"]/pdg[\"NR_BR\"]*pdg[\"NR_auc\"]/inte_fl))\n",
"\n",
"# print(\"psi2s:\", inte_fl)\n",
"# print(\"Increase am by factor:\", np.sqrt(pdg[\"psi2s_BR\"]/pdg[\"NR_BR\"]*pdg[\"NR_auc\"]/inte_fl))\n",
"# print(\"New amp:\", pdg[\"psi2s\"][3]*np.sqrt(pdg[\"psi2s_BR\"]/pdg[\"NR_BR\"]*pdg[\"NR_auc\"]/inte_fl))\n",
"\n",
"\n",
"\n",
"# 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": 14,
"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": 15,
"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": 16,
"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": 17,
"metadata": {},
"outputs": [],
"source": [
"# print(36000*(1+ pdg[\"jpsi_BR\"]/pdg[\"NR_BR\"] + pdg[\"psi2s_BR\"]/pdg[\"NR_BR\"]))"
]
},
{
"cell_type": "code",
"execution_count": 18,
"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": 19,
"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": 20,
"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": 21,
"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')\n",
"# 1-(0.21+0.62)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Toys"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"\n",
" \n",
"# print(list_of_borders[:9])\n",
"# print(list_of_borders[-9:])\n",
"\n",
"\n",
"class UniformSampleAndWeights(zfit.util.execution.SessionHolderMixin):\n",
" def __call__(self, limits, dtype, n_to_produce):\n",
" # n_to_produce = tf.cast(n_to_produce, dtype=tf.int32)\n",
" low, high = limits.limit1d\n",
" low = tf.cast(low, dtype=dtype)\n",
" high = tf.cast(high, dtype=dtype)\n",
"# uniform = tfd.Uniform(low=low, high=high)\n",
"# uniformjpsi = tfd.Uniform(low=tf.constant(3080, dtype=dtype), high=tf.constant(3112, dtype=dtype))\n",
"# uniformpsi2s = tfd.Uniform(low=tf.constant(3670, dtype=dtype), high=tf.constant(3702, dtype=dtype))\n",
"\n",
" list_of_borders = []\n",
" _p = []\n",
" splits = 10\n",
"\n",
" _ = np.linspace(x_min, x_max, splits)\n",
"\n",
" for i in range(splits):\n",
" list_of_borders.append(tf.constant(_[i], dtype=dtype))\n",
" _p.append(tf.constant(1/splits, dtype=dtype))\n",
" \n",
"# mixture = tfd.MixtureSameFamily(mixture_distribution=tfd.Categorical(probs=_p[:(splits-1)]),\n",
"# components_distribution=tfd.Uniform(low=list_of_borders[:(splits-1)], \n",
"# high=list_of_borders[-(splits-1):]))\n",
" mixture = tfd.MixtureSameFamily(mixture_distribution=tfd.Categorical(probs=[tf.constant(0.007, dtype=dtype),\n",
" tf.constant(0.95, dtype=dtype),\n",
" tf.constant(0.07, dtype=dtype),\n",
" tf.constant(0.025, dtype=dtype),\n",
" tf.constant(0.006, dtype=dtype)]),\n",
" components_distribution=tfd.Uniform(low=[tf.constant(x_min, dtype=dtype), \n",
" tf.constant(3093, dtype=dtype),\n",
" tf.constant(3683, dtype=dtype), \n",
" tf.constant(3070, dtype=dtype),\n",
" tf.constant(3660, dtype=dtype)], \n",
" high=[tf.constant(x_max, dtype=dtype),\n",
" tf.constant(3099, dtype=dtype), \n",
" tf.constant(3690, dtype=dtype),\n",
" tf.constant(3110, dtype=dtype), \n",
" tf.constant(3710, dtype=dtype)]))\n",
"# mixture = tfd.Uniform(tf.constant(x_min, dtype=dtype), tf.constant(x_max, dtype=dtype))\n",
"# sample = tf.random.uniform((n_to_produce, 1), dtype=dtype)\n",
" sample = mixture.sample((n_to_produce, 1))\n",
"# sample = tf.random.uniform((n_to_produce, 1), dtype=dtype)\n",
" weights = mixture.prob(sample)[:,0]\n",
"# weights = tf.broadcast_to(tf.constant(1., dtype=dtype), shape=(n_to_produce,))\n",
" # sample = tf.expand_dims(sample, axis=-1)\n",
"# print(sample, weights)\n",
" \n",
"# weights = tf.ones(shape=(n_to_produce,), dtype=dtype)\n",
" weights_max = None\n",
" thresholds = tf.random_uniform(shape=(n_to_produce,), dtype=dtype)\n",
" return sample, thresholds, weights, weights_max, n_to_produce"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"total_f._sample_and_weights = UniformSampleAndWeights"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [],
"source": [
"# psi2s_mass"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [],
"source": [
"# zfit.settings.set_verbosity(10)"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"6/6 of Toy 2/2\n",
"Time taken: 1 min, 27 s\n",
"Projected time left: \n"
]
}
],
"source": [
"# zfit.run.numeric_checks = False \n",
"\n",
"nr_of_toys = 2\n",
"nevents = int(pdg[\"number_of_decays\"])\n",
"nevents = pdg[\"number_of_decays\"]\n",
"event_stack = 1000000\n",
"# zfit.settings.set_verbosity(10)\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",
" \n",
" print(\"{0}/{1} of Toy {2}/{3}\".format(c, calls, toy+1, nr_of_toys))\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*(toy))*((nr_of_toys-toy)*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": 27,
"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": 28,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Time to generate full toy: 87 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": 29,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(5404696,)\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAD8CAYAAACVZ8iyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAW00lEQVR4nO3dfaxcdZ3H8fdnodIFulLKhTQt0EIKtopc6hWILAQF5cmUB3UpMVqUWNmFRNaH3bKalY0hQQm6Ie5iykIoCfIgD6VaXG1YEI3ycCttKZRKi1UuNO21LFBSYHn47h/zu3S4zO2de88583DO55VM5sxvzpn5nd+98/uc8ztnzigiMDOzavurdlfAzMzaz2FgZmYOAzMzcxiYmRkOAzMzw2FgZmY0EQaSDpR0n6R1kh6X9JVUvq+kFZKeSveTU7kkXS1pg6Q1kuYWvRJmZpZNM3sGbwBfi4jZwLHARZLmAIuAeyNiFnBvegxwGjAr3RYC1+ReazMzy9WoYRARmyPi92l6O7AOmAacCSxJsy0BzkrTZwI3Rs2DwD6SpuZeczMzy83uY5lZ0gzgKOAh4ICI2Ay1wJC0f5ptGvBM3WIDqWzzsNdaSG3Pgb322utD73vf+8ZRfbNqeezZF0ed54hp721BTXZtpHp2Qt3KZOXKlX+JiJ48XqvpMJC0N3AHcElEvCRpxFkblL3rmhcRsRhYDNDX1xf9/f3NVsWssmYsWj7qPP1XnNGCmuzaSPXshLqViaQ/5fVaTZ1NJGkCtSC4KSLuTMVbhoZ/0v3WVD4AHFi3+HTguXyqa2ZmRWjmbCIB1wHrIuL7dU8tAxak6QXA3XXln09nFR0LvDg0nGRmZp2pmWGi44DPAY9JWpXK/gW4ArhN0gXAn4HPpOfuAU4HNgA7gC/kWmMzM8vdqGEQEb+h8XEAgJMazB/ARRnrZWbWlNdff52BgQFeffXVdlelMBMnTmT69OlMmDChsPcY09lEZmadZmBggEmTJjFjxgx2cWJL14oItm3bxsDAADNnzizsfXw5CjPraq+++ipTpkwpZRAASGLKlCmF7/k4DMys65U1CIa0Yv0cBmZm5mMGZlYuzXwxbyw2jfGLcpdddhl77703X//61xs+v3TpUg477DDmzJmTR/Vy4z0DM7MWWrp0KU888US7q/EuDgMzs4wuv/xyDj/8cE4++WTWr18PwLXXXsuHP/xhjjzySD71qU+xY8cOfvvb37Js2TK+8Y1v0Nvby8aNGxvO1w4OAzOzDFauXMktt9zCo48+yp133skjjzwCwDnnnMMjjzzC6tWrmT17Ntdddx0f+chHmDdvHldeeSWrVq3i0EMPbThfO/iYgZlZBr/+9a85++yz2XPPPQGYN28eAGvXruVb3/oWL7zwAi+//DKnnHJKw+Wbna9oDgMzs4wanfp5/vnns3TpUo488khuuOEG7r///obLNjtf0TxMZGaWwQknnMBdd93FK6+8wvbt2/npT38KwPbt25k6dSqvv/46N91009vzT5o0ie3bt7/9eKT5Ws17BmZWKmM9FTSruXPncu6559Lb28vBBx/M8ccfD8B3vvMdjjnmGA4++GCOOOKItwNg/vz5fOlLX+Lqq6/m9ttvH3G+VlPtunLt5R+3MWtOM+fQt7ozbGSkehZRt3Xr1jF79uzcX7fTNFpPSSsjoi+P1/cwkZmZOQzMyibvb+BaNTgMzKzrdcJwd5FasX4OAzPrahMnTmTbtm2lDYSh3zOYOHFioe/js4nMrKtNnz6dgYEBBgcH212Vwgz90lmRRg0DSdcDnwS2RsQHUtmtwOFpln2AFyKiV9IMYB2wPj33YERcmHelzcyGTJgwodBfAKuKZvYMbgB+CNw4VBAR5w5NS7oKeLFu/o0R0ZtXBc2sxgeGrUijhkFEPJC2+N9Fte9g/x3wsXyrZWZmrZT1APLxwJaIeKqubKakRyX9StLxGV/fzMxaIOsB5POAm+sebwYOiohtkj4ELJX0/oh4afiCkhYCCwEOOuigjNUwM7Msxr1nIGl34Bzg1qGyiHgtIral6ZXARuCwRstHxOKI6IuIvp6envFWw8zMcpBlmOhk4MmIGBgqkNQjabc0fQgwC3g6WxXNzKxoo4aBpJuB3wGHSxqQdEF6aj7vHCICOAFYI2k1cDtwYUQ8n2eFzcwsf82cTXTeCOXnNyi7A7gje7XMzKyVfDkKMzNzGJiZmcPAzMxwGJiZGQ4DMzPDYWBmZjgMzMwMh4GZmeEwMDMzHAZmZobDwMzMcBiYmRkOAzMzw2FgZmY4DMyshWYsWt7uKtgIHAZmZuYwMDMzh4GZmeEwMDMzmggDSddL2ippbV3ZZZKelbQq3U6ve+5SSRskrZd0SlEVNzOz/DSzZ3ADcGqD8h9ERG+63QMgaQ4wH3h/WuY/Je2WV2XNzKwYo4ZBRDwAPN/k650J3BIRr0XEH4ENwNEZ6mdmZi2Q5ZjBxZLWpGGkyalsGvBM3TwDqexdJC2U1C+pf3BwMEM1zMwsq/GGwTXAoUAvsBm4KpWrwbzR6AUiYnFE9EVEX09PzzirYWZmeRhXGETEloh4MyLeAq5l51DQAHBg3azTgeeyVdHMzIo2rjCQNLXu4dnA0JlGy4D5kvaQNBOYBTycrYpmZla03UebQdLNwInAfpIGgG8DJ0rqpTYEtAn4MkBEPC7pNuAJ4A3gooh4s5iqm5lZXkYNg4g4r0HxdbuY/3Lg8iyVMrPym7FoOZuuOKPd1bDE30A2s5bylUs7k8PAzMwcBmZm5jAwMzMcBmZmhsPAzMxwGJiZGQ4DMzPDYWBmZjgMzErJX+yysXIYmJmZw8DMzBwGZtYGHsbqPA4DM2srB0NncBiYWcdxQLSew8DMzBwGZmbmMDAzM5oIA0nXS9oqaW1d2ZWSnpS0RtJdkvZJ5TMkvSJpVbr9qMjKm5lZPprZM7gBOHVY2QrgAxHxQeAPwKV1z22MiN50uzCfappZmfmAcfuNGgYR8QDw/LCyX0bEG+nhg8D0AupmZiXXzhBwAL1THscMvgj8vO7xTEmPSvqVpONHWkjSQkn9kvoHBwdzqIaZdQJ3st0pUxhI+ibwBnBTKtoMHBQRRwFfBX4s6W8aLRsRiyOiLyL6enp6slTDrNTcuVorjDsMJC0APgl8NiICICJei4htaXolsBE4LI+Kmlk15RGGDtTRjSsMJJ0K/DMwLyJ21JX3SNotTR8CzAKezqOiZmZWnGZOLb0Z+B1wuKQBSRcAPwQmASuGnUJ6ArBG0mrgduDCiHi+4QubWeG8RWzN2n20GSLivAbF140w7x3AHVkrZWaWxYxFy9l0xRntrkZX8TeQzUqum/YOuqmuZeMwMDMzh4GZmTkMzLrCeIZPyjzk0sy6lXn9i+AwMLOuN1LH70BonsPAzMwcBmbWGbwV314OAzOrHAfPuzkMzKxr1Hfi7tDzNeo3kM3M2iGvzt6h0RzvGZiZmcPArGq8pWyNOAzMKsAB0JjbZSeHgZlVhjv/kTkMzMzMYWBmY9NJW9c+4yg/DgMzM3MYmFXV8K1hXwm02poKA0nXS9oqaW1d2b6SVkh6Kt1PTuWSdLWkDZLWSJpbVOXNzBxQ+Wh2z+AG4NRhZYuAeyNiFnBvegxwGjAr3RYC12SvppmZFampMIiIB4DnhxWfCSxJ00uAs+rKb4yaB4F9JE3No7JmVZL3Fm8rtqDzeI/xvob3ELLJcszggIjYDJDu90/l04Bn6uYbSGXvIGmhpH5J/YODgxmqYWZV4k6/GEUcQFaDsnhXQcTiiOiLiL6enp4CqmFmzeiGzrXZOnbDunSqLGGwZWj4J91vTeUDwIF1800HnsvwPmZmTcszEKoULlnCYBmwIE0vAO6uK/98OqvoWODFoeEkMzPrTE39noGkm4ETgf0kDQDfBq4AbpN0AfBn4DNp9nuA04ENwA7gCznX2czMctZUGETEeSM8dVKDeQO4KEulzMystfwNZLMO1Oqx6tG+jVylsfOqchiYWSU40HbNYWBmb+vkDrNVdZuxaHlHt0NRHAZmFVbFTs8acxiYdTB31tYqDgMzs6TK4eswMOsQreyIqtzpjaTqZ1A5DMwsk27rNLutvq3iMDBrs07unDq5bpYvh4GZmTkMzMzMYWDWVh6GaT//DWocBmZm5jAws/Gp+qmYZeMwMGuhbu4wG9W9m9enWVVYR3AYmNkIxtIJVqXDLDOHgVmbuAO1TuIwMLNdcmhV47LW4w4DSYdLWlV3e0nSJZIuk/RsXfnpeVbYzKydyhoK4w6DiFgfEb0R0Qt8CNgB3JWe/sHQcxFxTx4VNetmZe1ArDzyGiY6CdgYEX/K6fXMKqlbQqNb6mnNyysM5gM31z2+WNIaSddLmtxoAUkLJfVL6h8cHMypGmbdzx1t5yvj3yhzGEh6DzAP+EkqugY4FOgFNgNXNVouIhZHRF9E9PX09GSthlnHGE9HUcbOxbpLHnsGpwG/j4gtABGxJSLejIi3gGuBo3N4D7OO1G1nmXRTXTtRmdsvjzA4j7ohIklT6547G1ibw3uYmVmBds+ysKQ9gY8DX64r/p6kXiCATcOeMzOzDpRpzyAidkTElIh4sa7scxFxRER8MCLmRcTm7NU063xlGUIoy3rY2PgbyGZm5jAwMzOHgVXYeIdDxrqch12sGzgMzAriELBu4jAwG0XenbpDohzK9nd0GJgVqGwdhpWXw8CsgVYdTzDrFA4DM7NxKlP4Owys9Ir6wNa/rs8wsm7/mzoMzMwy6vYgAIeBmZnhMLCKKcMWnFkRHAZmwzgwrIocBmbj4MCwsnEYmJmZw8C6z0hb5bvaWs9jS34sr+E9h2rrxr+/w8C62ngCYDxhkodu7CCsOhwG1lWa6VDd6Vq7dPP/XuYwkLRJ0mOSVknqT2X7Sloh6al0Pzl7Vc0ay/IBnLFoeeHDS2bdIK89g49GRG9E9KXHi4B7I2IWcG96bJabvId63Olb1RU1THQmsCRNLwHOKuh9rCLa8ZsCjeZxaFhZ5REGAfxS0kpJC1PZARGxGSDd7z98IUkLJfVL6h8cHMyhGtbNWtHJuiO3IpTl/2r3HF7juIh4TtL+wApJTzazUEQsBhYD9PX1RQ71MDOzccq8ZxARz6X7rcBdwNHAFklTAdL91qzvY9Uxni2tsmydmbVLpjCQtJekSUPTwCeAtcAyYEGabQFwd5b3sfIZz3i8O3yz4mTdMzgA+I2k1cDDwPKI+G/gCuDjkp4CPp4eW8W5MzfrXJnCICKejogj0+39EXF5Kt8WESdFxKx0/3w+1bVukvX8fzNrHX8D2QpVlp+D7NR6Wefp1v8Vh4HlrtkPQ7PHCFp5Ybpu/SCbZeUwsDFpx5e/zKx4DgMzM3MYWI230M3y1W2fKYeBmZk5DKyxvA4Ct+o1zCwbh4GZmTkMzMzMYVBa7Rp6yfLdADNrH4dBRY3WaTe7vJmVg8PAmtaOaw05dKzbdcv/sMOgwsYylFNf1i3/3NBddbXy6ab/P4dBRWQdFjKzcnMYlEy3bsGbWXs5DEpgeKefVwgUFSYOKbPO4zCokKI7YZ9Oata9HAYdoF2dtJnZkHGHgaQDJd0naZ2kxyV9JZVfJulZSavS7fT8qmv1Zixa3hE/Iu+wMet+WfYM3gC+FhGzgWOBiyTNSc/9ICJ60+2ezLUskSIv7NbOc/kdCGYj64bPx+7jXTAiNgOb0/R2SeuAaXlVzMzMWieXYwaSZgBHAQ+looslrZF0vaTJebxHFXXC1kQn1MHMipc5DCTtDdwBXBIRLwHXAIcCvdT2HK4aYbmFkvol9Q8ODmathrWBg8KsPDKFgaQJ1ILgpoi4EyAitkTEmxHxFnAtcHSjZSNicUT0RURfT09Plmp0jTw7z1Z2xO70zcovy9lEAq4D1kXE9+vKp9bNdjawdvzVs1Zzx29WTVn2DI4DPgd8bNhppN+T9JikNcBHgX/Mo6JV5c7ZzFohy9lEvwHU4CmfStoGDg0zy8LfQC5AEeftt+Kqow4Us+pyGOQgjx99aeY1Ov0CdGbWvRwGuzDWTrPZ+d0Zm1VPM5ePaSeHwTg1uyU/1ktHdPI/i5nloxM/5w6DUeRxvR93/GY2pFM/95UJg7yGfNpxldBO/ecxs/KoTBiYmdnIKhEGo/0s5GgHdpo54yfP00m9J2BmrTbuL52ZO20zK49K7BnkwR2/mZVZacMgyzn/7vjNrGpKOUw03rF9X+rBzKqqtHsGQ9rV+ZuZdZPShUErLgnhEDGzsilVGDRzCqmZmb1bqcLAzMzGx2FgZmYOAzMzcxiYmRkFhoGkUyWtl7RB0qKi3gd8YNjMuk+n9VuFhIGk3YD/AE4D5gDnSZpTxHsN6bSGNTPrJkXtGRwNbIiIpyPi/4BbgDMLei8zM8uoqMtRTAOeqXs8ABxTP4OkhcDC9PBlSduAvxRUn26zH26LIW6LGrfDTqVpC3030+L7AQfnU5PiwkANyuIdDyIWA4vfXkDqj4i+gurTVdwWO7ktatwOO7ktalI7zMjr9YoaJhoADqx7PB14rqD3MjOzjIoKg0eAWZJmSnoPMB9YVtB7mZlZRoUME0XEG5IuBn4B7AZcHxGPj7LY4lGerxK3xU5uixq3w05ui5pc20ERMfpcZmZWav4GspmZOQzMzKxDwqCVl65oF0nXS9oqaW1d2b6SVkh6Kt1PTuWSdHVqjzWS5tYtsyDN/5SkBe1YlywkHSjpPknrJD0u6SupvFJtIWmipIclrU7t8G+pfKakh9I63ZpOwEDSHunxhvT8jLrXujSVr5d0SnvWKDtJu0l6VNLP0uNKtoWkTZIek7RKUn8qK/7zERFtvVE7wLwROAR4D7AamNPuehWwnicAc4G1dWXfAxal6UXAd9P06cDPqX1f41jgoVS+L/B0up+cpie3e93G2A5TgblpehLwB2qXLKlUW6T12TtNTwAeSut3GzA/lf8I+Ps0/Q/Aj9L0fODWND0nfWb2AGamz9Ju7V6/cbbJV4EfAz9LjyvZFsAmYL9hZYV/Pjphz6ASl66IiAeA54cVnwksSdNLgLPqym+MmgeBfSRNBU4BVkTE8xHxv8AK4NTia5+fiNgcEb9P09uBddS+sV6ptkjr83J6OCHdAvgYcHsqH94OQ+1zO3CSJKXyWyLitYj4I7CB2meqq0iaDpwB/Fd6LCraFiMo/PPRCWHQ6NIV09pUl1Y7ICI2Q62TBPZP5SO1SanaKu3eH0Vtq7hybZGGRVYBW6l9WDcCL0TEG2mW+nV6e33T8y8CUyhBOyT/DvwT8FZ6PIXqtkUAv5S0UrXL9kALPh9FXY5iLEa9dEUFjdQmpWkrSXsDdwCXRMRLtQ27xrM2KCtFW0TEm0CvpH2Au4DZjWZL96VtB0mfBLZGxEpJJw4VN5i19G2RHBcRz0naH1gh6cldzJtbW3TCnkGVL12xJe3Ske63pvKR2qQUbSVpArUguCki7kzFlWwLgIh4Abif2pjvPpKGNtLq1+nt9U3Pv5fasGMZ2uE4YJ6kTdSGiT9GbU+him1BRDyX7rdS20g4mhZ8PjohDKp86YplwNBR/gXA3XXln09nChwLvJh2DX8BfELS5HQ2wSdSWddIY7vXAesi4vt1T1WqLST1pD0CJP01cDK14yf3AZ9Osw1vh6H2+TTwP1E7UrgMmJ/OsJkJzAIebs1a5CMiLo2I6VG76Np8auv2WSrYFpL2kjRpaJra//VaWvH5aPeR87oj4n+gNmb6zXbXp6B1vBnYDLxOLbUvoDbOeS/wVLrfN80raj8OtBF4DOire50vUjswtgH4QrvXaxzt8LfUdlfXAKvS7fSqtQXwQeDR1A5rgX9N5YdQ68A2AD8B9kjlE9PjDen5Q+pe65upfdYDp7V73TK2y4nsPJuocm2R1nl1uj0+1B+24vPhy1GYmVlHDBOZmVmbOQzMzMxhYGZmDgMzM8NhYGZmOAzMzAyHgZmZAf8P1OS/Tz3IMnQAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.clf()\n",
"\n",
"bins = int((x_max-x_min)/7)\n",
"\n",
"# calcs = zfit.run(total_test_tf(samp))\n",
"print(total_samp[:nevents].shape)\n",
"\n",
"plt.hist(total_samp[:nevents], bins = bins, range = (x_min,x_max), label = 'data')\n",
"# plt.plot(test_q, calcs_test*nevents , label = 'pdf')\n",
"\n",
"# plt.plot(sam, calcs, '.')\n",
"# plt.plot(test_q, calcs_test)\n",
"# plt.yscale('log')\n",
"plt.ylim(0, 200)\n",
"# plt.xlim(3080, 3110)\n",
"\n",
"plt.legend()\n",
"\n",
"plt.savefig('test2.png')"
]
},
{
"cell_type": "code",
"execution_count": 30,
"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": "code",
"execution_count": 31,
"metadata": {},
"outputs": [],
"source": [
"# jpsi_width"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [],
"source": [
"# plt.hist(sample, weights=1 / prob(sample))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Fitting"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"------------------------------------------------------------------\n",
"| FCN = -8.724E+05 | Ncalls=153 (153 total) |\n",
"| EDM = 8.38E-06 (Goal: 5E-06) | up = 0.5 |\n",
"------------------------------------------------------------------\n",
"| Valid Min. | Valid Param. | Above EDM | Reached call limit |\n",
"------------------------------------------------------------------\n",
"| True | True | False | False |\n",
"------------------------------------------------------------------\n",
"| Hesse failed | Has cov. | Accurate | Pos. def. | Forced |\n",
"------------------------------------------------------------------\n",
"| False | True | True | True | False |\n",
"------------------------------------------------------------------\n",
"jpsi_p: ^{+0.016484670160195978}_{-0.016715939969484755}\n",
"psi2s_s: ^{+3.5429799872792804}_{-3.5286564682010337}\n",
"psi2s_p: ^{+0.01844073930972661}_{-0.018643934945039613}\n",
"jpsi_s: ^{+28.041395536069153}_{-27.882719703773358}\n",
"Function minimum: -872362.5577928417\n"
]
}
],
"source": [
"nll = zfit.loss.UnbinnedNLL(model=total_f, data=data2, fit_range = (x_min, x_max))\n",
"\n",
"minimizer = zfit.minimize.MinuitMinimizer(verbosity = 5)\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": 34,
"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": 35,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ8AAAD8CAYAAACo9anUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXxc1ZXg8d+pKpV2WbI2L5ItL5JlecFg2YbgEAMmGAgYgtOYhIQOpJmkoZdkOgNMuukODdNDTz6hl5AmBJIh6YkNYbPDDjHQBAi2jLzJq7xgLbYWa99VVXf+qCdZyPWkkqxaJJ3v56OPq169d++t51Id3XvPu0+MMSillFLh5Ih0A5RSSk0+GnyUUkqFnQYfpZRSYafBRymlVNhp8FFKKRV2GnyUUkqFXVDBR0TWicghESkXkfsCvB4rIs9Yr38sInkDXrvf2n5IRK4erkwRmWOVccQq0z1UHSKSJyKdIrLL+nl8tCdDKaVUeAwbfETECTwGXAMUAbeKSNGg3e4EGo0x84FHgUesY4uAjcAiYB3wUxFxDlPmI8Cjxph8oNEq27YOy1FjzDLr59sjOgNKKaXCLpiez0qg3BhzzBjTA2wG1g/aZz3wtPX4OeBKERFr+2ZjTLcx5jhQbpUXsEzrmCusMrDKvHGYOpRSSo0zriD2mQlUDHheCayy28cY4xGRZiDd2v7HQcfOtB4HKjMdaDLGeALsb1cHwBwRKQVagL81xrw/+E2IyF3AXQCJiYnLCwsLh3/nSk1y7d0ejtW3MycjkaTYob8u9lY1k5UcS3ZKXJha51fZ2Elbt4fCacn92zw+w4FTLcxIjSc90R3W9kxkO3furDfGZI5FWcEEn0C9i8Fr8tjtY7c9UI9rqP2HquMUMMsYc0ZElgMvicgiY0zLZ3Y05gngCYDi4mJTUlISoDil1EAfltfz1Sc/5um7LmbV3PQh9y164HW+tmoWP7hu8Kh8aP33Z3fzx2Nn+OC+K/q3NXf2csEP3+QH1y3kW5+fG9b2TGQi8ulYlRXMsFslkDvgeQ5QbbePiLiAKUDDEMfaba8HUq0yBtcVsA5rSO8MgDFmJ3AUKAjifSmlhuG11n50OIYf4Y6PcdLR4w11k85hjMEx6Jss1uXf0O3xhb09KjjBBJ8dQL6VhebGn0CwddA+W4HbrccbgG3Gv2LpVmCjlak2B8gHttuVaR3zjlUGVplbhqpDRDKtBAZEZK5Vx7HgT4FSyo7PGndwBDG9GhfjpLM3/MHHZ8w57XM7NfhEu2GH3az5lXuANwAn8AtjTJmIPAiUGGO2Ak8BvxaRcvw9no3WsWUi8iywH/AAdxtjvACByrSqvBfYLCIPAaVW2djVAVwGPCgiHsALfNsY0zD6U6KU6uOzok8QHR8S3E66IhJ8zg2ODofgdjro0eATtYKZ88EY8yrw6qBtDwx43AV8xebYh4GHgynT2n4Mfzbc4O0B6zDGPA88P+ybUEqNmM8adnMGM+zmjsywm88YAnXM3K7wBp/e3l4qKyvp6uoKW52hEhcXR05ODjExMSGrI6jgo5SanLz9PZ8gh90iMucTuH2xLgfdnvC1p7KykuTkZPLy8hjPV4EYYzhz5gyVlZXMmTMnZPXo8jpKKVsjmfNJcEdyzufc7eHu+XR1dZGenj6uAw+AiJCenh7yHpwGH6WULV9/ttvw+0Yq2y1QwgH4g0+4Ew7Ge+DpE473ocFHKWWrf84niC+jxFgXHd2eYfcbaz4T+MsyzuUM67CbGhkNPkopW33DbsH8JZwU66I1AsHH2Ay7RSoBItq9++67fOlLXwKgu7ubtWvXsmzZMp555pmwtkMTDpRStkaSap0c56Kt24MxJqzDT4FSrcGag9LgM6TS0lJ6e3vZtWtX2OvWno9SytZIUq2TYl0YQ9h7G3YJBwmTsOdz4sQJCgsLuf3221m6dCkbNmygo6OD119/ncLCQlavXs0LL7wAQG1tLbfddhu7du1i2bJlHD16NKxt1Z6PUsrWSFKtk+L8Xydt3R4Sh1mEdCzZzfnEu1109IR/GBDgh78rY391y/A7jkDRjBT+/vpFw+536NAhnnrqKS699FLuuOMOfvzjH/Ozn/2Mbdu2MX/+fG655RYAsrKyePLJJ/nRj37Eyy+/PKZtDYb2fJRStkxfqnWQPR+A1q7wfuHbzfkkRCj7LtJyc3O59NJLAbjtttsoKSlhzpw55OfnIyLcdtttEW6hn/Z8lFK2+hcWDXLOB/w9n3CyS7WOj+CcTzA9lFAZ3Atsbm6OyhRw7fkopWyNJNU6Kda/FEtbmHs+Pp99wkFHrxfT132bJE6ePMlHH30EwKZNm1i7di3Hjx/vn9PZtGlTJJvXT4OPUspWX7ZbsKnWAG3dvSFt02B2a7sluJ14fYYe7+RaXHThwoU8/fTTLF26lIaGBr773e/yxBNPcN1117F69Wpmz54d6SYCOuymlBpC33U+wWS79Q27hX/OJ3D74t3+9nT2eIl1OcPapkhyOBw8/vjjn9m2bt06Dh48eM6+a9asYc2aNWFq2Wdpz0cpZcs7gut8+no+7ZGY8wnwTZbo9gecyZh0MB5o8FFK2eqb8wlm2C0xNvoSDmByBZ+8vDz27dsX6WYERYOPUsqWGcGwm9vlINblCPsSO3bX+SQMGHYLl4mS3BCO96HBRyllaySp1mAtsRMt1/lYPZ/2MF1oGhcXx5kzZ8Z9AOq7n09cXFxI69GEA6WUrf5bKgR5nUhSrCsCw26B29c37Baunk9OTg6VlZXU1dWFpb5Q6ruTaShp8FFK2fKNYHkdgOS4GJo7w59qHajnk9KXfRemYBgTExPSO39ONDrsppSyNZJUa4DUhBiaOsIdfALP+aTE+S96bQlzMFTB0eCjlLI1klRrgNQEd9h7PnZzPinxVvDp0uATjTT4KKVsGWv1gGDXBkuNj6GpoyfErfosu1TrWJcDt9NBS2dkVrZWQ9Pgo5Sy5bX5YreTmuCf8+mbKwoHu4QDESEl3qU9nyilwUcpZctngltUtE9qghufCe8SOz5f4LXdwD/vo3M+0UmDj1LKlt2inXZSrXmWxjAOvXmNsU2ISI4Pf/adCo4GH6WULZ9vZMNuaYn+4NMUxi98j9c++EyJj6ElzBe9quBo8FFK2fLZrBhtZ0q8GyCsSQc+Y3DZtDElzkWr9nyikgYfpZQt7xDzKYGkJlg9nzBe6+Px2fd8UuJjNOEgSmnwUUrZMkPMpwSSlhCBns9QwcdacWG8r7c2EWnwUUrZGmmqdd+SNo3h7vnYtDE90U2v14R9pW01PA0+SilbdtfQ2HE5HaTEucKa7ebv+QT+KstI9vfE6lu7w9YeFRwNPkopW/5st5Edk5Ecy5m28AUf/5yPTVuSYgGoD2N7VHA0+CilbPlGOOcDkJkUS10YexreIXo+6Yl9wUd7PtFGg49SypbXN7JhN4DM5Fjqwvhl7x0i1bp/2E2DT9TR4KOUsmVGuMIBWMEnTD0fYwxen8FhE3ymJrgR0WG3aBRU8BGRdSJySETKReS+AK/Hisgz1usfi0jegNfut7YfEpGrhytTROZYZRyxynQPV4f1+iwRaRORvxnpSVBKBTaqYbfkWNq6PXSE4fbVfeuX2vV8XE4HUxPc2vOJQsMGHxFxAo8B1wBFwK0iUjRotzuBRmPMfOBR4BHr2CJgI7AIWAf8VEScw5T5CPCoMSYfaLTKtq1jgEeB14J940qp4XlHmO0G/jkfgPrW0Pc2PD4fMPQqDOlJbs12i0LB9HxWAuXGmGPGmB5gM7B+0D7rgaetx88BV4r/BiDrgc3GmG5jzHGg3CovYJnWMVdYZWCVeeMwdSAiNwLHgLLg37pSajh2t6geSmayP/jUtXWFoEWfZcWeIYPPtCnxnG4JfVvUyAQTfGYCFQOeV1rbAu5jjPEAzUD6EMfabU8HmqwyBtcVsA4RSQTuBX441JsQkbtEpERESurq6oZ5y0opGPnCojAg+ISht9Hf8xmijTNT46ls7Ax5W9TIBBN8Av2vDl6rwm6fsdo+VB0/xD9M1xbg9bM7GvOEMabYGFOcmZk51K5KKYt3iKVr7IQz+PTd5nuoNuakxdPQ3hOWOSgVPFcQ+1QCuQOe5wDVNvtUiogLmAI0DHNsoO31QKqIuKzezcD97epYBWwQkX8GUgGfiHQZY34SxHtTSg3B6zO4nCMLPumJsTgEalrCF3yGamNOWjwA1U2dzM9KDnmbVHCC6fnsAPKtLDQ3/gSCrYP22Qrcbj3eAGwz/pX8tgIbrUy1OUA+sN2uTOuYd6wysMrcMlQdxpjPG2PyjDF5wL8A/0sDj1JjwzPEBZx2nA5hWkoc1U2hH+rqCz5DDQ3OTPUHnwodeosqw/Z8jDEeEbkHeANwAr8wxpSJyINAiTFmK/AU8GsRKcffG9loHVsmIs8C+wEPcLcxxgsQqEyrynuBzSLyEFBqlY1dHUqp0PH6DCPs+AAwMy2eynAEH2u1artUa4CctAQAqjT4RJVght0wxrwKvDpo2wMDHncBX7E59mHg4WDKtLYfw58NN3i7bR0D9vmHoV5XSo2Mx+fDNcKeD/h7GyWfNoagRZ/l8Vo9nyGCT1ZyLG6Xg0/PtIe8PSp4usKBUsqWzzeyO5n2mZEaz+nmrv5hsVDxBdHzcTiEeZlJlNcOmZOkwkyDj1LKlsfnG3HCAfiH3Tw+Q02Ir6/xBJHtBpCflcThGg0+0USDj1LK1mhSreHsJH9ViOd9gkm1Bn/wqWrqpF1vKhc1NPgopWwNdZfQofSlN4d6kr8/1Xq44JOdBKBDb1FEg49SytZoez4zrJ5PZWPHWDfpM4JJtQYoyPZf33PgVEtI26OCp8FHKWVrNBeZAiS4XWSnxHKsPrQZZj1e//I6btfQX2V56YlMiY9hV0VTSNujgqfBRylla6i7hA5nbkYSx+pCG3x6PVbwsbuPtsXhEJblplJ6UoNPtNDgo5Sy5RnlRaYAczMTOVbXhjGhS7futa7ziRmm5wNw4axUDte20trVG7L2qOBp8FFK2Tqvnk9mEi1dHs60h+6+Pr3WsFvMMD0fgOWz0zCGsFz8qoanwUcpZcvrM8NmktmZm5kIENKht57+4DN8G1fkTSXW5eC9Q3pLlWigwUcpZcvjMzhHOe42L8Of3ny8PnTpzX09n+HmfADiYpxcMi+d9w5r8IkGGnyUUra8Pt+orvMB/yoHsS5HSFcWGMmwG8CagkyO17dzrE6v94k0DT5KKVueUV7nA/5VBwqnp1BW3TzGrTqr1xN8wgHA1YunIQJbdg2+JZkKNw0+SilbvvOY8wFYNCOF/dUtIct4G8mcD8D0KfFcOi+DF0urQpqFp4anwUcpZet85nzAH3xaujxUhmiZnZHM+fT58kUzOdnQwYdHz4SkTSo4GnyUUrbOJ9sNYNGMKQAhG3ob6ZwPwLVLppORFMvj7x0NSZtUcDT4KKVsjXZh0T6F05JxOoS9VaEKPtaczwiCT1yMkztW5/H+kXr2VOqKB5GiwUcpFZCv/3YFo/+aiItxUjgtmZ0hurCzxzOyOZ8+t108m/REN//48n6d+4kQDT5KqYD6btQ2moVFB1o5ZyqlJ5v6A8VY6vX6iHEKMsLeWUpcDH9z9QJ2nGjkpV1VY94uNTwNPkqpgIK9UdtwVuZNpdvjC8nQW4/Hh2uUPbM/Kc7lwlmpPLClLOS3flDn0uCjlArI4/P3VM5nzgegOG8qADtONJx3mwbr9viIdztHdazTIfzrLRdiDNzzm1K6er1j3Do1FA0+SqmArNhz3j2fzORY5mYm8lEIUps7e73Ex4wu+ADMSk/gR19Zyu7KJv5qc2l/b0+FngYfpVRAfT2f853zAfhCQSYfHTtDR4/nvMsaqLPXS1zM+X2NrVs8nb+7rog3ymr43rO7+tO3VWhp8FFKBTRWcz4AVxZm0+Px8WH52PZ+unq8ox52G+iO1XO4d10hW3ZV8+1f76S9e2yDpDqXBh+lVEB92W7nO+cD/oy3RLeTbYdqz7usgTp7vcS5zj/4AHxnzTweunEx7xyq5cbHPuB4iG8BPtlp8FFKBdTX83GMQc/H7XJwWUEmb+2vGdN5lc7esen59Lnt4tn86o5V1Ld1c8NP/sDLe3QB0lDR4KOUCmg066YN5YYLZlDX2s0H5fVjUh5AZ4+XuPNIOAhkdX4GW+9ZzbzMJO75TSnffWYXLXrr7TGnwUcpFdBolq4ZyuWFWaTEuXipdOwu6uw6z2w3O7lTE3ju25fw12vz2bq7mmv+5X3+eEwXIh1LGnyUUgH1jvB2BcOJi3Fy3dLpvF52mtYx6kmcb6r1UFxOB3+9toDnvn0JbpeDW3/+R/7ptQN0e/R6oLGgwUcpFVB/8AnyRm3B2LhiFh09Xn5bUjkm5XWMUbbbUC6clcbLf7GajStm8bP3jnHTYx9yuKY1pHVOBhp8lFIB9Q+7ncfCooNdkJtK8ew0fvnh8fNOPPD5DO3dHpLjXGPUOnuJsS7+6ctLePIbxdS0dPGlf/8Dv/zgeP/iq2rkNPgopQLyjPGwW587V8+hoqGT1/adOq9y2ns8+Ix/kdBwWVuUzet/fRmfn5/BD3+3n9t/uZ36tu6w1T+RaPBRSgXUE4JhN4AvLppGQXYSP37zcH+AG42WLv+FoCnxoe/5DJSZHMuTtxfz8E2L2X68gRv+/Q/srQzN/YomMg0+SqmAQjHsBv4VE75/dSHH6tt5pqRi1OW0dPqTFsLZ8+kjInxt1Wye/87nANjw+Ie8WDo281iTRVCfKhFZJyKHRKRcRO4L8HqsiDxjvf6xiOQNeO1+a/shEbl6uDJFZI5VxhGrTPdQdYjIShHZZf3sFpGbRnsylFJnnU04GNthN4C1C7NYkZfGj944xJlRDlv1B5/48AefPotnTmHrX6xmWW4q331mNz99t1xvThekYYOPiDiBx4BrgCLgVhEpGrTbnUCjMWY+8CjwiHVsEbARWASsA34qIs5hynwEeNQYkw80WmXb1gHsA4qNMcusOn4mIuHthys1AZ1NtR77ARIR4eGbltDW7eHBl/ePqoz+YbcI9HwGykiK5T+/tYr1y2bwz68f4uFXDmgiQhCC+VStBMqNMceMMT3AZmD9oH3WA09bj58DrhT/rQXXA5uNMd3GmONAuVVewDKtY66wysAq88ah6jDGdBhj+lYBjAP0f12pMdA37DZWKxwMVpCdzD2X57NlVzW/2z3yZWwa2v09ptSEyAYf8AfoR/9kGX/6uTye/MNxfvDSXu0BDSOYT9VMYODAbKW1LeA+ViBoBtKHONZuezrQNCCYDKzLrg5EZJWIlAF7gW8POL6fiNwlIiUiUlJXVxfE21Zqcuvr+YzFLRXs/Pnl8yienca9z++hvHZk187UtPiDT1ZKbCiaNmIOh/D31xfx52vmsWl7Bf/48gENQEMIJvgE+uQNPqN2+4zV9iHbYYz52BizCFgB3C8icefsaMwTxphiY0xxZmZmgKKUUgOFctitT4zTwU++ehEJbid/9qudI0pbrmnpIi0hhtgxWtV6LIgI3796AX/6uTx+8cFxHn37SKSbFLWC+VRVArkDnucAg/vI/ftY8y1TgIYhjrXbXg+kDpizGViXXR39jDEHgHZgcRDvSyk1hLFe283OtClx/OzryznV3Mk3f7mDtiDvpVPb2k12yjl/Z0aciPDAl4rYsDyHf/v9EV74RLPgAgnmU7UDyLey0Nz4Ewi2DtpnK3C79XgDsM34+5tbgY1WptocIB/Yblemdcw7VhlYZW4Zqg6rDBeAiMwGFgAngj4DSqmAxnpV66Esnz2Vn37tIvafauGbv9we1CrSlY2dTJ8SfcEH/ENw//TlJVw8dyr3Pb+XnZ82DH/QJDPsp8qaP7kHeAM4ADxrjCkTkQdF5AZrt6eAdBEpB74H3GcdWwY8C+wHXgfuNsZ47cq0yroX+J5VVrpVtm0dwGpgt4jsAl4E/twYM3Zrtis1SfV6Qj/nM9AVhdn868ZllJ5s4qs//+OQKdhen+FoXRvzs5LC0rbRiHE6+I+vLWdGahx3/WonNS1dkW5SVJHJOCFWXFxsSkpKIt0MpaLaj988xL9tK+f4P12LjMHdTIP1zsFavv2fO5mRGs/Pv1EcMMCcqG9nzY/e5ZGbl3DLillha9tolNe2cv2/f8Cy3FT+81urxuS25JEiIjuNMcVjUZaucKCUCqjXZ3A7HWENPOC/78//+9YqWjp7uemnH/BOgFtvl3zaCMCSmalhbdtozM9K5ofrF/HRsTP8x7vlkW5O1NDgo5QKqNfjG/NFRYNVnDeVLfdcSm5aAnf+3x386I1D9HjOrgP3zsFapia6KZyWHJH2jdRXluewftkMfvzWYXZVNEW6OVFBg49SKqBerw9XGJIN7OSkJfDcdy7h5oty+Mk75ax/7APe3l/D+0fqeKPsNDddOBPHOBnCEhH+8cbFZCXHce9zez4TSCcrDT5KqYB6vCbkadbDSXC7+D9fuYCff6OY5o4evvWrEr7+1HayU+K4+/L5EW3bSKXExfDQjYs5VNPK4+8djXRzIk7XQFNKBdTrjdyw22BXFWWzZkEm7x+po7G9lysXZpGa4I50s0ZsbVE2N1wwg3/fdoRrl0yP6my9UNOej1IqoG6Pj7iY6Fk9IMbp4IrCbG5enjMuA0+fB64vIs7l5KFXRreg6kShwUcpFVB3r5fYMb6RnPKvgv2XV+bz7qG6gJl8k4V+spRSAXV5fMRGUc9nIrn9c3nkpSfw8CsH+leSmGw0+CilAtKeT+i4XQ5+cF0R5bVtPHsed3Mdz/STpZQKqCvK5nwmmrULs7hwViqPbSun2+ONdHPCToOPUiog7fmElojwvasKqG7u4tmSybfytX6ylFIBRVu220S0en4GK/LSeGxbOV29k6v3o8FHKRWQ9nxCT0T47toCTrd08fwku++PfrKUUgH553z0KyLULpmXztKcKTz1/nF8vslzlwH9ZCmlAvL3fHTYLdREhG99fi7H6tvZdnDyXPejwUcpFZD2fMLn2sXTmJkazxPvH4t0U8JGP1lKqXN4vD68PqM9nzBxOR1889I8th9vYG9lc6SbExYafJRS5+iylvzXnk/4/MmKXOJjnPxm+6eRbkpY6CdLKXWObivtV3s+4ZMSF8P1F0xny65qWrt6I92ckNPgo5Q6R1/PR1Otw+urq2bT0eNl6+7qSDcl5PSTpZQ6R2ePv+cT79aeTzhdkDOFhdNT+M3HJzFmYqdda/BRSp2jvdsDQFKs3m8ynESEr67Mpay6hX1VLZFuTkhp8FFKnaO9xx98EtwafMLthmUzcTsdvFA6sVc80OCjlDpHe7d/2E17PuE3JT6Gywsz+d3uU3gm8L1+NPgopc7R0dfzidU5n0i46cKZ1Ld18+HRM5FuSsho8FFKnUN7PpG1ZkEWyXEuXtpVFemmhIwGH6XUOfoSDhI02y0i4mKcXLt4Om/sO92feTjRaPBRSp1DEw4ib/2yGbT3eHnn0MRcbFSDj1LqHO3dHuJjnDgdEummTFor50wlLSGGN8pOR7opIaHBRyl1jvYeL4mabBBRLqeDq4qy2Xawlh7PxMt60+CjlDpHe7dHh9yiwNWLptHa5eGjYxMv602Dj1LqHC2dvUyJj4l0Mya9S+dnkOh2TsihNw0+SqlzNHX2kpqgwSfS4mKcrCnM4s2yGrwT7BbbGnyUUudo7tCeT7S4etE06tu62VXRGOmmjCkNPkqpc2jPJ3p8oSATp0N491BdpJsypoIKPiKyTkQOiUi5iNwX4PVYEXnGev1jEckb8Nr91vZDInL1cGWKyByrjCNWme6h6hCRq0Rkp4jstf69YrQnQykFPp+hqaOH1Hh3pJui8K/1dmFuKu8dnmTBR0ScwGPANUARcKuIFA3a7U6g0RgzH3gUeMQ6tgjYCCwC1gE/FRHnMGU+AjxqjMkHGq2ybesA6oHrjTFLgNuBX4/sFCilBmrr8eAzaM8niqxZkMmeymbq27oj3ZQxE0zPZyVQbow5ZozpATYD6wftsx542nr8HHCliIi1fbMxptsYcxwot8oLWKZ1zBVWGVhl3jhUHcaYUmNM323/yoA4EYkN9gQopT6rucN/C2ed84keXyjIAuD9IxOn9xNM8JkJVAx4XmltC7iPMcYDNAPpQxxrtz0daLLKGFyXXR0D3QyUGmPO+fNARO4SkRIRKamrmzj/gUqNtcaOHgBSE3TYLVosmpFCRpJ7Qs37BBN8Aq2vMTjnz26fsdo+bDtEZBH+obj/FmA/jDFPGGOKjTHFmZmZgXZRSgF1rf6/3bKSdQAhWjgcwmX5mfzX4boJk3IdTPCpBHIHPM8Bqu32EREXMAVoGOJYu+31QKpVxuC67OpARHKAF4FvGGOOBvGelFI2alqs4JOiwSeafGFBJo0dveytao50U8ZEMMFnB5BvZaG58ScQbB20z1b8k/0AG4Btxhhjbd9oZarNAfKB7XZlWse8Y5WBVeaWoeoQkVTgFeB+Y8wHI3nzSqlz1bR0IQIZSRp8osnn5mUA8NEEucHcsMHHml+5B3gDOAA8a4wpE5EHReQGa7engHQRKQe+B9xnHVsGPAvsB14H7jbGeO3KtMq6F/ieVVa6VbZtHVY584G/E5Fd1k/WKM+HUpNebWsX6YluYpx6GWA0yUyOJT8racKs8xbUyoHGmFeBVwdte2DA4y7gKzbHPgw8HEyZ1vZj+LPhBm8PWIcx5iHgoWHfhFIqKDUt3WQlx0W6GSqAS+al89zOSnq9vnH/x8H4br1SasxVN3UyfYoGn2h0ydx0Onq87KlsinRTzpsGH6VUP2MMn57pYHZ6YqSbogJYNdd/dclEmPfR4KOU6lfb2k1nr5e8jIRIN0UFMDXRTeG05Akx76PBRynV73h9OwB52vOJWpfMS6fkRCPdHm+km3JeNPgopfqV17YBMCdDg0+0umRuOt0eH7srxvf1Php8lFL99lU1k5oQQ05afKSbomwU500FYOen4/v+Php8lFL99lQ2s2TmFPxr/KpoNDXRzdyMRA0+SqmJobmzl0M1rSzLTY10U9QwLpqdxicnG/EvCjM+afBRSgHwYXk9Xp/hsgJdeDfaFc9OowlRWyIAABUISURBVKG9pz9BZDzS4KOUAuD1stOkxLm4UHs+UW/57DRgfM/7aPBRSnGmrZvX9p3mpgtn4hrny7ZMBvMyk0iJc/HJSQ0+Sqlx7F9/fwSP18fXL8mLdFNUEBwOYfnsNEpOaPBRSo1Tvy2p4Fcffco3LsljflZSpJujgrR8dhpHatv6b3s+3gS1qrVSamLw+QytXR6aO3s5eLqFF0ureG3faS6Zm8791xZGunlqBC6y5n1KKxpZs2D83UVGg49SE0Bnj5eTDR18eqadT890UN3cSVNHLw3tPTR19NDU2Uuz9TMwOzc51sVfXZnP3ZfPx+3SgZDxxH89lv/aLA0+SqmQ6/X62FPZzJ7KJvZWNrO3qpnyurbPBJVEt5OpSW7SEtykJrjJy0hkSnwMqfExpMTHkJrgJjctngtnpWnQGaeS42KYm5HInsrxucyOBh+lopzPZ9h/qoUPj9bz4dEzbD/eQEePf1HJrORYluZM4dol05mXlUReegKzpyYyJSEmwq1W4bA0J5UPj9ZHuhmjosFHqSjU1evlD0fqeWt/Db8/WEN9Ww8A87OS2LA8h0vmpnPR7DSyU/Smb5PZ0pwpvFhaRU1L17j7LGjwUSpKtHT18lZZDW+Uneb9I/V09npJinXxhQWZXFmYxaXzM8bdF4wKraU5UwD/vM9VRePrs6HBR6kI6ur1su1gLVt3VbPtUC09Hh/TUuLYsDyHq4qyuXhuus7JKFtF06fgdAh7K5u4qig70s0ZEQ0+SoWZx+vjg6Nn2LKrijfLamjr9pCRFMtXV87ihmUzuDA3VVeVVkGJdzvJz0pi9zhMOtDgo1SYHK9vZ/OOkzy/s5L6th6S41xcu2QaN1wwk0vmpeN0aMBRI3dBTipvHajBGDOu/mjR4KNUCHV7vLy+7zSbt1fw0bEzOB3CFYVZbFiew5oFmcS6nJFuohrnluRM4ZmSCiobO8mdmhDp5gRNg49SIVBe28qm7RW88EkljR295KTF8/2rF7BheY4mDagxdUGOfxXyvVXNGnyUmoy6er28uvcUm7afZMeJRlwO4YuLsrl15SwunZeBQ4fVVAjkZyfhcgj7q1u4dsn0SDcnaBp8lDpPB0+3sNnq5bR0echLT+D+awq5eXkOGUmxkW6emuDiYpzMy0ziwKmWSDdlRDT4KDUKHT0eXt59ik07TlJ6sgm308G6xdO4deUsLp47dVxN/Krxb+H0ZD4+3hDpZoyIBh+lRmBfVTObtp9ky65q2ro9zM9K4m+vW8iXL8phaqI70s1Tk1TRjBRe2lVNY3sPaePkc6jBR6lhtHb1snV3NZu3V7C3qplYl4Prlk7n1pWzKJ6dpr0cFXELp6cAcOBUC5+bnxHh1gRHg49SARhj2F3ZzObtJ9m6u5qOHi+F05L54Q2LuHHZTF24U0WVvuCzX4OPUuNTc2cvW3ZVsWl7BQdOtRAf4+T6C/y9nGW68oCKUhlJsWQlx7J/HCUdaPBRk54xhp2fNrJpewWv7K2mq9fH4pkpPHTjYtYvm0FynPZyVPRbOD2FA6daI92MoGnwUZNWY3sPL5RWsXn7SY7UtpEU6+LLF+Vw64pZLLFWC1ZqvCiakcKH7x+jx+MbF4vRavBRk4oxho+PN7Bp+0le23eaHo+PZbmpPHLzEr60dAaJsforocanhdNT6PUaymvbKJqREunmDCuo3zQRWQf8K+AEnjTG/O9Br8cCvwKWA2eAW4wxJ6zX7gfuBLzAXxpj3hiqTBGZA2wGpgKfAF83xvTY1SEi6cBzwArg/xpj7hnluVAT2Jm2bp7/pJLN2ys4Vt9OcpyLW1fksnHlrP7JWqXGs6IBGW8TIviIiBN4DLgKqAR2iMhWY8z+AbvdCTQaY+aLyEbgEeAWESkCNgKLgBnA2yJSYB1jV+YjwKPGmM0i8rhV9n/Y1QF0AX8HLLZ+lAL8t5/+8OgZNu04yZtlp+n1Gopnp3H35fO5dsl04t26qKeaOOZkJBIX42D/qRZujnRjghBMz2clUG6MOQYgIpuB9cDA4LMe+Afr8XPAT8SfFrQe2GyM6QaOi0i5VR6ByhSRA8AVwFetfZ62yv0PuzqMMe3AH0Rk/gjet5rAalq6+G1JBc+UVFDR0ElqQgzfuCSPjStyyc9OjnTzlAoJp0OYn5XE4ZrxkXQQTPCZCVQMeF4JrLLbxxjjEZFmIN3a/sdBx860HgcqMx1oMsZ4AuxvV0d9EO8BEbkLuAtg1qxZwRyixhGP18e7h+rYvKOCdw7V4vUZLpmbzt98cQFXL5pGXIz2ctTEV5CdzIflZyLdjKAEE3wCXdhggtzHbnugVIyh9g+2HbaMMU8ATwAUFxcHfZyKbhUNHTyzo4Lf7qygpqWbjKRY7rpsLrcU55KXkRjp5ikVVgXZybzwSRXNnb1MiY/uSwSCCT6VQO6A5zlAtc0+lSLiAqYADcMcG2h7PZAqIi6r9zNwf7s61CTT7fHy1v4antlRwftH6nEIfKEgkwfXz+KKwixinNGfZqpUKBRkJwFwpKaV4rypEW7N0IIJPjuAfCsLrQp/AsFXB+2zFbgd+AjYAGwzxhgR2Qr8RkR+jD/hIB/Yjr8Xc06Z1jHvWGVstsrcMlQdo3vbajwqr23137qgtIqG9h5mpsbz3bUFfKU4hxmp8ZFunlIRV2DNaR6uaRv/wceaX7kHeAN/WvQvjDFlIvIgUGKM2Qo8BfzaSihowB9MsPZ7Fn9ygge42xjjBQhUplXlvcBmEXkIKLXKxq4Oq6wTQArgFpEbgS8OysZT41Rnj5dX9p5i8/aTlHzqv0HbVUXZ3LIil8/nZ+LUG7Qp1W9majyJbue4SDqQydh5KC4uNiUlJZFuhrJhjKG0oonfllTy8u5qWrs9zMlIZOOKXL58UQ6ZyXqDNqXsrH/sAxLdTn7zZxePedkistMYUzwWZenl3Cpq1LR08cInVTy3s4Kjde3ExTi4dvF0/mRFLqvm6A3alArGguwkth2si3QzhqXBR0VUt8fL2/treW5nBe8drsNnoHh2Go/cPJdrl0zXRT2VGqGC7GSeLamkob0nqm9wqMFHhZ0xhrLqFn5bUsGW3dU0dfQyLSWO76yZx80X5TA3MynSTVRq3MrvTzpo5eK56RFujT0NPipszrR182JpFc/trOTg6VbcLgdfLMrmK8W5rJ6fockDSo2BBVbwOaLBR01mPR4f7xyq5fmdlWw7WIvHZ7ggZwr/eONiblg6Q+8IqtQYy06JJTnOxaEoz3jT4KPGnDGGT0428sInVbyy9xRNHb1kJMVyx+o5bFie038tglJq7IkIC7KTOVzTFummDEmDjxozx+vbebG0ipdKqzjZ0EFcjIMvFk3jpotm8vn5Gbh05QGlwiI/O5nX953CGBO1WaIafNR5aWjv4eU91bzwSRW7KpoQgUvnZfCXV+azbvE0kvTmbEqFXUF2Epu291LX1k1WclykmxOQfjOoEevq9fL2gRpeKq3i3UN1eHyGwmnJ/M9rC7nhgplMmxKdH3alJouC/qSDNg0+anzz+gwfHz/DS6VVvLb3NK3dHrJTYrlz9RxuvHCm3g1UqSiSby0werimlUvnZ0S4NYFp8FG2+pa5+d3ual7Zc4ra1m4S3U7WLZ7Oly+aycVz0zU9WqkolJkUS2pCDEdqozfpQIOP+gxjDAdOtbJ1dzUv76mmsrETt8vB5Qsyuf6CGVxZmK23n1YqyokIBVnJHInidGsNPgqAY3VtbN1dze92V3O0rh2nQ1g9P4Pvri3gqkXZpOgyN0qNK/nZSby8J3oz3jT4TGKVjR28vOcUv9tdTVl1CyKwas5U7lg9h2sWT4/qdaGUUkMryE6mufMkda3dZKVEX9KBBp9Jprqpk9f3neaVvafY+WkjAMtyU/m7LxVx3ZLpmqmm1ARxNumgTYOPioyKhg5e33eaV/edovRkEwCF05L5/tULuH7pDGalJ0S4hUqpsVYwYIHR1fnRl/GmwWeCOlHfzmv7TvPavlPsqWwGYPHMFL5/9QKuWTxNV45WaoLLSIplaqKbI7XRmXSgwWcCOVrXxmt7T/Hq3tPsP9UCwAW5qdx/TSHXLJ6uPRylJpn5WUlRu8abBp9xzBjDwdOtvFF2mtf2nu5fxXb57DT+9rqFrFs8jZw0DThKTVYF2Uls2VUdlRlvGnzGmR6Pj+3HG3j7QA1v7a+hqqkTEViRN5V/uL6IdYs1aUAp5VeQnUxrl4ealu6o+17Q4DMONHX08O6hOt4+UMN7h+po7fYQF+Ng9fwM/uKK+VyxMCtq129SSkVOftbZpAMNPioon55p5639Nbx9oIYdJxrx+gwZSbFct3Q6axdmc+n8DF1pQCk1pIIBa7xdVpAZ4dZ8lgafKOH1GXZVNPH2gRre3l/TvybTguxkvv2FuaxdmM0FOak4dC01pVSQ0pNiSU90cyQKkw40+ERQR4+H94/U8/sDNWw7WEt9Ww8uh7Bq7lRuXTmLtQuzNUNNKXVe8rOTOByF6dYafMKspqWL3x+o5e0DNfyhvJ4ej4/kOBeXL8hibVE2XyjIZEq8rqOmlBobBdnJvPhJVdRlvGnwCTFjDEdq23iz7DRv7a9ht3XBZ+7UeL62ahZXLcxmxZypxOgtppVSIZCfnUxrt4dTzV3MSI2PdHP6afAJAZ/P8MnJRt7cX8ObZac5caYD8K+h9v2rF7B2YTYF2UlR9VeIUmpiKsg6m3SgwWeCKq9t5flPqniptIpTzV3EOIXPzcvgzy7zJwxkR+Hifkqpia1vjbfy2jbWLMiKcGvO0uBznowxfFB+hsffO8ofyutxOoTL8jO475pCLi/M0vvgKKUiKi3RTUZSLIej7MZyGnzOQ2VjBw9sKWPbwVqykmO5d10hG5bnkJkcG+mmKaVUv4Ls6FvjTYPPKG0/3sB/+3UJPR4ff3vdQr5+yWxiXXrRp1Iq+hRkJ/PczsqoynjT4DMK+6tb+OYvt5M9JY6nbl/BnIzESDdJKaVs5Wcn0dbtobq5i5lRknSg+b0j5PH6+O+/3U1irIvffOtiDTxKqag3cI23aKHBZ4Re3XeaA6da+PvrF0XdQn1KKRXIgmn+4FNW1RzhlpylwWeEnt9ZSU5aPNcsnhbppiilVFCmxMdQkJ1EyaeNkW5Kv6CCj4isE5FDIlIuIvcFeD1WRJ6xXv9YRPIGvHa/tf2QiFw9XJkiMscq44hVpnu0dYw1j9dHyYkGrijM0gU+lVLjSnHeVHZ+6l8hPxoMG3xExAk8BlwDFAG3ikjRoN3uBBqNMfOBR4FHrGOLgI3AImAd8FMRcQ5T5iPAo8aYfKDRKnvEdYz0RASjprWb9h4vhdNSQlG8UkqFzIq8NFq7PJRVR8fQWzA9n5VAuTHmmDGmB9gMrB+0z3rgaevxc8CV4s/nWw9sNsZ0G2OOA+VWeQHLtI65wioDq8wbR1nHmGvq6AFgaqI7FMUrpVTIrCnIwu1y8JNt5ZFuChBcqvVMoGLA80pgld0+xhiPiDQD6db2Pw46dqb1OFCZ6UCTMcYTYP/R1NFPRO4C7rKetonIGaDe9l0P4ZpHRnNUVMtglOdiAtJz4afn4awJdS6OAD+/fVSHZgCzx6odwQSfQJMbgwcN7fax2x6oxzXU/qOp47MbjHkCeKLvuYiUGGOKAxw76ei5OEvPhZ+eh7P0XPhZ5yFvrMoLZtitEsgd8DwHqLbbR0RcwBSgYYhj7bbXA6lWGYPrGmkdSimlolQwwWcHkG9lobnxT+5vHbTPVqCvI7cB2GaMMdb2jVam2hwgH9huV6Z1zDtWGVhlbhllHUoppaLUsMNu1vzKPcAbgBP4hTGmTEQeBEqMMVuBp4Bfi0g5/t7IRuvYMhF5FtgPeIC7jTFegEBlWlXeC2wWkYeAUqtsRlPHMJ4YfpdJQ8/FWXou/PQ8nKXnwm9Mz4P4Ow9KKaVU+OgKB0oppcJOg49SSqmwm5TBZ7jlgiYCEfmFiNSKyL4B26aKyFvW0kVviUiatV1E5N+s87FHRC4acMzt1v5HRGR0VwdEkIjkisg7InJARMpE5K+s7ZPqXIhInIhsF5Hd1nn4obU9apezCjVrtZVSEXnZej4pz4WInBCRvSKyS0RKrG2h//0wxkyqH/wJDkeBuYAb2A0URbpdIXiflwEXAfsGbPtn4D7r8X3AI9bja4HX8F8zdTHwsbV9KnDM+jfNepwW6fc2wvMwHbjIepwMHMa/pNOkOhfW+0myHscAH1vv71lgo7X9ceA71uM/Bx63Hm8EnrEeF1m/M7HAHOt3yRnp9zfKc/I94DfAy9bzSXkugBNAxqBtIf/9mIw9n2CWCxr3jDH/hT8rcKCBSxQNXrroV8bvj/ivtZoOXA28ZYxpMMY0Am/hXz9v3DDGnDLGfGI9bgUO4F8BY1KdC+v99N1HOcb6MUTxclahJCI5wHXAk9bzqF7aKwJC/vsxGYNPoOWCzlmOZ4LKNsacAv+XMpBlbbc7JxPqXFnDJRfi/6t/0p0La5hpF1CL/8vhKEEuZwUMXM5qXJ8Hy78A/wPwWc+DXtqLiXcuDPCmiOwU/zJkEIbfj8l4G+2gluOZZM5r6aLxQESSgOeBvzbGtIj9fewn7Lkw/uvflolIKvAisDDQbta/E/Y8iMiXgFpjzE4RWdO3OcCuE/5cWC41xlSLSBbwlogcHGLfMTsXk7HnM5mX46mxushY/9Za20e6DNK4IiIx+APP/zPGvGBtnpTnAsAY0wS8i3/MfjIuZ3UpcIOInMA/7H4F/p7QZDwXGGOqrX9r8f9RspIw/H5MxuATzHJBE9XAJYoGL130DSuT5WKg2epqvwF8UUTSrGyXL1rbxg1rbP4p4IAx5scDXppU50JEMq0eDyISD6zFP/816ZazMsbcb4zJMf5FMjfif29fYxKeCxFJFJHkvsf4P9f7CMfvR6QzLSLxgz9j4zD+Me8fRLo9IXqPm4BTQC/+v0ruxD9O/Xv8q6r/Hphq7Sv4b+53FNgLFA8o5w78E6nlwDcj/b5GcR5W4+/+7wF2WT/XTrZzASzFv1zVHuvL5QFr+1z8X5jlwG+BWGt7nPW83Hp97oCyfmCdn0PANZF+b+d5XtZwNttt0p0L6z3vtn7K+r4Pw/H7ocvrKKWUCrvJOOymlFIqwjT4KKWUCjsNPkoppcJOg49SSqmw0+CjlFIq7DT4KKWUCjsNPkoppcLu/wNIdT85YgImXQAAAABJRU5ErkJggg==\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, 5e-6)\n",
"# plt.yscale('log')\n",
"# plt.xlim(3080, 3110)\n",
"plt.savefig('test3.png')\n",
"# print(jpsi_width)"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [],
"source": [
"# _tot = 4.37e-7+6.02e-5+4.97e-6\n",
"# _probs = []\n",
"# _probs.append(6.02e-5/_tot)\n",
"# _probs.append(4.97e-6/_tot)\n",
"# _probs.append(4.37e-7/_tot)\n",
"# print(_probs)"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [],
"source": [
"# dtype = 'float64'\n",
"# # mixture = tfd.Uniform(tf.constant(x_min, dtype=dtype), tf.constant(x_max, dtype=dtype))\n",
"# mixture = tfd.MixtureSameFamily(mixture_distribution=tfd.Categorical(probs=[tf.constant(0.007, dtype=dtype),\n",
"# tf.constant(0.917, dtype=dtype),\n",
"# tf.constant(0.076, dtype=dtype)]),\n",
"# components_distribution=tfd.Uniform(low=[tf.constant(x_min, dtype=dtype), \n",
"# tf.constant(3080, dtype=dtype),\n",
"# tf.constant(3670, dtype=dtype)], \n",
"# high=[tf.constant(x_max, dtype=dtype),\n",
"# tf.constant(3112, dtype=dtype), \n",
"# tf.constant(3702, dtype=dtype)]))\n",
"# # for i in range(10):\n",
"# # print(zfit.run(mixture.prob(mixture.sample((10, 1)))))"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1.2311725935705207\n",
"1.4004132078629088\n"
]
}
],
"source": [
"print((zfit.run(jpsi_p)%(2*np.pi))/np.pi)\n",
"print((zfit.run(psi2s_p)%(2*np.pi))/np.pi)"
]
},
{
"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
}
