{
"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": [
"## C_q,qbar constraint"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"\n",
"\n",
"\n",
"# r = rho_scale * rho_width/rho_mass * np.cos(rho_phase)*(1-np.tan(rho_phase)*rho_width/rho_mass)\n",
"# o = omega_scale*np.cos(omega_phase)*omega_width/omega_mass\n",
"# p = phi_scale*np.cos(phi_phase)*phi_width/phi_mass\n",
"\n",
"# phi_s = np.linspace(-500, 5000, 100000)\n",
"\n",
"# p_ = phi_s*np.cos(phi_phase)*phi_width/phi_mass\n",
"\n",
"# p_y = r+o+p_\n",
"\n",
"# plt.plot(phi_s, p_y)\n",
"\n",
"# # print(r + o + p)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Build pdf"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"class total_pdf(zfit.pdf.ZPDF):\n",
" _N_OBS = 1 # dimension, can be omitted\n",
" _PARAMS = ['rho_mass', 'rho_scale', 'rho_phase', 'rho_width',\n",
" 'jpsi_mass', 'jpsi_scale', 'jpsi_phase', 'jpsi_width',\n",
" 'psi2s_mass', 'psi2s_scale', 'psi2s_phase', 'psi2s_width',\n",
" 'p3770_mass', 'p3770_scale', 'p3770_phase', 'p3770_width',\n",
" 'p4040_mass', 'p4040_scale', 'p4040_phase', 'p4040_width',\n",
" 'p4160_mass', 'p4160_scale', 'p4160_phase', 'p4160_width',\n",
" 'p4415_mass', 'p4415_scale', 'p4415_phase', 'p4415_width',\n",
" 'omega_mass', 'omega_scale', 'omega_phase', 'omega_width',\n",
" 'phi_mass', 'phi_scale', 'phi_phase', 'phi_width'] # the name of the parameters\n",
"\n",
" def _unnormalized_pdf(self, x):\n",
" \n",
" x = x.unstack_x()\n",
" \n",
" def rho_res(q):\n",
" return resonance(q, _mass = self.params['rho_mass'], scale = self.params['rho_scale'],\n",
" phase = self.params['rho_phase'], width = self.params['rho_width'])\n",
" \n",
" def omega_res(q):\n",
" return resonance(q, _mass = self.params['omega_mass'], scale = self.params['omega_scale'],\n",
" phase = self.params['omega_phase'], width = self.params['omega_width'])\n",
" \n",
" def phi_res(q):\n",
" return resonance(q, _mass = self.params['phi_mass'], scale = self.params['phi_scale'],\n",
" phase = self.params['phi_phase'], width = self.params['phi_width'])\n",
"\n",
" def jpsi_res(q):\n",
" return resonance(q, _mass = self.params['jpsi_mass'], scale = self.params['jpsi_scale'],\n",
" 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'],\n",
" phase = self.params['psi2s_phase'], width = self.params['psi2s_width'])\n",
" \n",
" def p3770_res(q):\n",
" return resonance(q, _mass = self.params['p3770_mass'], scale = self.params['p3770_scale'],\n",
" phase = self.params['p3770_phase'], width = self.params['p3770_width'])\n",
" \n",
" def p4040_res(q):\n",
" return resonance(q, _mass = self.params['p4040_mass'], scale = self.params['p4040_scale'],\n",
" phase = self.params['p4040_phase'], width = self.params['p4040_width'])\n",
" \n",
" def p4160_res(q):\n",
" return resonance(q, _mass = self.params['p4160_mass'], scale = self.params['p4160_scale'],\n",
" phase = self.params['p4160_phase'], width = self.params['p4160_width'])\n",
" \n",
" def p4415_res(q):\n",
" return resonance(q, _mass = self.params['p4415_mass'], scale = self.params['p4415_scale'],\n",
" phase = self.params['p4415_phase'], width = self.params['p4415_width'])\n",
" \n",
"\n",
" funcs = rho_res(x) + omega_res(x) + phi_res(x) + jpsi_res(x) + psi2s_res(x) + p3770_res(x) + p4040_res(x)+ p4160_res(x) + p4415_res(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": 7,
"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": 8,
"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": [
"#rho\n",
"\n",
"rho_mass, rho_width, rho_phase, rho_scale = pdg[\"rho\"]\n",
"\n",
"rho_m = zfit.Parameter(\"rho_m\", ztf.constant(rho_mass), floating = False)\n",
"rho_w = zfit.Parameter(\"rho_w\", ztf.constant(rho_width), floating = False)\n",
"rho_p = zfit.Parameter(\"rho_p\", ztf.constant(rho_phase), floating = False)\n",
"rho_s = zfit.Parameter(\"rho_s\", ztf.constant(rho_scale))\n",
"\n",
"#omega\n",
"\n",
"omega_mass, omega_width, omega_phase, omega_scale = pdg[\"omega\"]\n",
"\n",
"omega_m = zfit.Parameter(\"omega_m\", ztf.constant(omega_mass), floating = False)\n",
"omega_w = zfit.Parameter(\"omega_w\", ztf.constant(omega_width), floating = False)\n",
"omega_p = zfit.Parameter(\"omega_p\", ztf.constant(omega_phase), floating = False)\n",
"omega_s = zfit.Parameter(\"omega_s\", ztf.constant(omega_scale))\n",
"\n",
"\n",
"#phi\n",
"\n",
"phi_mass, phi_width, phi_phase, phi_scale = pdg[\"phi\"]\n",
"\n",
"phi_m = zfit.Parameter(\"phi_m\", ztf.constant(phi_mass), floating = False)\n",
"phi_w = zfit.Parameter(\"phi_w\", ztf.constant(phi_width), floating = False)\n",
"phi_p = zfit.Parameter(\"phi_p\", ztf.constant(phi_phase), floating = False)\n",
"phi_s = zfit.Parameter(\"phi_s\", ztf.constant(phi_scale))\n",
"\n",
"#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), floating = False)\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), floating = False)\n",
"psi2s_s = zfit.Parameter(\"psi2s_s\", ztf.constant(psi2s_scale))\n",
"\n",
"#psi(3770)\n",
"\n",
"p3770_mass, p3770_width, p3770_phase, p3770_scale = pdg[\"p3770\"]\n",
"\n",
"p3770_m = zfit.Parameter(\"p3770_m\", ztf.constant(p3770_mass), floating = False)\n",
"p3770_w = zfit.Parameter(\"p3770_w\", ztf.constant(p3770_width), floating = False)\n",
"p3770_p = zfit.Parameter(\"p3770_p\", ztf.constant(p3770_phase), floating = False)\n",
"p3770_s = zfit.Parameter(\"p3770_s\", ztf.constant(p3770_scale), floating = False)\n",
"\n",
"#psi(4040)\n",
"\n",
"p4040_mass, p4040_width, p4040_phase, p4040_scale = pdg[\"p4040\"]\n",
"\n",
"p4040_m = zfit.Parameter(\"p4040_m\", ztf.constant(p4040_mass), floating = False)\n",
"p4040_w = zfit.Parameter(\"p4040_w\", ztf.constant(p4040_width), floating = False)\n",
"p4040_p = zfit.Parameter(\"p4040_p\", ztf.constant(p4040_phase), floating = False)\n",
"p4040_s = zfit.Parameter(\"p4040_s\", ztf.constant(p4040_scale), floating = False)\n",
"\n",
"#psi(4160)\n",
"\n",
"p4160_mass, p4160_width, p4160_phase, p4160_scale = pdg[\"p4160\"]\n",
"\n",
"p4160_m = zfit.Parameter(\"p4160_m\", ztf.constant(p4160_mass), floating = False)\n",
"p4160_w = zfit.Parameter(\"p4160_w\", ztf.constant(p4160_width), floating = False)\n",
"p4160_p = zfit.Parameter(\"p4160_p\", ztf.constant(p4160_phase), floating = False)\n",
"p4160_s = zfit.Parameter(\"p4160_s\", ztf.constant(p4160_scale), floating = False)\n",
"\n",
"#psi(4415)\n",
"\n",
"p4415_mass, p4415_width, p4415_phase, p4415_scale = pdg[\"p4415\"]\n",
"\n",
"p4415_m = zfit.Parameter(\"p4415_m\", ztf.constant(p4415_mass), floating = False)\n",
"p4415_w = zfit.Parameter(\"p4415_w\", ztf.constant(p4415_width), floating = False)\n",
"p4415_p = zfit.Parameter(\"p4415_p\", ztf.constant(p4415_phase), floating = False)\n",
"p4415_s = zfit.Parameter(\"p4415_s\", ztf.constant(p4415_scale))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup pdf"
]
},
{
"cell_type": "code",
"execution_count": 9,
"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",
" p3770_mass = p3770_m, p3770_scale = p3770_s, p3770_phase = p3770_p, p3770_width = p3770_w,\n",
" p4040_mass = p4040_m, p4040_scale = p4040_s, p4040_phase = p4040_p, p4040_width = p4040_w,\n",
" p4160_mass = p4160_m, p4160_scale = p4160_s, p4160_phase = p4160_p, p4160_width = p4160_w,\n",
" p4415_mass = p4415_m, p4415_scale = p4415_s, p4415_phase = p4415_p, p4415_width = p4415_w,\n",
" rho_mass = rho_m, rho_scale = rho_s, rho_phase = rho_p, rho_width = rho_w,\n",
" omega_mass = omega_m, omega_scale = omega_s, omega_phase = omega_p, omega_width = omega_w,\n",
" phi_mass = phi_m, phi_scale = phi_s, phi_phase = phi_p, phi_width = phi_w) \n",
" \n",
"# print(total_pdf.obs)\n",
"\n",
"# print(calcs_test)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Test if graphs actually work and compute values"
]
},
{
"cell_type": "code",
"execution_count": 10,
"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": 23,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ8AAAD8CAYAAACo9anUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJztnXl8nNV1979nZrRLlixZ8iLJloxlvNsYYQgQwo5ZTQIkTkJKlre0KfR9G9q3QLM0oeFtSdOQtCUhJCShJMEQIMEBh30JCYsXbLyBsSxvso0la99nu+8f88xoNJpNsjSbzvfz0cczd+5z7p3H0vzmnHvuuWKMQVEURVESiS3ZE1AURVEmHyo+iqIoSsJR8VEURVESjoqPoiiKknBUfBRFUZSEo+KjKIqiJJy4xEdEVovIHhFpEJE7wryeIyKPWq+/LSI1Qa/dabXvEZHLYtkUkVrLxl7LZnYcYywTkTdFZJeI7BCR3LHcDEVRFCUxxBQfEbED9wGXA4uAT4vIopBuXwLajTHzgHuBe6xrFwFrgcXAauCHImKPYfMe4F5jTB3QbtmONoYD+CXw18aYxcD5gGuU90FRFEVJIPF4PquABmNMozHGCawD1oT0WQM8ZD1+HLhIRMRqX2eMGTTG7AcaLHthbVrXXGjZwLJ5bYwxLgW2G2PeBTDGtBpjPPHfAkVRFCXROOLoUwkcDnreBJwZqY8xxi0inUCZ1f5WyLWV1uNwNsuADmOMO0z/SGPMB4yIPAeU4xO774S+CRG5GbgZoKCg4PQFCxbE8dYVZfJyomeQY50DLJ41BZvIiNc/7BqgpXuQpZXFSZjdEDuOdFJRlMP0Kb5oe2NLLwBzywuSOa2MZMuWLSeMMeXjYSse8Rn5WwehNXki9YnUHs7jitY/2hgO4FzgDKAPeElEthhjXhrW0ZgHgAcA6uvrzebNm8OYUxTFz49e3cc9z77PG/+ymtws+4jXv/f8Hv7z5QY2/esVSBhxSgQer+GUf9rAbZfM539fVAfA2gfexGvgsb/6SFLmlMmIyMHxshVP2K0JqA56XgUcjdTHWoMpBtqiXBup/QRQYtkIHSvaGK8ZY04YY/qADcDKON6XoihR8Hi9ANht4YXFZrV7k1ge0mvVpgyeok0ErVmZ+sQjPpuAOisLLRtfAsH6kD7rgZusx9cDLxvf//56YK2VqVYL1AEbI9m0rnnFsoFl86kYYzwHLBORfEuUPgbsjv8WKIoSDrelKvYIXo3D+sT3JFF9/OIT7HnZRJIqiEp8xAy7Wesrt+L7kLcDPzPG7BKRu4DNxpj1wIPAwyLSgM8bWWtdu0tEHsMnBm7gFn8yQDib1pC3A+tE5NvAVss2UcZoF5Hv4RM0A2wwxjxzUndFURTcHoNNhjycUGwpID5+Byd4TUpkSJSU1CWeNR+MMRvwhbOC274R9HgAuCHCtXcDd8dj02pvxJcNF9oebYxf4ku3VhRlnHB7DQ5b5OCI3yPyJPGDPlLYLRl66HK5aGpqYmBgIPGDjzO5ublUVVWRlZU1YWPEJT6Kokw+PF4vDnvkRAJ7Cng+3jCej01IyppPU1MTRUVF1NTUJC0BYzwwxtDa2kpTUxO1tbUTNo6W11EUJSxur4mYbABD4uNNovj4hU9CPJ9kCOLAwABlZWVpLTzgWz8rKyubcA9OxUdRlLB4vCaQVBAOv/i4k7rmY4bNBXxrUcmaUroLj59EvA8VH0VRwuLyGOxR1nz8oa5kLu6nUthNGR0qPgoer6G5O/0XSZXxxeP1RvV8UinVemTCgYpPJF599VWuuuoqAAYHB7n44otZsWIFjz76aELnoQkHCt959n1+/MdGNn/tYqYV5iR7OkqK4PaaqAkHqZBqHRAZ3eczJrZu3YrL5WLbtm0JH1s9H4UX3zsOQEefM8kzUVIJtydGwoGkgPh4h88FJvc+nwMHDrBgwQJuuukmli1bxvXXX09fXx/PPvssCxYs4Nxzz+XJJ58EoLm5mRtvvJFt27axYsUK9u3bl9C5quejKEpY3F4v2fbI30/9XlEy9/m4LfUJDg/6yuska0Y+vvX7Xew+2jWuNhfNmsI/X704Zr89e/bw4IMPcs455/DFL36R733ve/z4xz/m5ZdfZt68eXzqU58CoKKigp/+9Kd897vf5emnnx7XucaDej6KooTF6TY4oohPIOEgBTwfm214wkEyvbFkU11dzTnnnAPAjTfeyObNm6mtraWurg4R4cYbb0zyDH2o56OMKFGuKAAuj5fseDaZJtHN8ARSrYfafKnWyf2tjsdDmShC06Q7OztTMgVcPR9FUcLi9nrJisPzcXuSv8nUJqkVdksmhw4d4s033wTgkUce4eKLL2b//v2BNZ1HHnkkmdMLoOKjKEpYXG4TVXwctlTY5xNmk+kkTjgAWLhwIQ899BDLli2jra2Nr3zlKzzwwANceeWVnHvuucyZMyfZUwQ07KYoSgScHi9FWZE/IlKhtpsnzLEPk32fj81m4/777x/Wtnr1at5///0Rfc8//3zOP//8BM1sOOr5KIoSlljZbrYU8HwCYTdbcKq17vNJB1R8FEUJS7xht2Su+QTCblpeB4Camhp27tyZ7GnEhYqPoihhcXm8ZDliJxwkNdvNG27NJ3meT6aIXiLeh4qPEkTqpWMqycPp8ZIV15EKiZrRSMKLT3LWoXJzc2ltbU17AfKf55Obmzuh42jCgaIoYXF7oofd/C+lnOeTpH0+VVVVNDU10dLSkvCxxxv/SaYTiYqPoihh8YXdonk+PvXxJNH18ZjU2eeTlZU1oSd/ZhoadlMUJSxOjxdHlPN8hgqLJmpGIwkUFtV9PmmHio+iKGFxebxkR0s48IfdkrnPJ1x5nUm+zyddUPFRFCUsvjWfaIfJ+cNuySwsOjLspvt80gMVH0VRRuD1GtzeNE04mMT7fNIJFR9FUUbgshZT4iksmswjFSIlHEzmIxXSBRUfRVFG4LKqFkQLu6VCbTdvxFTrZM1IiRcVH0UP9FFG4PbE9nxSQXz8no8jJOwGGnpLdVR8FEUZgTMO8fG/5krmPp8whUUD4UDVnpQmLvERkdUiskdEGkTkjjCv54jIo9brb4tITdBrd1rte0Tkslg2RaTWsrHXspkdbQwRqRGRfhHZZv0MryWuKMqoiSfslgqFRcMfqeD7V9OtU5uY4iMiduA+4HJgEfBpEVkU0u1LQLsxZh5wL3CPde0iYC2wGFgN/FBE7DFs3gPca4ypA9ot2xHHsNhnjFlh/fz1qO6AoigjcLljez4Ov+eTxF2m4bLdJOD5qPikMvF4PquABmNMozHGCawD1oT0WQM8ZD1+HLhIfL8Ba4B1xphBY8x+oMGyF9amdc2Flg0sm9fGGEM5WfQuKiG448h283tF7mQmHJjIYTfVntQmHvGpBA4HPW+y2sL2Mca4gU6gLMq1kdrLgA7LRuhYkcYAqBWRrSLymoh8NI73pChKFJzueMJuvo8Pd1I9H9+/GnZLP+IpLBruty/0fzVSn0jt4UQvWv9oYxwDZhtjWkXkdOB3IrLYGNM1bIIiNwM3A8yePTuMKUVR/LjiSjgQq2/ys91sIeV1ILlZeEps4vF8moDqoOdVwNFIfUTEARQDbVGujdR+AiixbISOFXYMK6TXCmCM2QLsA+aHvgljzAPGmHpjTH15eXkcb3sSoX+jSgjxiI+IYLdJIESXDLzhEg5smu2WDsQjPpuAOisLLRtfAsH6kD7rgZusx9cDLxtfkv16YK2VqVYL1AEbI9m0rnnFsoFl86loY4hIuZXAgIjMtcZojP8WKIoSitNKOMiJUlgUfBlvSfV8IpTXAd3nk+rEDLsZY9wicivwHGAHfmaM2SUidwGbjTHrgQeBh0WkAZ/Hs9a6dpeIPAbsBtzALcYYD0A4m9aQtwPrROTbwFbLNpHGAM4D7hIRN+AB/toY0zb2W6IoyqBffLLsUftl221JzXaLlnCgnk9qE9dhcsaYDcCGkLZvBD0eAG6IcO3dwN3x2LTaG/Flw4W2hx3DGPME8ETMN6EoStwMuj1AHJ6PXZK6z8efaReuwoGu+aQ2WuFAUZQRDMYbdrPbkrrm4wlzpIL/hFXNdkttVHwURRnBoCu+sFtWktd8whUWDVReUM8npVHxURRlBPGH3WzJ3edjRma7BQqeJlEUldio+CgBtF6E4if+sJvgSmZVa69BZHjCgT3g+SRPFJXYqPgoijKCIfGJFXZLrufj8hiybMM/xuyBfT7q+aQyKj6Kooxg0OVBJHp5HYAsR3LXfNweL46QOYau+RhjkpoOroRHxUdRlBEMur3kOGyBCtGRcNiSu8/H7TXD0qwhKOxmieJ3nttD3Vf/kPC5KdFR8VEUZQQ+8YkecgOfZ5TMfT4uj3dECaDQE1Z/9Oo+YCiJYqL51z+8x/de+CAhY6UzKj6Kooxg0O2JmWwAPs8nmQv7bo8ZEXazR0i17h2cePHxeg0/fq2R/3xpr5b3iYGKj6IoIxh0ecnJikN87Mld83F5vYGjHfw4Imwy7R10M9Ec6egPPG7pGZzw8dIZFR9Fi1orI4g/7JZ8zyc0KSJ0zcdPn3PiPZ/WXmfgcXOXik80VHwURRlB/GG3ZNd28waO8/YTuubjp9c58Z5PsHcVLETKSFR8FEUZgT/bLRZZSa5q7fJEyXYL8cj6E+D5dA8Mic+JbvV8oqHiowTQ9VHFz6Ar/my3ZO/zCc12c0TYZJoI8Rnu+aj4REPFR1GUEQy6PXEmHCS3woHbGyXbLUQU+10TLz49QeLT0eea8PHSGRUfJQh1fRQf8YfdklvbzeXxjiiv4xcj/5pPtuUZJVJ8CnMcdA2o+ERDxUdRlBEMuDxkxxF2cyS5tpvbY4YdpwBDFa79+3z82XCJCLv1DLrJsgvTCrPp6p/4BId0RsVHCaBrPoqfPqeH/Bhn+UDyTzJ1RQm7+dd8shwJ9HwG3BTmOJiSl6WeTwxUfBRFGUG/y0Nednz7fJzJPM/HGy7hwPfcL4r+BIREJRwU5jqYkptFV7+KTzRUfJQA6vgofvqdHvLjEJ8chy/VOlmlZNzhUq1D1nz8S1KJ8Hy6B90UZDsozsuiU8UnKio+iqIMw+n24vaauMXHa5J3ZHXYwqIhaz7+fUgJWfMZcFOU62BKnoOuAV3ziYaKj6Iow/B/SOdlO2L2zbbWU/yHzyWaaKnW/iO2/R5QIjyfXqebghwNu8WDio8SQBMOFIA+l+8be3yej6/PYAI+2MPhC7uF32TqsTwedwLFJzjhYNDtZSBJ9yUdUPFRFGUY/gKc8YbdIHmejy/sFn7Nxy86Ac8nQanWRbkOpuT6vMZuDb1FRMVHCWA05UAhKOwWR6q1vwqCM5XCbjKUcGCMSbj4FGT7PB9A062joOKjKMowhjyf2Gs+gbBbEj2f0LBb8JpPcCLERIfdPF5Dn9PjS7W2xEcz3iKj4qME0DUfBaDPOnognn0+/tI1iTqiOpRB98hD74bWfMywYxUmev3Ff2RDYY4v1RpUfKKh4qMoyjD6R7Pmk5W8NR9jDM4wh94FH6Md7PlM9GFyPQNhxEeLi0YkLvERkdUiskdEGkTkjjCv54jIo9brb4tITdBrd1rte0Tkslg2RaTWsrHXspkdawzr9dki0iMi/zDam6D4UM9HgdEmHPiz3RIvPv7KCqEFUEUEm/jCYB5P4sJu/uMUCnMdlKjnE5OY4iMiduA+4HJgEfBpEVkU0u1LQLsxZh5wL3CPde0iYC2wGFgN/FBE7DFs3gPca4ypA9ot2xHHCOJe4A/xvnFFUcLT5/Lv84k/283pSXzYzZ/kEK76tsNms9Z8fH2y7MLABHs+3Zb4FOTomk88xOP5rAIajDGNxhgnsA5YE9JnDfCQ9fhx4CIREat9nTFm0BizH2iw7IW1aV1zoWUDy+a1McZARK4FGoFd8b91xU+ySqMoqUm/07/PJ46EA3/YLQmejz/Ulx1GfOw2wRMUdivKzaLP5ZnQ33W/51OU4yDLbqMg265n+kQhHvGpBA4HPW+y2sL2Mca4gU6gLMq1kdrLgA7LRuhYYccQkQLgduBb0d6EiNwsIptFZHNLS0uMtzw50VRrBYbCbvGkWg8lHCQh7BbV8xFcHm9AfApzHHi8ZkJPXfWv+RTk+ES7JD9bPZ8oxCM+EqYt9H8wUp/xao82xrfwhel6wrw+1NGYB4wx9caY+vLy8mhdFWVS0+f0kJtlG3FOTjhysvyp1okPuw0GxGekSDrsMmzNp8ja9DmR6z7dQQfJAUzR4qJRie1X+7yP6qDnVcDRCH2aRMQBFANtMa4N134CKBERh+XdBPePNMaZwPUi8h2gBPCKyIAx5r/jeG9KEBp9UwC6B1wU5WbF1Tew5pMEz8cveOHCbll2m1Ug1TcvvyAMuDyBTLTxJhB2s4SuOM9BZ79zQsbKBOLxfDYBdVYWWja+BIL1IX3WAzdZj68HXja+4Op6YK2VqVYL1AEbI9m0rnnFsoFl86loYxhjPmqMqTHG1ADfB/6fCo+ijJ0uqzJzPCSzvE60sJv/nCGPd7jnM5Hp1qFhNz1WIToxf8OMMW4RuRV4DrADPzPG7BKRu4DNxpj1wIPAwyLSgM8bWWtdu0tEHgN2A27gFmOMByCcTWvI24F1IvJtYKtlm0hjKIoyvnQPuOP2fJJZ1Tpa2M13ztDQGo///UxkiZ0ep5tshy1wxENJXjad/R0TNl66E9fXG2PMBmBDSNs3gh4PADdEuPZu4O54bFrtjfiy4ULbI44R1Oeb0V5XFCU23QOuQGHMWAQSDpJQvdkZJdsty27D5R7p+Uzoms+Ae9h9K87P0my3KGiFAyWArvko4Pd84hMfESHHYUuS5+MTkrBhN4fg9Ixc85lQz8c6TsFPsR6rEBUVH0VRhtE94Br2IRqLZIlPTM8naM2nMAGeT8+gOzAOEEhs0EPlwqPioyjKMHpGseYDkJtlT8hxBaEMRkk4yA5ku1lht5wEiM+Am6KcofvmF58OFZ+wqPgoAXSTqeLxGnqdnrjDbuDL7upLQmjJX1UhJ8xm2GyHz/NxhyQcTGSJne4Ino9mvIVHxUdRlAD+dOHReD752Xb6BhN/YmfUNR8r1Tp0zcd/XMRE0DPoCnhYACX5luejSQdhUfFRAmjCgeI/eXM0nk9+tn3CjysIR7Tq29l2Gy63CbPmM3FrUz0Dwz2fqfnZALT1Dk7YmOmMio+iKAG6Lc8n3lRr8BUgnUiPIhK9Tg8ikBtmn0+WP+xmiU9B9sSu+RhjfAkHQZ5PeVEOACd6tMpBOFR8lADq+ChDns8ow25J8Hz6nW7ysuzYwtSgy7ILg0H7fLIcQl6WPVCxe7wZdHtxecwwzyc3y05RjoOWbvV8wqHioyhKgI4+37d0f8goHnyeT+LFp9fpiXjgXXZIqrXDJuRNoEj6PcaikBT1aUU5tPSo+IRDxUcJoOf6KG29Ps9nasFoPZ/Eh936nZ6IZw75s9384mO32Sa01lq35TEWhoQrpxVmc0I9n7Co+CiKEqB9LJ5PTnLCbr2D7oiej2+T6dBhcg6bUDKB5W7aLbuh9628KIcT6vmERcVHCaB+j9Le6yQvy05uHAfJ+cnPcgxbX0kU/a7IYTf/kQoeK9XabhOm5mcHxHW8ae8NL9rTCnM04SACKj6KogRo63NSWhC/1wNQkOMTgESH3nyeT4Swm91X281f1drhF5/eiRGCNkvUQu/dtMIcOvtdSTlsL9VR8VEUJUBHn2tU6z0Aedl+8UnsB2xftISDkKMefJ5PViA8Nt4EPJ+CkWE3gFb1fkag4qME0HwDpa3XOar1HhjaQ5No8YkVdgMCqdUOm42pBdn0uzwTUmW6vc9Ftt1GQch8KizxOd41MO5jpjsqPoqiBGjvG734+D2f3gSX2Okd9JAXIewWEB9LaOx2CbyviUg6aO91UpKfhcjwPUezSvIAONLRP+5jpjsqPkoQ6vpMdtp7nUzNH13Yzb+rP9Hi0zXgYkpe5FRrGPLGHFbYDXze3XjT2jsYdq2scqolPu0qPqGo+CiKAvgKdXYNuCktyBnVdVOsaghdA4kTnwGXB6fbGxg7lLys4d6Y3SaB9Zfm7vEPgR3tGGBmce6I9im5WRTlOtTzCYOKjxJA13wmN/4yMNOnjE58knF0gH8s/9ih+EOBPX7xEWGmFQL7sHPs4tPZ7wq7ZvRh10DAfiiVJXnq+YRBxUdRFACaA+Iz8ht8NJIhPl2xxMfyfLoH3NgEbDahoigHETg6RvFpaO5h1d0v8okfvoHbM1Qde8Dloa3Xyawwng9A1dQ89XzCoOKjBFDHZ3LTbGVk+cNT8eIvKZNKnk9ukPg4bL6PuSy7jfLCHD7sHJsQ/G7rEQbdXnYf6+KZHccC7UctYZlZrJ7PaFDxURQFgONdY/N87DahKNcR8EYSwWjCbvagqtczS/I4NkbPZ/uRThbPmkJNWT4PvXEg0L6vpReA2vKCsNdVTc2ne9A9YRtc0xUVH0U9HgXw7UWx24SyUVY4AJ8IJEN8psQIu/UMuHEEi8+U3ICnMlqa2vuYU5bPjWfN4Z1DHew80gnAB8e7AairKAx73SkVPlHa19IzpnEzFRUfJYAmHExumrsHKS/MCXs+TiwmsmJ0OGJ6PllBno996P3MKcvncFv/sDWbeDDGcKS9n+qp+dxwejW5WTYefvMgALuPdlFZkhfxDKRTyn2ipOIzHBUfRVEAn+cz2kw3P1NyEys+7b1ORCKfuJqb7fto6xkc7vmcUlGI0+Pl8CjXYFp6Bhl0e6mamkdxfhbXrqjkqXeP0N7r5M3GVs6sLY14bdXUfLIdtkB4TvGh4qME0PN8JjfHOgdGvd7jJ9GeT0vPIGUFOTjs4T/C8oKqcgev+QS8kObReSGH23xiVTU1H4DPfWQOAy4vn//5Rtp6nVy4sCLitXabMHdaAQ2jHDPTUfFRFAVjDE3tfVSX5o/p+oSLT/dg1Ky84CMhchxDj+dZ6zINowyBNbX3Ab60aYDFs4q5evks3m3qZMGMIi5bPCPq9adUFGrYLYS4xEdEVovIHhFpEJE7wryeIyKPWq+/LSI1Qa/dabXvEZHLYtkUkVrLxl7LZna0MURklYhss37eFZGPj/VmTHbU75m8tPQMMuDyUj01fLpwLMoKs2nrdeJN0Jk+scQny24jy1rr8ZfaAZ9ITp+Sw/vHukY1XpMVpqsMuj/f/9QKnvjyR3j8y2cHaslFoq6ikENtfQkvQZTKxBQfEbED9wGXA4uAT4vIopBuXwLajTHzgHuBe6xrFwFrgcXAauCHImKPYfMe4F5jTB3QbtmOOAawE6g3xqywxvixiIQPBCuKEhZ/WGmsnk95UQ5ur6EjQd7PiR4n5YXR16f83k+OY/jH3PKqErYd7hjVeE3t/ZQVZA87P8huE06fUxqobReNpZXFGAO7Ryl6mUw8ns8qoMEY02iMcQLrgDUhfdYAD1mPHwcuEl951zXAOmPMoDFmP9Bg2Qtr07rmQssGls1ro41hjOkzxvi/TuSiX+DHjC75TF78YaWTER8YKtEzkRhjYno+MLTukx0iPivnTOVAa9+oCow2tfcFQm5jYWllMQDbmzrHbCPTiEd8KoHDQc+brLawfSwh6ATKolwbqb0M6AgSk+CxIo2BiJwpIruAHcBfB10fQERuFpHNIrK5paUljretKJMHf1hprB+wfi/kRM/Ei09brxOnxxs4KycSfo8kOyQkdlp1CQDvHGyPe8ym9v5AssFYqJiSy/QpOYG9QZHoHnBxxQ9e5xd/3j/msdKFeMQnXNJ/6HfkSH3Gqz3qPIwxbxtjFgNnAHeKyIiUHWPMA8aYemNMfXl5eRhTijJ5OdzWx7TC7IjHUscikZ6PP016dgwvzb8BNSdr+AFvy6tLyMuy8/re+L6Eer2+PT5VpWP3fMDn/Wxvih7ue33vCXYf6+Kbv9+d8dmn8YhPE1Ad9LwKOBqpj7XeUgy0Rbk2UvsJoCRozSZ4rEhjBDDGvAf0AkvieF9KCEYjlpOWfS09zJ0Wfod+PCRSfA61xRciLLHO7wn1fHKz7Jwzbxovvtcc1wd8S88gTo/3pDwfgKWVJTSe6KVrIPK6WHD1haYMrwcXj/hsAuqsLLRsfAkE60P6rAdush5fD7xsfP+r64G1VqZaLVAHbIxk07rmFcsGls2noo1h2XAAiMgc4FTgQNx3QFEmOcYYPjjew7zpYxefwhwHuVk2WhIQdjscEJ/onkhJwPMZ+TF30cIKjnT0896x7pjjhaZZj5UzaqdiDGw+0BaxT3DduT0fxp5bOhNTfKz1k1uB54D3gMeMMbtE5C4Rucbq9iBQJiINwG3AHda1u4DHgN3As8AtxhhPJJuWrduB2yxbZZbtiGMA5wLvisg24LfA3xhjToztdkxy1PGZlLT0DNLZ74pYmyweRISKotyTOisnXuINEZZYx2YXhclGu2zxDLLswuNbmuIYz8oEPEnxWTl7Ktl2G281RhOf/sBa1p7jmS0+cQV4jTEbgA0hbd8IejwA3BDh2ruBu+OxabU34suGC20PO4Yx5mHg4ZhvQlGUsDQc921+rKsoOik71aV5HLa8hImkobmH2mnhK0gH46/7VhBGfEoLsrl08Qye3NrEP64+ddim1FD8Yb6TDbvlZtlZXl3M242tEfsc6xygbnohWXZboGBppqIVDpQA6vhMTvZaZV/qTiLsBr4EAH9IbKLweg3vf9jNwplTYvYtsuq+9TlHnjwK8NkzZ9PR5+I3mw+Hfd3PobY+ZkzJjSpQ8XLW3DJ2Hu0KnLAayrGOAWYW5zF/eiF7j2d2RQQVH0X390xy3v+wmym5jpipy7GomprPiR7nhO7iP9LRT8+gmwUzYovPOfOmAXDW3PBFPz8yt4z6OVO575V9YY/G9nOorS9mZl28nDW3DI/XsHH/SO/H7fHS3D3ArOJc6qYXsa+lB0+CKkYkAxUfJYCK0ORke1MHy6pK8O3xHjv+D+iJDL3tOuqrELBgZuwQ4cKZU3jn65ewZkXotkQfIsJtl87nw64BfvxaY0Q7B070jnnzbSinz5lKXpadV/eMTPNu7h7Ea2BGcR7zKgoZdHsDyQ6ZiIqPoinWk5gBl4c9H3ZkbJPSAAAgAElEQVSzrKr4pG35P6APtk7cB+amA21kO2wsiiPsBr61nWicfco0rlk+i/teaaCheeQay4meQZq7B1kYh9jFgy/Nu4yX3x+Z5u3PdJtZkhtI/sjk0JuKjxJARWjysftYF26vYVlVyUnbmmsdIz2RRwe8vb+V06pLxmX9xc/Xr1pEQY6dW3+9lT7n8JDhbsvTWjQrPrGLhwsWVNDU3j+iyvWxTl9W3czi3ED17b0ZfAyDio+i4bZJzHarwOby6pP3fKbkZlE1NW/Cimd29DnZfbSLM+eWjavd8qIcfrD2NPYc7+b2J3YMq8y9cX8bdpuwpPLk74+f80/1nf3zyvvDQ2/+NPWZxb5TUWcW57I3jDeWKaj4KCo+k5hNB9uZWZzLjDEeIhfKwplTRn1cQbw8v/s4XgMXRzm4baycN7+cf7xsAb9/9yjfWL8Tr9dgjOGF3cdZObuEKRGOyB4LlSV5nDq9iJfePz6s/WjHAPnZ9sDprPMqMjvjTY8eUAKoCE0uvF7DW/ta+dip5SedbOBn4cwpvPTecfqdHvKyxy80BrBhxzEqS/ICFaLHm7/+2Fw6+13c/9o+jrT3s6J6KnuOd/Nvn1g67mNdung6973SMKw697HOfmYW5wb+L+oqinhk4yG8XoPNNj7/P6mEej6KMkn5oLmb1l4nZ58ybdxsLq8qxmsY9Xk5sTjc1sdrH7Rw7Wmzxk0oQxERbl99Kt+6ZjFv7Gvl3hc/4CNzy7ju9KpxH+vKZTPxGnh214eBtiMd/VQGbWStm15Iv8vDkY7MrPGmno8SyLpRx2dy8eY+316Tj5wyfmso9TWl2ATebGwdV7s///MB7CJ87qyacbMZDhHhprNruPa0Sg639bFw5hTsE+B1nDq9iFPKC3hm+1E+d9YcAI6097M4KLHBn/HW0NwzbqneqYR6PoqKziTl5febmTutgMqSk6tZFkxxXhaLZxXz1r7IJWRGy6HWPn751kGuPa2SGcXjszYVi+K8LJZUFk+I8IBP5K5cNouN+9to7h6g3+mhtdc57P9iKOMtM5MOVHyUAJl+fogyRGe/izf3tXLJ4unjbvujddPYcqid1nGocO31Gr721E7sNuH/XnbqOMwudbjaCr2t33aUA629wPBjIkrysykvysnYpAMVH0UTDSYhr+5pxu01XLZ4xrjbvmrZLDxew4adH8buHIMfvbaPP37Qwj9duZDp45SRlyrUTS/i9DlT+eVbBwMnnC4O2U9UV1GYsXt9VHyUAKpBk4cNO45RUZTDinHYXBrKwplFzJ9eyKObDp2UN/3IxkP8+3N7uHr5LG48c/Y4zjB1uOnsGg609nHHkzuYkuugNuRAv7qKQvYe7x629yhTUPFRtLLBJKO1Z5CX3mtmzYpZE5LCKyJ86dxadh7p4rUP4juqOhiv1/Dd5/Zw55M7+Nj8cr57w7IJy3BLNlcsmcHyqmI8XsOaFZUj1pgWzyqm1+lhvxWWyyQ0203RsNsk47dbj+D2Gm6or47deYx8/LQq/vuVBr71+92cWVsW956fhuZu7nhiB5sPtvOp+mr+5dolZDsy9zuyw27joS+u4o19rVy4YOTm2WVW5YntTR2cUn5yR16kGpn7v6qMHhWhjMcYw6ObDrOiuoT508enWGY4sh027vnEMvaf6OXvHt2K0+2N2v9Qax93PLGd1d9/nb3NPXz3huX823VLM1p4/JTkZ3PF0plh69XNKy8kN8vG9qbOJMxsYlHPR1HNmUS8+kFL4MN9ojl73jT++epFfOv3u/nEj/7MbZfM5+xTppGbZcft8XKgtY+3Glt5Zvsx3mxsJdth47NnzubWC+sCu/4nOw67jSWzilV8lMxG134ynx+/to+Zxblcs3xWQsb7wjm1zCzO5Z/X7+KLv9iMTaAg20Gv041/DX1OWT63XTKfT51RnXEZbePBsqoSfr3xIG6PF4c9czxBFR9F13wmCVsOtvFWYxtfu3JhQsNZq5fM5MIF03l9bwvbmzrpGnBRlOOgqjSfVTWlzCnLz9iEgvFgeXUxP/uzlw+O94zr0Q7JRsVHwR94UxHKXIwx3P3Me1QU5fCZJKQtZztsXLRwOhctHP9NrZmO/6yld5s6Mkp8MseHU8aMik7m8+zOD3nnUAe3XTKf/Gz9zplO1JTlU1qQzZaD7cmeyrii4qMEUBHKTLoHXPzL07s5dXrRhKZXKxODiHD6nKlsPtCW7KmMKyo+iqYZZDj3PPs+x7oG+Nfrlk5YoUxlYjmjZioHWvto7h5I9lTGDRUfRQuKZjB/2nuCX751iC+eU8vK2VOTPR1ljNTXlAKw5UDmhN5UfJQAKkGZxYedA/yfdVupqyjk7y+dn+zpKCfBklnF5DhsbM6gdR8VH0VFJwNxur387SPv0O/y8KMbV2qSQZqT7bCxoroko9Z94hIfEVktIntEpEFE7gjzeo6IPGq9/raI1AS9dqfVvkdELotlU0RqLRt7LZvZ0cYQkUtEZIuI7LD+vXCsN2Oyo+G3zMAYw+1PbGfTgXb+7bplzKuYuDI6SuI4o6aUnUe76HO6kz2VcSGm+IiIHbgPuBxYBHxaRBaFdPsS0G6MmQfcC9xjXbsIWAssBlYDPxQRewyb9wD3GmPqgHbLdsQxgBPA1caYpcBNwMOjuwWKak5m8Z3n9vDbrUf4h0vnJ6ySgTLxnF4zFY/XsO1QR7KnMi7E4/msAhqMMY3GGCewDlgT0mcN8JD1+HHgIvFtWV4DrDPGDBpj9gMNlr2wNq1rLrRsYNm8NtoYxpitxpijVvsuIFdEtDDUGFANSn/+66W9/OjVfXzmzNnccsG8ZE9HGUdWzp6KCGzMkNBbPOJTCRwOet5ktYXtY4xxA51AWZRrI7WXAR2WjdCxIo0RzHXAVmPMiPN7ReRmEdksIptbWkZ/xkgmo+G29McYw/ee38N/vPABnzitkruuWawlazKM4rwsTp1exOYMyXiLR3zC/QaHflpF6jNe7THnISKL8YXi/ipMP4wxDxhj6o0x9eXl5eG6THpUg9ITr9dXOuc/X27gk/VV/PsNyzOqAKUyxKraUt451I7bE/2IinQgnt/QJiB4W3QVcDRSHxFxAMVAW5RrI7WfAEosG6FjRRoDEakCfgv8hTFmXxzvSQlCNSd96Xd6+PKvtvDTP+3n82fX8G+fWKYbSTOYM2pK6XN62H2sK9lTOWniEZ9NQJ2VhZaNL4FgfUif9fgW+wGuB142vljOemCtlalWC9QBGyPZtK55xbKBZfOpaGOISAnwDHCnMebPo3nzio8hj0dlKJ043jXApx54k+d3H+efr17EN69ZPCHHYiupwxnWZtON+9N/3Sem+FjrK7cCzwHvAY8ZY3aJyF0ico3V7UGgTEQagNuAO6xrdwGPAbuBZ4FbjDGeSDYtW7cDt1m2yizbEcew7MwDvi4i26yfkefRKkoG8ae9J7jiB6/T0NzDTz5XzxfOqU32lJQEMKM4l+rSPDZlQNJBXDvPjDEbgA0hbd8IejwA3BDh2ruBu+OxabU34suGC20PO4Yx5tvAt2O+CSUimnCQPni8hv96eS8/eGkv88oL+dGNK3UfzyTjjJpSXtvTgjEmrZNKdFVSCaAalNoc6ejncw++zfdf3MvHV1Ty1K3nqPBMQlbVlNLa66TxRG+yp3JSaM0NRVd6UhxjDL/Z3MRdT+/GGMM91y3lk/XVaf2tVxk7/iKjm/a3cUp5YZJnM3ZUfJQAKkKpR3PXAHc+uYOX3m/mrLml/Pv1y6kuzU/2tJQkckp5AWUF2Ww80MbaVYk/lXa8UPFRNNyWgni8hl+9fZB/f24PTreXb1y1iM+fXaPZbAoiQn3N1LTfbKrio+C11EdFKDXY0dTJV3+3g+1NnZw7bxp3rVnM3DQOryjjzxk1pTy36zjHuwaYPiU32dMZEyo+SkB8lOTSNeDie89/wP+8eYDSghx+sHYF1yyfpWs7ygiC9/tcnabFY1V8FLyW9hhd9UkKbo+XRzYd5t4XPqC9z8nnzprD3196KsV5WcmempKiLJ41hfxsO5sOqPgoaYzHq6KTLF7d08zdz7zH3uYeVtWW8o2rFrGksjjZ01JSHIfdxsrZU9mUxus+Kj6KkgQ+ON7N3c+8x2sftDCnLJ/7bzydyxZP1xCbEjf1NVP5wUt76ex3paWXrOKjBNCln4mnqb2PH7y4lyfeaaIgx8HXrlzI5z4yhxyHPdlTU9KMVTWlGAPvHGznggXpV1FMxUdREkBz9wD3vdzArzceQhBuOruGv72wjtKC7GRPTUlTTps9FYdN2HigTcVHSW/U8Rl/Ovqc3P9aI794Yz8uj+GT9VX87YV1zCrJS/bUlDQnL9vOkspiNqdpkVEVH0WZALoGXPzizwf4yR8b6XG6uWb5LL5y8XxqphUke2pKBnFGzVQeeuMgAy4PuVnpFbpV8VECaHXrk6e918nP/ryfX7xxgO4BN5csms7fXzqfBTOmJHtqSgZyRk0pP3l9P9ubOllVW5rs6YwKFR9FGQdaugf56euNPPzWQfqcHi5fMoNbLpinadPKhBIoMnqgTcVHUSYTH3YOcP9r+3hk4yFcHi9XL5/FLRfMY/50PepAmXhKC7KZV1GYlofLqfgoyhjYf6KXn7zeyOObm/Aaw8dPq+TL55+iNdiUhHNGTSlPbz+K12vSqvCsio+ijIItB9t54I/7eH73cbJsNq6vr+LLHztFjzlQksZps0t4ZOMhGk/0pNXhgio+SgDNNwiP12t44b3j/OSPjWw+2E5xXha3nD+Pvzh7DhVF6VlRWMkcTqsuAWDb4U4VH0XJBAZcHp585wg/fb2RxhO9VE3N45tXL+KG+moKcvRPR0kNTikvpDDHwbbD7Vx/elWypxM3+hekBNCq1j5aewb51duH+J83D3Cix8nSymL+69OncfmSGTjstmRPT1GGYbMJy6qK2Xa4I9lTGRUqPhnCWDeZ6d6eIXYe6eQXbxxg/btHcbq9nH9qOTefN5ePzC3Tgp9KSrOiuoQH/tiYVptNVXwygA07jvE3v3qH5/7uPE6dMbqYb/BxCpNRh9weL8/tOs4v3tjPpgPt5Gfb+VR9NTedPSet4ufK5GZFdQlur2Hnkc7A3p9UR8UnA3jxveMA7DjSOWrxcXq8EzGllKet18kjGw/xy7cOcqxzgOrSPL525UJuqK9Oy/L0yuRmxWx/0kGHio+SHrjck8vz2X20i4feOMDvth1h0O3l3HnT+Jc1S7hgQQX2NNojoSjBVBTlUlmSl1brPio+GYDL41ONLPvoPzwHPZ7xnk7KMeDy8Iedx/jVW4fYfLCdvCw7159exU1n12glAiVjWF6dXkkHKj4ZgNPtE5Acx+gzsfzClYk0tvTwyMZD/GZLEx19LmqnFfDVKxbyyfpqivM1tKZkFiuqS9iw40NO9AwyrTAn2dOJSVyfViKyWkT2iEiDiNwR5vUcEXnUev1tEakJeu1Oq32PiFwWy6aI1Fo29lo2s6ONISJlIvKKiPSIyH+P9UakM4Nu37qNwzZ68XG6h9Z8MkGGXB4vG3Yc47M/fYsL/+M1fv7nA5xzyjR+/b/O5OW//xh/ed5cFR4lI1lRPRWAd9PE+4np+YiIHbgPuARoAjaJyHpjzO6gbl8C2o0x80RkLXAP8CkRWQSsBRYDs4AXRWS+dU0km/cA9xpj1onI/ZbtH0UaAxgAvg4ssX4mHS4racA7hkUbV4YkHDS197Fu42Ee3XyYlu5BKkvy+L+XncoN9VVahUCZFCytLMZuE7Yd7uCihdOTPZ2YxBN2WwU0GGMaAURkHbAGCBafNcA3rcePA/8tvo0Ra4B1xphBYL+INFj2CGdTRN4DLgQ+Y/V5yLL7o0hjGGN6gT+JyLxRvO+Mwu+9jEV8hnk+aZZx4PZ4eXVPC7/eeIhX9jQjwIULKvjsmXM4b365JhAok4q8bDt1FYVsb+pM9lTiIh7xqQQOBz1vAs6M1McY4xaRTqDMan8r5NpK63E4m2VAhzHGHaZ/pDFOxPEeMhq/gIzFiRl0p1/Cwf4TvTy66TBPvNNES/cg5UU53HrBPNaumk2lHk+tTGKWVBbzyvvNGGNSfmN0POIT7h2EfkWO1CdSe7jFiWj9451HRETkZuBmgNmzZ8d7WVoweBKeT/eAO/A4lf2ePqebDTs+5LFNh9l4oA27Tbjg1HI+WV/NBQsqyNKyN4rC0spiHt/SxIddA8wsTu0vYvGITxNQHfS8CjgaoU+TiDiAYqAtxrXh2k8AJSLisLyf4P6RxogLY8wDwAMA9fX1qfw5O2pOJuzWO5i6no8xhnebOnl002F+/+5Regbd1E4r4B9Xn8p1K6uYPkXXchQlGP/JuTuaOjNCfDYBdSJSCxzBl0DwmZA+64GbgDeB64GXjTFGRNYDvxaR7+FLOKgDNuLzYkbYtK55xbKxzrL5VLQxxva2M4s+p09AgkvlxEvv4JDnkyquT1uvk99uPcJjmw6z53g3uVk2rlg6k0/VV7OqtjTlwwmKkiwWzZyCTWDn0S4uXTwj2dOJSkzxsdZXbgWeA+zAz4wxu0TkLmCzMWY98CDwsJVQ0IZPTLD6PYYvOcEN3GKM8QCEs2kNeTuwTkS+DWy1bBNpDMvWAWAKkC0i1wKXhmTjZTS9Tp+AjEV8eoLFJ4l4vIbX97bwm81NPL/7Q1wew/LqEu7++BKuXj6LKbmaHq0oscjLtjOvopCdR1I/6SCuTabGmA3AhpC2bwQ9HgBuiHDt3cDd8di02hsZyogLbo82Rk3UN5Dh+L2XsfiBw9d8Eu/6HGrt4/Eth3l8SxNHOweYmp/F586q4ZNnVLFgxpSEz0dR0p0ls4p5vSH187C0wkEG4Hd4PGNQn5aegXGeTWz8yQO/2XyYt/e3IQIfrSvnq1cu4uJFFeQ40qMkvKKkIksqi3ly6xGauwaoSOF1URWfNCc41OYeQ9ituWuQohwH3RMcfjPGsPlgO7/ZfJhnth+j1+mhpiyff7h0Pp9YWcUsTZFWlHFhaZWVdHCkk4tUfJSJorl7yHPxjkV8ugepmJJDd4t7QqpaH+vs58l3jvD4lib2n+glP9vOlUtnckN9NWfUTNXkAUUZZxbNnIII7DzSldKVDlR80pyjHUPiM1rPxxjDgdZeVlSXsK+ld9zmNODy8MLu4/xmSxN/2tuC18CZtaXccsE8Ll8yg4Ic/bVTlImiIMfB3GkF7EjxpAP9FEhzDpwYEg2Pd3QlDlp6Bunoc3HqjCJe3dPCGBynYexo6uSxzYd5atsRugbcVJbkcesF87ju9CrmlBWcnHFFUeJmSWUxbzfGvQ0yKaj4pDk7jnQi4st0G63n8+5h3zejpdbGtLEkLPQMunlq2xEe2XiInUe6yHHYWL1kBjecXs3Zp5Rh0/pqipJwllYW89S2o4HyU6mIik8aY4zhtQ9aWFVTytv72/CM8mye1z5oJj/bzirr2F33KIrD7Wjq5NcbD7F+2xF6nR4WzCjirjWLWbOiUo+hVpQks3iW7wvlzqOdXHBqRZJnEx4VnzTm9b0n2H+il5vPm8vb+9tG5fn0Oz2s33aUCxZUkJftS22OtUnV4zU8t+tDfvJ6I1sPdZCbZeOqZbP4zJmzOa26RJMHFCVFWFzp2yO3s0nFRxlnOvqcfO13O6mamscnVlby1d/uGFWFg/tf20fXgJsvnF0TOIQu0qmmg24Pj206zE9e38+htj5ml+bzjasWcd3pVerlKEoKMiU3i9oUTzpQ8UlD2nudfP4Xm/iwc4B1f3UWOQ47Dpstbs/n9b0t3PdKA9csn0V9TWnQkQzDw27GGH637Qj/8fwHNLX3s3J2CXdevoBLF8/Qs3IUJcVZUlnMOwfbkz2NiKj4pBmNLT381cNbONjWx39/5jRWzvYdnWu3SVzZbht2HOO2x7Yxr6KQuz/uO/jVYQlJsHgdbuvj9ie288a+VhbPmsL/+/hSPlo3TUNripImLK2cwu/fPUpbr5PSguxkT2cEKj5pxO/fPcqdT+7AYRd+8YUzOPuUaYHXHDaJ6vm0dA/ynWff5zdbmlg5u4QH/qKeIqtYp80m2ATcVtht84E2/vJ/NuPyGO7++BI+fcZszVpTlDRjyayhSgcfm1+e5NmMRMUnDTja0c831+/i+d3HWTm7hP/6zMoRJ3ba7RK2woHT7eWhNw7wny/tZcDt4a/Om8tXLplPbtbw+mn+sN17x7q46WcbqZiSy88+fwa103R/jqKkI4utLRQ7VXyU0dI14OKnr+/nwdcb8RjD7asX8L8+Whv21M5wns+re5q56+ndNLb0cuGCCr525ULmlheGHcthFwZcHv72ka0U5jr49V+emfKHUSmKEpnivCzmlOWzoyk1kw5UfFKQ9l4nv3r7ID95fT+d/S4uXzKDf7piIdWl+RGv8a35+MTnUGsfdz29ixffa6Z2WgE///wZXLAgerql3Sb8euMhnG4vP/2LehUeRckAllQWs+1QR7KnERYVnxTivWNdPPTGAX679QiDbi8XLajgK5fMDxyNGw1/2Ox3W49wx5PbsYtw5+UL+MI5tWQ7RnpKI68Xut1eFswo4qKFqbkvQFGU0bG0sphnth+jvdfJ1BRLOlDxSTIt3YM8s/0o6989yjvWxs1PrKziprPnjOowNbtNeHNfK0+808SqmlK+v3bFqLyX9j4XAJ+sr9aMNkXJEPyls3Yc6eS8FFv3UfFJAl0DLp7d+SG/f/cof244gdfAghlF/NMVC/hkfTUl+aP/huKwCYfa+qgsyePnXziD/OzR/dfOLS+gsaWXNStmjXpsRVFSk+CMNxWfScqAy8NL7zWz/t0jvLKnBafby+zSfP7m/Hlcs2IW86cXnZT9lp5BAD69qnrUwgNwz3XLaOkepKwwNYsQKooyeorzs5hdms/OFKx0oOIzgbg8Xv7UcILfbzvKc7s+pNfpobwoh8+eOZtrls9ixTjWQ/NXHLhy2dg8lzOs4qKKomQWSyuLebcp9ZIOVHzGGWMMWw938LutR3h6+zHaep1MyXVw1bJZrFkxizPnlk1IaZqvX7mILYfadV+OoijDWFJZzDM7Ui/pQMVnnGjtGeSRjYf4zZYmDrb2keOwcfGi6Vy7opLz5k8jx2GPbeQkuO70Kq47vWpCx1AUJf3wJx3sPNrJR+tSZ91Hxeckae4a4N4X9/LEO0043V7OPqWMWy+Yx+olMwLlaxRFUZLFEut4hR1HVHwyAmMMv3zrIP/6h/dxebzcUF/NF8+pYV7FySUOKIqijCcl+dlUl+ax60hXsqcyDBWfMeD1Gv7ptztYt+kwH5tfzl1rFjOnTNdaFEVJTZZWFqfc2T6xt74rI/j+ix+wbtNhbrngFH7++TNUeBRFSWmWVBZzqK2Ptl5nsqcSQMVnlHxwvJv7Xt3HJ1ZW8g+XnqpHDSiKkvKcWevbSvF2Y2uSZzKEis8oefD1/WTbbXz9ykVahkZRlLRgWVUJ+dl23kw38RGR1SKyR0QaROSOMK/niMij1utvi0hN0Gt3Wu17ROSyWDZFpNaysdeymT3WMcYbYwwvvX+cSxZNT6l8eUVRlGhk2W3U15Tyxr40Eh8RsQP3AZcDi4BPi8iikG5fAtqNMfOAe4F7rGsXAWuBxcBq4IciYo9h8x7gXmNMHdBu2R71GKO9EfFwvGuQEz1O6mumToR5RVGUCePceWU0NPew/0RvsqcCxOf5rAIajDGNxhgnsA5YE9JnDfCQ9fhx4CLxxaTWAOuMMYPGmP1Ag2UvrE3rmgstG1g2rx3jGONOa6+vflpFkdY/UxQlvbh2RSV2m/Afz+9J9lSA+FKtK4HDQc+bgDMj9THGuEWkEyiz2t8KubbSehzOZhnQYYxxh+k/ljECiMjNwM3W0x4RaQVORHzXUbj8nrFcldJMY4z3IgPRe+FD78MQGXUv7gPu++yYLp0GzBmvecQjPuFW1U2cfSK1h/O4ovUfyxjDG4x5AHjA/1xENhtj6sNcO+nQezGE3gsfeh+G0Hvhw7oPNeNlL56wWxNQHfS8CjgaqY+IOIBioC3KtZHaTwAllo3QsUY7hqIoipKixCM+m4A6KwstG9/i/vqQPuuBm6zH1wMvG2OM1b7WylSrBeqAjZFsWte8YtnAsvnUGMdQFEVRUpSYYTdrfeVW4DnADvzMGLNLRO4CNhtj1gMPAg+LSAM+b2Stde0uEXkM2A24gVuMMR6AcDatIW8H1onIt4Gtlm3GMkYMHojdZdKg92IIvRc+9D4MoffCx7jeB/E5D4qiKIqSOLTCgaIoipJwVHwURVGUhDMpxSdWuaBMQER+JiLNIrIzqK1URF6wShe9ICJTrXYRkf+07sd2EVkZdM1NVv+9InJTuLFSGRGpFpFXROQ9EdklIv/Hap9U90JEckVko4i8a92Hb1ntKVvOaqKxqq1sFZGnreeT8l6IyAER2SEi20Rks9U28X8fxphJ9YMvwWEfMBfIBt4FFiV7XhPwPs8DVgI7g9q+A9xhPb4DuMd6fAXwB3x7ps4C3rbaS4FG69+p1uOpyX5vo7wPM4GV1uMi4AN8JZ0m1b2w3k+h9TgLeNt6f48Ba632+4EvW4//BrjferwWeNR6vMj6m8kBaq2/JXuy398Y78ltwK+Bp63nk/JeAAeAaSFtE/73MRk9n3jKBaU9xpg/4ssKDCa4RFFo6aL/MT7ewrfXaiZwGfCCMabNGNMOvICvfl7aYIw5Zox5x3rcDbyHrwLGpLoX1vvpsZ5mWT+GFC5nNZGISBVwJfBT63lKl/ZKAhP+9zEZxSdcuaAR5XgylOnGmGPg+1AGKqz2SPcko+6VFS45Dd+3/kl3L6ww0zagGd+Hwz7iLGcFBJezSuv7YPF94B8Br/U87tJeZN69MMDzIrJFfGXIIAF/H5PxGO24yvFMMk6qdFE6ICKFwBPA3xljuiTyWUwZey+Mb//bChEpAX4LLAzXzfo3Y95nRXgAAAGtSURBVO+DiFwFNBtjtojI+f7mMF0z/l5YnGOMOSoiFcALIvJ+lL7jdi8mo+czmcvxHLdcZKx/m6320ZZBSitEJAuf8PzKGPOk1Twp7wWAMaYDeBVfzH4ylrM6B7hGRA7gC7tfiM8Tmoz3AmPMUevfZnxfSlaRgL+PySg+8ZQLylSCSxSFli76CyuT5Syg03K1nwMuFZGpVrbLpVZb2mDF5h8E3jPGfC/opUl1L0Sk3PJ4EJE84GJ861+TrpyVMeZOY0yV8RXJXIvvvX2WSXgvRKRARIr8j/H9Xu8kEX8fyc60SMYPvoyND/DFvL+a7PlM0Ht8BDgGuPB9K/kSvjj1S8Be699Sq6/gq7S+D9gB1AfZ+SK+hdQG4AvJfl9juA/n4nP/twPbrJ8rJtu9AJbhK1e13fpw+YbVPhffB2YD8Bsgx2rPtZ43WK/PDbL1Vev+7AEuT/Z7O8n7cj5D2W6T7l5Y7/ld62eX//MwEX8fWl5HURRFSTiTMeymKIqiJBkVH0VRFCXhqPgoiqIoCUfFR1EURUk4Kj6KoihKwlHxURRFURKOio+iKIqScP4/oF+3jiatpawAAAAASUVORK5CYII=\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(770, 785)\n",
"plt.savefig('test.png')\n",
"# print(jpsi_width)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"\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": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"# 0.213/(0.00133+0.213+0.015)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Adjust scaling of different parts"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"# total_f.update_integration_options(draws_per_dim=2000000, mc_sampler=None)\n",
"# inte = total_f.integrate(limits = (1000, 1040), 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": null,
"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",
"# name = \"phi\"\n",
"\n",
"# print(name+\":\", inte_fl)\n",
"# print(\"Increase am by factor:\", np.sqrt(pdg[name+\"_BR\"]/pdg[\"NR_BR\"]*pdg[\"NR_auc\"]/inte_fl))\n",
"# print(\"New amp:\", pdg[name][3]*np.sqrt(pdg[name+\"_BR\"]/pdg[\"NR_BR\"]*pdg[\"NR_auc\"]/inte_fl))\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": null,
"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": null,
"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": null,
"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": null,
"metadata": {},
"outputs": [],
"source": [
"# print(36000*(1+ pdg[\"jpsi_BR\"]/pdg[\"NR_BR\"] + pdg[\"psi2s_BR\"]/pdg[\"NR_BR\"]))"
]
},
{
"cell_type": "code",
"execution_count": null,
"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": null,
"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": null,
"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": null,
"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": null,
"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.05, dtype=dtype),\n",
"# tf.constant(0.93, dtype=dtype),\n",
"# tf.constant(0.05, dtype=dtype),\n",
"# tf.constant(0.065, dtype=dtype),\n",
"# tf.constant(0.05, dtype=dtype)]),\n",
"# components_distribution=tfd.Uniform(low=[tf.constant(x_min, dtype=dtype), \n",
"# tf.constant(3090, dtype=dtype),\n",
"# tf.constant(3681, 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(3102, dtype=dtype), \n",
"# tf.constant(3691, dtype=dtype),\n",
"# tf.constant(3110, dtype=dtype), \n",
"# tf.constant(3710, dtype=dtype)]))\n",
" dtype = tf.float64\n",
" mixture = tfd.MixtureSameFamily(mixture_distribution=tfd.Categorical(probs=[tf.constant(0.04, dtype=dtype),\n",
" tf.constant(0.90, dtype=dtype),\n",
" tf.constant(0.02, dtype=dtype),\n",
" tf.constant(0.07, dtype=dtype),\n",
" tf.constant(0.02, dtype=dtype)]),\n",
" components_distribution=tfd.Uniform(low=[tf.constant(x_min, dtype=dtype), \n",
" tf.constant(3089, dtype=dtype),\n",
" tf.constant(3103, dtype=dtype), \n",
" tf.constant(3681, dtype=dtype),\n",
" tf.constant(3691, dtype=dtype)], \n",
" high=[tf.constant(3089, dtype=dtype),\n",
" tf.constant(3103, dtype=dtype), \n",
" tf.constant(3681, dtype=dtype),\n",
" tf.constant(3691, dtype=dtype), \n",
" tf.constant(x_max, 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": null,
"metadata": {},
"outputs": [],
"source": [
"total_f._sample_and_weights = UniformSampleAndWeights"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"0.00133/(0.00133+0.213+0.015)*(x_max-3750)/(x_max-x_min)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# zfit.settings.set_verbosity(10)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"# zfit.run.numeric_checks = False \n",
"\n",
"nr_of_toys = 1\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": null,
"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": null,
"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_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": null,
"metadata": {},
"outputs": [],
"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": null,
"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": null,
"metadata": {},
"outputs": [],
"source": [
"# jpsi_width"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# plt.hist(sample, weights=1 / prob(sample))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Fitting"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"start = time.time()\n",
"\n",
"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": null,
"metadata": {},
"outputs": [],
"source": [
"print(\"Time taken for fitting: {}\".format(display_time(int(time.time()-start))))\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": 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-6)\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": [
"# _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": null,
"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": null,
"metadata": {},
"outputs": [],
"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
}
