{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Import"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"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": null,
"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": null,
"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": null,
"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": null,
"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",
" '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",
" ] # 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'],\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 = 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": null,
"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": [],
"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), 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))\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))\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))\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": 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",
" 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",
" \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": 15,
"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": 18,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ8AAAD8CAYAAACo9anUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJztvXt4XNV57/95Z0ZXW9bN8lWSJdvCV2yDjTFgCOVOIDFNyIlJIaThhKYN5/Q0p23CSW+hSX8hp79w0pPQhAIJJQ2G0BBcQiAJGCjE2Mj4ho0v8l2WbVl3Wde5rPPH3jMejWY0W7JmtqR5P8+jx3vWXvtda29L+zvvWu96lxhjUBRFUZR04nG7A4qiKErmoeKjKIqipB0VH0VRFCXtqPgoiqIoaUfFR1EURUk7Kj6KoihK2nEkPiJyi4jsF5E6EflqnPM5IvKsfX6LiFRFnXvQLt8vIjcnsyki1baNg7bNbAdtLBORzSKyR0R2i0juSB6GoiiKkh6Sio+IeIHvA7cCi4G7RGRxTLX7gFZjzHzgEeBh+9rFwHpgCXAL8KiIeJPYfBh4xBhTA7Tatodqwwf8BPiiMWYJcC3gH+ZzUBRFUdKIE89nNVBnjDlsjOkHNgDrYuqsA56yj58HrhcRscs3GGP6jDFHgDrbXlyb9jXX2Tawbd6RpI2bgF3GmJ0AxphmY0zQ+SNQFEVR0o3PQZ3ZwImoz/XA5YnqGGMCItIOlNrl78ZcO9s+jmezFGgzxgTi1E/UxkWAEZFXgTIssft27E2IyP3A/QCTJk1auXDhQge3riiZS2evn6PN3cwvm0xetnfQ+ZNtPXT0+Fk0c4oLvTvP3lMdFOZlMbsoD4C6xnP4vEJV6aS0tN/ZG+BocxcAF88uTEubbrFt27YmY0zZaNhyIj4Spyw2J0+iOonK43lcQ9Ufqg0fsBa4DOgGXhORbcaY1wZUNOYx4DGAVatWmdra2jjmFEUJ8+qe0/zR09t47r+vZcmswS/Vr72wm1c+OE3tX9/oQu/Os+KhX7Nu+Sy+vm4pAOu+9zbFk7L58R+uTkv7m/Y38oc/eg+A2m/dlpY23UJEjo2WLSfDbvVARdTncqAhUR17DqYQaBni2kTlTUCRbSO2raHaeNMY02SM6QZeBi51cF+KogxBIGh97/N54r8mvB4hOAZyQ4ZCBmsE3kJECKWxW5ofc2Q4EZ/3gBo7Ci0bK4BgY0ydjcC99vGdwOvG+h/ZCKy3I9WqgRpgayKb9jWbbBvYNl9M0sarwDIRybdF6SPAXuePQFGUeARCIQB83niDDuARIZjOt3wCjLH6EsYj6RUE+zEpwyTpsJs9v/IA1kveCzxpjNkjIg8BtcaYjcATwNMiUofljay3r90jIs9hiUEA+FI4GCCeTbvJrwAbROQbwHbbNkO00Soi38ESNAO8bIz55QU9FUVRojyf+OLj9QihMSA+IWOI7qJHhFA6xUc9nxHhZM4HY8zLWMNZ0WV/E3XcC3wqwbXfBL7pxKZdfhgrGi62fKg2foIVbj1i/H4/9fX19Pb2XoiZMUNubi7l5eVkZWW53RVlnBL2anze+AMkPo8QGBPiAx5PtOcjafVGxsAjGJc4Ep9MoL6+noKCAqqqqgaMH49HjDE0NzdTX19PdXW1291Rxin+8LBbAs/H40mvh5GIkDFE/8mKpNcb0TmfkaHpdWx6e3spLS0d98ID1oRraWnphPHiFHcIez7eRMNuY3bOR0inHkQ3pULkHBWfKCaC8ISZSPeiuEN4zicrQbSb5fm4/8IdNOfjSa/nE93WWBDj8YKKj6IocQlHu3kTRLt57S84br9vLfGJmfNJq/icPx4LoefjBRWfccYbb7zB7bffDkBfXx833HADK1as4Nlnn3W5Z8pEIxxMkGjOJxyCHXA51jhkGLTOJ+jSOh8Nu3aOBhyMY7Zv347f72fHjh1ud0WZgCQLtQ57G26+cMMv/oGh1mle5xPVliXEg1MRKYNRz2cMcfToURYuXMi9997LsmXLuPPOO+nu7uaVV15h4cKFrF27lp///OcANDY2cvfdd7Njxw5WrFjBoUOHXO69MtEIJAs4sN8ebg41hYe8oofdvOkedgvFP1aGRj2fOHz9P/awt6FjVG0unjWFv/3YkqT19u/fzxNPPMFVV13F5z//eb7zne/wwx/+kNdff5358+fz6U9/GoBp06bx+OOP84//+I+89NJLo9pXRQEIBEP4PJIweCX8wndzkj0Ux/ORtK/ziQo40Dkfx6jnM8aoqKjgqquuAuDuu++mtraW6upqampqEBHuvvtul3uoZArBkEno9cB5j8jNLAfhF7/EpNdJ6zqfqGO357/GE+r5xMGJh5IqYr9ltre3a9i04gqBkCErQXYDOD8X5GaWg/C73tV1PhpwMCLU8xljHD9+nM2bNwPwzDPPcMMNN3DkyJHInM4zzzzjZveUDCIQDA3p+YRT2riZ5SDesFv61/mcP9ZhN+eo+IwxFi1axFNPPcWyZctoaWnhz/7sz3jssce47bbbWLt2LXPmzHG7i0qGYHk+Qwy7jYE5n/DLPlokxcXEosF0xniPc3TYbYzh8Xj4wQ9+MKDslltuYd++fYPqXnvttVx77bVp6pmSaQSCQ8/5hD0fVwMOQmHPx71hN/V8RoZ6PoqixCUQMgk3koPoDAcuej5xwsHTHnCg6XVGhIrPGKKqqooPPvjA7W4oCmBFbiXaSA6iMxy4P+w2aEuFdHo+UY2NhSzf4wUVnyjcTpA4mkyke1HcIZAk1Pp8hgP3o928A9LrpNcDGRBqrXM+jlHxscnNzaW5uXlCvLTD+/nk5ua63RVlHBMMmoQZreH8UJeb8xznAw7Ol1lzPu5Eu6nn4xwNOLApLy+nvr6es2fPut2VUSG8k6mijJRAKEmo9RiIdosXcOBN87CbzvmMDBUfm6ysLN31U1GiSBpqHclwkK4eDSZuwIGL+/mMhW3Fxws67KYoSlyShVqfz3DgnvokXueTvj7osNvIUPFRFCUu/cEQ2b7Er4gxkeEgQai1W1sq6LCbc1R8FEWJS38gNGRut/MZDtLVo8FEtn1wcSfT6KbcjPwbb6j4KIoSF38wRM6Qno/1r6vpdcIBBy6u8zE65zMiVHwURYmLP+jM8xkLiUVj1/loYtGxj4qPoihxSTbsNiYyHMSd80l3brfoLRVUfJyi4qMoSlz8waH38xkTGQ7iptdxz/PRYTfnqPgoihKXZNFu3jGwmVwwTnodj0h60+uo5zMiHImPiNwiIvtFpE5EvhrnfI6IPGuf3yIiVVHnHrTL94vIzclsiki1beOgbTN7qDZEpEpEekRkh/0zcD8CRVFGRH8gRPZQiUXtiIOgm+t8IgEH58vExWE3nfNxTlLxEREv8H3gVmAxcJeILI6pdh/QaoyZDzwCPGxfuxhYDywBbgEeFRFvEpsPA48YY2qAVtt2wjZsDhljVtg/XxzWE1AUJS7JAg7C2Q/8LibTjBdw4HVzJ1P1fBzjxPNZDdQZYw4bY/qBDcC6mDrrgKfs4+eB60VE7PINxpg+Y8wRoM62F9emfc11tg1sm3ckaUNRlBTgTzLs5rOFydUMBwkCDlzbyVTFxzFOxGc2cCLqc71dFreOMSYAtAOlQ1ybqLwUaLNtxLaVqA2AahHZLiJvisjVDu5JUZQhMMYkDTgIp9dx0/OJt59PutPrGPV8RoSTxKLxvIvYJ5yoTqLyeL/RQ9Ufqo1TQKUxpllEVgK/EJElxpiOAR0UuR+4H6CysjKOKUVRwvTbM/lDeT5hYXJzD5ug3bYvJtoNLAFNx+CIbiY3Mpx4PvVARdTncqAhUR0R8QGFQMsQ1yYqbwKKbBuxbcVtwx7SawYwxmwDDgEXxd6EMeYxY8wqY8yqsrIyB7etKJlL2JvJdrTOx/3Eop6YaDcgbd6PhlqPDCfi8x5QY0ehZWMFEGyMqbMRuNc+vhN43VjxhxuB9XakWjVQA2xNZNO+ZpNtA9vmi0O1ISJldgADIjLXbuOw80egKEos/oAlKENtqRDeaM7VgIMEiUUhfV6ILjIdGUmH3YwxARF5AHgV8AJPGmP2iMhDQK0xZiPwBPC0iNRheTzr7Wv3iMhzwF4gAHzJGBMEiGfTbvIrwAYR+Qaw3bZNojaAa4CHRCQABIEvGmNaRv5IFEUJD7tlDRlwYHs+LmYWTbSlArgjPjrn4xxHm8kZY14GXo4p+5uo417gUwmu/SbwTSc27fLDWNFwseVx2zDG/Dvw70lvQlEUx/RHPJ/xkV4n3rBbuqZfQsbadM8fNDrsNgw0w4GiKIPw297MUFmtzw+7uef5hOJ4PmG9TJcQhMz5BbcacOAcFR9FUQYRGXYbKrebR/CIy9FucdLreCOZF9LTL2NMxAt0c2+j8YaKj6Iog/AHrBf3UOID1kJTv4vRbqE46XXCMRLpEp9Q6HxUoJuphsYbKj6KogzCyTofgCyPuOv5xBt286bX8wma84tx1fNxjoqPoiiD6HcQag3WS9/VaLc422iHF5ymzfMxJiJ+mljUOSo+iqIMIhxEMNQiU7CG5fxjbD+fsBCla/GrMZYIez2i63yGgYqPoiiD8DscdvN53fV8AnHS64S9kHRNv4SMwSOW6GmotXNUfBRFGYSTdT5ghRi7+cKN5/mkO+1PyFhrizxp3sphvKPioyjKIHoDQQBys7xD1svyuhtwEBa+eJ5P+qLdDCKWEGuGA+eo+CiKMohev+U15GYlD7V2M7FoeMjPFxVrfX7OJ70BBx7R9DrDQcVHUZRB9PotzyfHN7Tn4/OIq4lFw21HR+Wl3fMxBo9YAQcqPs5R8VEUZRB9AWeeT5bX427AQSiE1yMD9u05n23AEoIn3z7C2odfT1kfQsZKZur1iIZaDwNHiUUVRckswp5PbjLPx+tuhFcgaAYsMIXziUXD/Xropb0A9AWCST25kWDC0W4aaj0s1PNRFGUQvf4Q2V7PgCiyeGR5PK4mFg2EDFkxfUyU5PNcbyAlfQiG7GE3DbUeFio+iqIMotcfHDKjdRif29FuwRC+mHDwsCcU26/u/mBK+hAydpJV9XyGhYqPoiiD6AuEyEkSZg3hxKIuBhyEzKAUQLFzPmG6+lPj+YQXmfp0zmdYqPgoijKIPn8wabABhBOLuhtqHR1mDdFzPgP71dWXGs/HRBaZ6rDbcFDxURRlEL2BYNIFpjAWht3O76UTJrzgNHbOp6svtZ6PV3TYbTio+CiKMohef8iR5+P2fj7WsJuzOZ8ef2o8n2DInA+1VvFxjIqPoiiDsAIOkns+bu/nYw27DfR8YheZhs/3pkh8jLG8HhWf4aHioyjKIPoCzjyfLK+7odb+oBkU7RbZz8cM3I01ddFuBo8HXWQ6TFR8FEUZRK8/mHSBKUBOlieSAdsNAqHQoGi3WM8nfL4nleIjgkfU8xkOKj6Kogyi1+8s4CDb642k4nGDQNAkHHYLDwdm2yKaqjmfcHodn0d0S4VhoOKjKMogev0hR4tMc7I89AVS81J3gn+IRaZhLySsTan1fKyFpm7Of403VHwURRmE00WmOT4P/qBxLcQ4EG+Rqb3uJzz/EhahVM75eO30Our5OEfFR1GUQThdZBreZrvfpaCDeItMI8NutuiE/03ZsFvIHnbzuru9xHhDxUdRlAEYY+jqDzApO3nS+3A4dp/fHfHxDzHnE7QFMez59KQ4vU6Oz93gi/GGio+iKAPoC4QIGcjPcRBwYHs+fUF35n2CocEZDiLiYzsh4VDwVHk+4fQ6OT6vax7geMSR+IjILSKyX0TqROSrcc7niMiz9vktIlIVde5Bu3y/iNyczKaIVNs2Dto2s5O1YZ+vFJFzIvLnw30IiqKcJ5yGxpnnY4uPW55PaHDAQWSdT2ig55OqOZ+gvc4nx+du8MV4I6n4iIgX+D5wK7AYuEtEFsdUuw9oNcbMBx4BHravXQysB5YAtwCPiog3ic2HgUeMMTVAq207YRtRPAL8yumNK4oSn/BLOj/bWcABuDnnM3g/n+g5H2NMZM4nVRkOwut8sn0e10R4POLE81kN1BljDhtj+oENwLqYOuuAp+zj54HrxdrXdh2wwRjTZ4w5AtTZ9uLatK+5zraBbfOOJG0gIncAh4E9zm9dUZR4hLcemJQzDjyfYChhbrdg0AxY9Jn6YTePDrsNAyfiMxs4EfW53i6LW8cYEwDagdIhrk1UXgq02TZi24rbhohMAr4CfH2omxCR+0WkVkRqz549m+SWFSVzCW894MzzsQMOXBpu6g+EyImJyvPK+fQ60VscpGzYLWQHHGR51fMZBk7EJ94+urHxhInqjFb5UG18HWuY7lyc8+crGvOYMWaVMWZVWVnZUFUVJaPpHoHn41aUV18gNCgBqscjeMQShWjPpzeF4uP1eCJzPkbX+jgi+W+X5X1URH0uBxoS1KkXER9QCLQkuTZeeRNQJCI+27uJrp+ojcuBO0Xk20AREBKRXmPM9xzcm6IoMQzH84lEu7kmPsFIH6LxeazFr9EZB7pTuKWCzyNkez2ETPyFr8pgnHg+7wE1dhRaNlYAwcaYOhuBe+3jO4HXjSX/G4H1dqRaNVADbE1k075mk20D2+aLQ7VhjLnaGFNljKkC/g/wDyo8ijJyIp7PcNb5uCA+oZDBHzRx0wBlea0dVsO7meZleVOWXicQMni9Ehn+07U+zkj622WMCYjIA8CrgBd40hizR0QeAmqNMRuBJ4CnRaQOyxtZb1+7R0SeA/YCAeBLxpggQDybdpNfATaIyDeA7bZtErWhKMro0hWOdhvGOh83Xrjhyf14+w5l2ZP/4WG3glwfjZ199hDZ6HolwZC1p1C0EE/KGdUmJiROht0wxrwMvBxT9jdRx73ApxJc+03gm05s2uWHsaLhYssTthFV5++GOq8oSnK6R7LOx4WAg/Dkfrxht/A+Q4EY8en1Bx3NZQ2HgC1obj6L8YhmOFAUZQBhzyfPSWJRF4eawi/5eMNu2V4P/YHzcz4FuVlAaiLewnM+Ouw2PFR8FEUZQHdfgPxsLx4Hw1PZXvcCDsJtxhUfX9jzseoU5FreTioWmgbsaLdsr3vzX+MRFR9FUQbQ1R90FOkGRLZdcGXYLTDUsJvgj5nzgRR7Pi4vuB1vqPgoijKAzl5/ZJgqGW6+cM8Pu8UJOLDnfMJbHBTkWPcz2lkOjDGRIIbwsJvO+ThDxUdRlAF09AaYkutsUt7nsbaPTlXqmqEIz63EZjgAS3z6o9LrnPd8RndbhbB9n0ci24678SzGIyo+iqIMoKPHz5Q8Z56PiJCX7U1Z6pqhiMz5eOMHHPgD0XM+1v2M9pxPOJrO65VIdGA4K7gyNCo+iqIMoKPXufiAlQlhtD0KJ/QN5fn4ZECo9WTb8+npH93hwWjPZ7Idwn2uTz0fJ6j4KIoygI6eAFMczvmAtR7IFc/HP/ScT38wFBVqnZpht4jn4/EwyV6Uq56PM1R8FEUZgOX5OF+ImZedutQ1QxHOcJBokWl/4Hy025QUhVpHez6TIp6Pio8TVHwURYnQ6w/SHwgNy/PJz/ZG9gBKJ91DLIbNDke72XM+k3NSs8g0PKcUznDg84h6Pg5R8VEUJUJHrx9gmHM+Plc8n54hdly1FpkagsGYOZ8Uej4ilvejno8zVHwURYnQ0WO9OJ2GWkM44CD94hP2tvLj5KALLzINz8lkez1MyvbS2TvKcz62uIWzQUxW8XGMio+iKBHaeyzPp3AYno9bodY9/UFEIDfBOp/oDAc+r1CUn01rd/+o9iHa8wGYlOPVYTeHqPgoihJhZMNu7oRad/UFyc/yIjI4B1044CB6TqYoP4u2bv+o9iFs32evNZqU44tsxqcMjYqPoigR2ruH7/m4FWrd4w+Qn2B7hPCcT3hYLMvjoTgFnk8kv5wtPlNysyICrgyNio+iKBGau6yXc+mkbMfX5GV76YsKa04XXX2JE6BmeWXAZnJeb2o8n/6YzNolk7Jp6RpdgZuoqPgoihKh+VwfPo8MO9QaRn8BZzK6+4Nxgw3AGnYLhgx99logn0cozs+mbZQ9HxWfkaPioyhKhJaufoonZTvayydMWADSHW7d3R8YwvOxXm09tiB6PUJxfhbtPX5Co+ihxW7rUDIpm+7+YEr2DZpoqPgoihKh6Vz/sIbcgEhamc40R3l1D7HvUNgTCc9FZXk8FOVnEzKM6pxMfxzxAdT7cYCKj6IoEVq6+iidPDzxCQ/RjfYammQM5fnk2eXn7D55vRIRhuZRFIZwip9wfjkVH+eo+CiKEqG5q5/SSTnDuiYcGdfRk94or3O9gUjanFjCKXfCCz59HmH6lFwAzrT3jlofYj2fUhUfx6j4KIoSoeVcf+Tbu1PCa4La0yw+7T3+hCHhYfEJDwV6PcLMQkt8To2i+IR3LQ2Lz9TJlnA3dvaNWhsTFRUfRVEA60Xa2Rdg6giH3dK5vsUfDNHVH0woPrkxw24+jzDDFp/THaPv+YTnmMJtnGztGbU2JioqPoqiAFawAUDp5JEOu6VvzqcjkgYofqh19LCb1076mZvlpTg/i1PtoycMfTHik5vlZVpBDifbuketjYmKio+iKACctoejwt/enZKb5SHLK2kddovkoMtPMuzW68cbFTY+ozAvcp+jQXhOaVLUeqPZxXmcbFPPJxkqPoqiAOfFZ+YwxUdEKMxLb1qZZAlQo6PdsqLEZ1ZhLvUjHBI71d7D7f/3P/nubw9Gys71BpiU7R2wLmp2UZ4OuzlAxUdRFIDIcNTMKXnDvnZKblZao92SiU+u73zAQbTnM7dsEkeauka00PSZrSf44GQHj/z2QMSzOdcXiOxgGqa8OJ+TbT0E7DBsJT6OxEdEbhGR/SJSJyJfjXM+R0Setc9vEZGqqHMP2uX7ReTmZDZFpNq2cdC2mT1UGyKyWkR22D87ReT3R/owFCWTOd3eS26WZ1hbaIcpyMtyZ9gtYcCB9Wo71xeIZJwGmFc2mb5AaETDYu/UNVFWYM2H/XJXQ8T+5Ji9jy6aPhl/0HC0uWvYbWQSScVHRLzA94FbgcXAXSKyOKbafUCrMWY+8AjwsH3tYmA9sAS4BXhURLxJbD4MPGKMqQFabdsJ2wA+AFYZY1bYbfxQRIb/16MoGc6pjl5mFubF3aIgGYV56fV8Wu11NEX58SPzwnM+xjDA85k3bTIAh86eG3abB890csuSGSydPYWXd58GrIW1k3NixacAgP2nh99GJuHE81kN1BljDhtj+oENwLqYOuuAp+zj54HrxfoNXgdsMMb0GWOOAHW2vbg27Wuus21g27xjqDaMMd3GmHCYTS6Q3tS6ijJBON3ey4wpw5vvCVM6KXtUMwck46ydALUkgfjkZp3PfBA95zN36iQA6hqHJwzt3X46egNUluRz69KZ7DjRRkNbD81dfYPSEc2fNhmvR9h/umNYbWQaTsRnNnAi6nO9XRa3ji0E7UDpENcmKi8F2qLEJLqtRG0gIpeLyB5gN/DFqOsjiMj9IlIrIrVnz551cNuKklmcausZdqRbmKmTs2k614cx6fnu19jRx9TJOQkToGZ5rQg8gJwoISqZlM3UyTnsPTU8YTjRaoVOV5TkcevSGQC88sFpznT0Ma1g4DPLzfJSVZrP3lOdw2oj03AiPvH+d2N/wxLVGa3yIfthjNlijFkCXAY8KCKD/oKMMY8ZY1YZY1aVlZXFMaUomUtfIMipjl4qSvJHdH1ZQQ69/lAk9DjVnD3XF5l/SUTY+wmvwQErMu+SyiJ2HG8bVnvHW8Lik8/cssksnFHAL3acpPlcH9PjCPYllcW8f7w1bWI8HnEiPvVARdTncqAhUR17vqUQaBni2kTlTUBR1JxNdFuJ2ohgjPkQ6AKWOrgvRVFsTrT0YAxUlY5MfMJpZcILVVPN2c7k4hNeexMtPgArKoo43NQ1rL19osUH4KMXz2RXfTshAwvsOZ5oVleV0NLVz6GzGnSQCCfi8x5QY0ehZWMFEGyMqbMRuNc+vhN43ViSvxFYb0eqVQM1wNZENu1rNtk2sG2+OFQbtg0fgIjMARYARx0/AUVROGZHZs0pnTSi68NCcDZNOc3OdvYxLYn4hCPhwhmnw1xSWQTAtmOtjts70dJNcX5WJJXQZy6vZEquj/xsL1fOKx1U/7LqEgDeO9oy6JxikTQqzBgTEJEHgFcBL/CkMWaPiDwE1BpjNgJPAE+LSB2WN7LevnaPiDwH7AUCwJeMMUGAeDbtJr8CbBCRbwDbbdskagNYC3xVRPxACPgTY0zTyB+JomQex5qtb/YX7vmkXnz6AyGazvUxLUlwRER8sgZ+x760spi8LC9v7D/L9YumO2rzeEv3gCHJqZNzeOV/XIMIFMdJxFpVms+MKbm8deAsd62udNRGpuEoJNkY8zLwckzZ30Qd9wKfSnDtN4FvOrFplx/GioaLLY/bhjHmaeDppDehKEpCjjV3UZDrG3ZG6zDp9HxOtvUQMjAnyfxUeL1Stneg+ORmeblq/lRe39fIQ8Y4Ci0/0dLNktmFA8pmFSVejCsiXL9oGi9sP0mvPzgg+k6x0AwHiqJwpLmbOaX5I1rjA1Ccn43XIzR2jl7etESE518qk3hp4a0esn2DX3PXL5rGybYe9jQkj3oLhgwn23qoHGYwxg2Lp9PdH2Tz4WbH1xhjeK72BEebJv5ckYqPoijUnelkftnkEV8f3i9npHnThsNxe34qmRiEh91iF4EC3Lp0Btk+Dz+rPTHoXCyn2nvwBw0VxcMTnyvnlTI5x8cvd51yfM3WIy385fO7+MK/1g6rrfGIio+iZDjt3X4a2ntZOHPKBdmpLMnnREvqtxI41txNjs9DWZKtH8LiE5t7DazMCDcvmcEvdjTQ0x8c0k7Y05ozzPmwHJ+X25fN5OXdpxyHoB+wF78ebDyX1kStbqDioygZzj57Jf7CGYNDhodDRXE+J9Lg+ew/00nN9MkJF5iGCW+bnSjn3GevmEN7j59nth4f0k5YUIc77AbwqVXldPcHedmh99MQlXNur4MhwfGMio+iZDj7Tlsr8RfOuDDPp6Ikj7OdffT6h/YkLgRjDHsbOljkoK9X10xlWkEOd64sj3v+sqoS1swt4QdvHhqyz8cfDZWRAAAe7UlEQVSau/FFbcM9HC6tLGZu2SQ2vDe0wIVpaOuJrEvaN8wsDOMNFR9FyXD2ne6kKD+L6VOGt4NpLOFQ5PrW1A29ne3so7mrn8WzkotPeXE+W792A1fNn5qwzp/dcBGNnX08uqkuYZ3jLd3MLs4bkB3bKSLCPWvm8P7xNrYfT76uqKGthxUVRZRMyubDCZ6eR8VHUTKcnSfaWDJryogj3cKEF6imclX/zvp2AJbMKkxS0xmXzy1l3YpZ/PObhxImG61rPBdJSDoSPrWqgoJcH0+8fSRp3Ya2XmYX57FwRkFkOHSiouKjKBlMV1+Afac7WFlZfMG2Lpo+GRHYl8Jv7JsPNZPj87C8YnTEB+Brty1iUo6PP92wfdDwW68/SF3juQsSu8k5Pu5aXcmvPjg95D5CgWCI0x29zC7KY8GMAg6cOTeiTe/GCyo+ipLB7KxvI2TgkjkXLj752T6qSiel9Bv75sPNrJxTPChlzoUwrSCXf7xzOXsaOvjGL/cOOLf/dCeBkGHRBUYC3ntlFQA/GsL7aezsIxgyzCqyPJ8efzCSTXsiouKjKBnM+3Z+s0srLlx8AHu4KDWeT31rNx+e6mBtTeI5nJFyw+Lp3H/NXH7y7nEe/8/DkfK366xMXZdVX9jzmV2Ux8eXz+KnW49HNsKLJRzpNqsoL7IhXaqe5VhAxUdRMpitR1upmTaZwvz421EPl0Uzp3C0uSslW2qHF2vefvGsUbcN8JVbFvLRi2fwjV9+yI/fOUIoZNi4o4Fl5YWD9uwZCX987Ty6+4P86HdH454PD8nNLsqN2g1VxUdRlAlGrz/IlsPNo+pJXFZVgjFQO8rZnEMhw8+21bO8vDBpWp2R4vUIj3x6BTcsms7f/cderv72Jvaf6eS+tdWjYv+i6QXcvGQ6P37nCJ1xFpCGxWdmYR6TcnxUluSz/4yKj6IoE4ytR1roC4S45qLR21zxksoisn0e3h1GPjMnvHngLHWN5/jcVVWjajeWHJ+Xx+5Zyd9+bDFzSvN58NaFfHz56HlaX/q9+XT0BvjJu4PX/TS09VCUnxXJyLBgRoF6PoqiTDzeOnCWbJ+HNdWD96MZKblZXi6pKOLtutETn0AwxMOv7GN2UR63L0vNkFs0Ho/wh1dV89MvrOGPPjLvgkPQo1lWXsTVNVN54u3DgyLrGtp6mVV4PlP2whkFHGnqoi+QukW7bqLioygZiDGGX+89w5q5peRlj266/5uWzODDUx0cPht/3cxwefKdI+w73cnXbltE1ggWeo41Hvi9+TSd6+fZ9wYmNT3R0k158XnxWTCjgGDIJFx/NN4Z//+TiqIMm90n2zne0s3tF88cddu3XTwTEXhxR8MF29pyuJlvv7KfmxZP59alM0ahd+6zurqEVXOK+eGbh+gPhADry0B9a8+ADevCufYm6tCbio+iZCD/sbOBLK9w85LRf6HPKMzlmpoy/m3LsQvK81Z7tIX7nqqlsiSf//2p5aM6/OUmIsKXrptPQ3svv9hxEoDmrn56/EEqojyfOaWTyPZ6JmzQgYqPomQYgWCIjTsbuLqmbNRCrGP542vn0XSun6cShBUPhTGGn245zmf+ZQtlBTn89AtrItsjTBSuvaiMJbOm8IM3DhEIhiJDa1VRaXyyvB7mTZuc0owRbqLioygZxm8/PMOZjj7uWl2ZsjYury7hhkXTeOS3Bzg4jG/udY2d3PPEVv7XC7tZM6+Uf//jK5kxgmzSYx0R4b9dV8Phpi42vHeCD07Gz1m3dNYUdp9sx5iJl2ZHxUdRMoyn3z3G7KI8rls4LWVtiAjfuONiCnKz+OyTWyMv13iEQoYth5v54tPbuPGRt9hV38bXP76EH33uMkomZaesj25z85LpXF5dwv//6/3825bjVJXmU1YwMLP4JZXFtHT1Rza0m0gM3uJPUZQJy96GDt6pa+Yvbl6AN8lmbBfKjMJc/vXzq/nDH73HHd9/h48vn8W1C6cxY0oufYEg9a097DjexhsHGjnT0UdhXhZ/cu08Pn9VNaVJdimdCIgID61byicefYcjTV385S0LBtW5pLIIgO3H2yJZwycKKj6KkkH802sHKcj1cfeaOWlpb9HMKbz8p1fz3d8e4Plt9fx8+8kB5wvzsri8uoTbls3khkXT4255PZFZMKOAFx9Yy8EzndwUJ/jjoukF5Gd72X68lTsume1CD1NHZv1PK0oGs7ehg1f2nOZPr69J6wR+yaRsvr5uKX91+2IOnOmkpaufHJ+XGVNyqSjJmzBRbCNl/rTJzJ82Oe45r0dYXl7E9hNtae5V6lHxUZQMwBjDQy/toTAvi8+PUq6y4ZLl9YzaJnCZxIrKIv7lLSsjQm7W6C4IdhMNOFCUDOClXad493ALf37zggkXtjzRubSymEDIsHOCeT8qPooywWnr7ucbv9zL0tlT+EwKw6uV1HBZVTEisOXI6GYKdxsVH0WZwBhj+F8v7Kalq59vfWJZyiPclNGnKD+bRTOmjHqmcLdxJD4icouI7BeROhH5apzzOSLyrH1+i4hURZ170C7fLyI3J7MpItW2jYO2zeyh2hCRG0Vkm4jstv+9bqQPQ1EmGj+rrefl3af58o0LWDpb51vGK2vmlvL+8dYJleE6qfiIiBf4PnArsBi4S0QWx1S7D2g1xswHHgEetq9dDKwHlgC3AI+KiDeJzYeBR4wxNUCrbTthG0AT8DFjzMXAvcDTw3sEijIx2Xaslb/6xQdcOa+U+6+Z63Z3lAvg8rkl9PpD7KpPvFh3vOHE81kN1BljDhtj+oENwLqYOuuAp+zj54HrxYqfXAdsMMb0GWOOAHW2vbg27Wuus21g27xjqDaMMduNMeH0uXuAXBGZ+CvUFGUITrb18EdPb2NmUS7f/8ylOtw2zrm8ugQRePfQxBl6cyI+s4HojSfq7bK4dYwxAaAdKB3i2kTlpUCbbSO2rURtRPNJYLsxpi/2JkTkfhGpFZHas2fPJrllRRm/NHb2cs/jW+jzB3ni3lUUT+AUNZlCUX42C2dM4d0jmSU+8b4yxWa5S1RntMqT9kNElmANxf1RnHoYYx4zxqwyxqwqKxu9bYMVZSzR0tXPPY9v5XRHLz/+/GXMn1bgdpeUUWLN3BK2HWuN7AE03nEiPvVARdTnciB2l6hIHRHxAYVAyxDXJipvAopsG7FtJWoDESkHXgA+a4w55OCeFGXCcaq9h7see5ejzV08/tlVrJxT4naXlFFkzdxSev0hdtZPjPU+TsTnPaDGjkLLxgog2BhTZyPWZD/AncDrxsoBvhFYb0eqVQM1wNZENu1rNtk2sG2+OFQbIlIE/BJ40BjzznBuXlEmCnWNnXzy0d9xsq2HH33uMq6cP9XtLimjTHjeZ/MEmfdJKj72/MoDwKvAh8Bzxpg9IvKQiHzcrvYEUCoidcCXga/a1+4BngP2Aq8AXzLGBBPZtG19BfiybavUtp2wDdvOfOCvRWSH/ZO6XPGKMsZ4+2ATd/5gM/1Bw4b716jwTFDC630mivjIRNykKBmrVq0ytbW1bndDUS4IYww/fOsw335lH/OnTebxz15GZWm+291SUsjfv7SXp989xq6/vcmVPG8iss0Ys2o0bGmGA0UZh7T3+PnST9/nW7/ax61LZ/LCn1ylwpMBXDG3lP5AiO3Hx/+8j2a1VpRxxuZDzfzP53ZwprOPB29dyP3XzM34bQkyhdVzS/AIbD7czBXzYleajC9UfBRlnNAXCPKd3xzgsbcOM6ckn+e/eAWXVBa73S0ljUzJzWLp7EJrsemNbvfmwlDxUZRxwLZjLTz4890cOHOOu1ZX8le3Lcq4XT8ViyvmlvLkO0fo6Q+Slz1+9/fROR9FGcO09/j52gu7+eQ/b+Zcb4AnP7eK/+8TF6vwZDBr5pXiDxq2HWt1uysXhP4GK8oYxBjDL3ef4uv/sZfmc33ct7aaL994kYqOwmVVJXg9wubDTaytGb9h9fqbrChjjN317fz9S3vZerSFJbOm8OS9l3FxuW6HoFhMzvGxrLxw3K/3UfFRlDFCY0cv//vV/Tz/fj0l+dn8w+9fzKcvq9CM1MogrphbymNvHaarLzBuveHx2WtFmUB09wd48u0jPPrGIfzBEF+4ei4PXDefKblZbndNGaNcMa+UR984xHtHW7h2wfhM6KLioygu0RcI8syW43xv0yGazvVx4+LpfO2ji6iaOsntriljnFVzSsjyCpsPN6v4KIrijEAwxM+3n+S7vz3IybYeLq8u4Yf3XKpZqBXH5GV7WVFRNK43l1PxUZQ0EQoZXv7gFN/5zQEOn+1ieXkh3/rkxaydP1UzFCjD5oq5pXxvUx0dvf5xOUSr4qMoKSYQDPHSrlN8f1MdBxvPcdH0yfzwnpXctHi6io4yYtbMK+WfXq/jvSMtXL9outvdGTYqPoqSIvoDIV7YXs+jbxziWHM3C6YX8N31K7h92SyNYFMumEsri8n2edh8qFnFR1EU6PUHea72BD944xAN7b1cPLuQH96zkhsXTcejoqOMErlZXlZWFrP58Pic91HxUZRRor3bz79tPcaP3jnK2c4+Vs4p5h8+cTEfuahMh9eUlHDFvFIe+e0B2rr7KcrPdrs7w0LFR1EukBMt3Tzx9hGeqz1Bd3+QtfOn8k/rL2HN3BIVHSWlXDGvlO/8BrYcaeHmJTPc7s6wUPFRlBGy/Xgrj//nEX71wSk8Inx8+Sz+69VzWTxrittdUzKE5eVF5GV52XyoWcVHUSYyoZDhtx+e4V/+8zDvHW2lINfHF66Zy+eurGJmYZ7b3VMyjGyfh1VVxbw7Dud9VHwUxQHtPX6e31bP05uPcrS5m9lFefz17Yv59GUVTB6nubWUicHKOcV897WDdPb6KRhH6330r0ZRhmDf6Q7+dfMxXnj/JD3+IKvmFPM/b1rArUtn4PPqdliK+6yaU4IxsONEG1fXlLndHceo+ChKDIFgiN/sPcOPf3eULUdayPF5WLdiFp+9ooqls3VrA2VssbyiEI9A7dFWFR9FGY80netjw9bj/NuW45xq76W8OI8Hb13If1lVQfGk8RXGqmQOBblZLJgxhfePj6+dTVV8lIzGGMPmw808s/UEr35wmv5giKtrpvLQuqVct3CaZiJQxgUr5xTxi+0NBENm3PzOqvgoGUnTuT6e31bPhq3HOdrcTWFeFp+5vJK718xh/rTJbndPUYbFqjkl/OTd4xw408mimeMj1F/FR8kYQiHDO4ea2LD1BL/eexp/0LC6uoQ/vaGGW5fOJDfL63YXFWVErJxTDEDtsVYVH0UZKzR29vKz2no2vHecEy09FOVnce8VVaxfXcH8aQVud09RLpjy4jzKCnJ4/1gr96yZ43Z3HOFIfETkFuC7gBd43BjzrZjzOcC/AiuBZuDTxpij9rkHgfuAIPDfjTGvDmVTRKqBDUAJ8D5wjzGmP1EbIlIKPA9cBvzYGPPACJ+FMoHoD4TYtL+R57fVs2lfI4GQYc3cEv78pgXcvGSGejnKhEJEWFlZzLZj4yfoIKn4iIgX+D5wI1APvCciG40xe6Oq3Qe0GmPmi8h64GHg0yKyGFgPLAFmAb8VkYvsaxLZfBh4xBizQUR+YNv+50RtAL3AXwNL7R8lg9nT0M7z2+p5cUcDLV39lBXkcN/aaj59WQVzy3QuR5m4rJxTzCt7TtPY2cu0gly3u5MUJ57PaqDOGHMYQEQ2AOuAaPFZB/ydffw88D2xMiquAzYYY/qAIyJSZ9sjnk0R+RC4DviMXecp2+4/J2rDGNMFvC0i84dx38oEoulcH7/YfpLnt9Wz73Qn2V4PNy6Zzp0ry7l6/lRdDKpkBCurrHmf94+1csvSmS73JjlOxGc2cCLqcz1weaI6xpiAiLQDpXb5uzHXzraP49ksBdqMMYE49RO10eTgHpQJRnhY7We19byx3xpWW15RxN/fsZSPLZs57tLLK8qFsmTWFLK8wo4T7RNGfOIFjRuHdRKVx/sqOlR9p/1IiIjcD9wPUFlZ6fQyZQxhjOH94628uKOBl3adoqWrn2kFOdx3dTV3XlpOzXQNHlAylxyfl0Uzp7DzRJvbXXGEE/GpByqiPpcDDQnq1IuIDygEWpJcG6+8CSgSEZ/t/UTXT9SGI4wxjwGPAaxatcqxaCnuc/BMJ7/YcZIXdzRQ39pDbpaHGxZN55M6rKYoA1heXsQL208SCpkxv2uuE/F5D6ixo9BOYgUQfCamzkbgXmAzcCfwujHGiMhG4Kci8h2sgIMaYCuWFzPIpn3NJtvGBtvmi0O1MbLbVsY6DW09/MfOBn6xo4EPT3Xg9QhXzZ/Kl2+8iJuWzNBM0ooSh+UVRTz97jEON50b88sIkv4F2/MrDwCvYoVFP2mM2SMiDwG1xpiNwBPA03ZAQQuWmGDXew4rOCEAfMkYEwSIZ9Nu8ivABhH5BrDdtk2iNmxbR4EpQLaI3AHcFBONp4wD2rr7eXn3aV7ccZKtR1swBlZUFPF3H1vMbctmUVaQ43YXFWVMs6LCSny740T7mBcfyUTnYdWqVaa2ttbtbihAZ6+f1z5s5KVdp3jzQCP+oGFu2STuWDGbdStmMad0kttdVJRxQyhkWPb1X/P7l8zm7+8Y/ZUnIrLNGLNqNGzp2IWSdsKC88vdp3jzwFn6AyFmTMnl3iuquOOS2SyZNQUrUl9RlOHg8QgXzy5kV/3YDzpQ8VHSQiLB+YPLK7l92UwuqSge8xOkijIeWF5RxBNvH6YvECTHN3Yzeaj4KCnjXF+A1z48Yw+pqeAoSjpYUVGIP2j48FQnKyqK3O5OQlR8lFGlvdvPa/vO8MoHp3kjRnBuu3gml1aq4ChKKllWbgnOrvo2FR9lYnO6vZdf7z3Nq3tOs+VwC4GQUcFRFJeYWZhLWUEOO0608dkr3O5NYlR8lBFR13iOV/ec5td7TrOzvh2AeWWT+MI1c7l5yQyWzS5UwVEUFxARlpcXjflMByo+iiNCIcPO+jZ+vfcMr+45zeGzXYA1ufkXN1vbFOgOoIoyNlhRUchr+87Q0etnSm6W292Ji4qPkpBef5AtR1r47d4z/Hrvac509OHzCGvmlvK5K6u4cfF0Zhbmud1NRVFiWFZehDGwu76dq+ZPdbs7cVHxUQbQ2NHLpv2NvPZhI2/XNdHdHyQvy8tHLirj5qXTuW7BdArzx+Y3KUVRLJaVW5kOdta3qfgoY5NQyLCnoYPX9p3h9X2N7LLnb2YV5vLJS8u5btE0rphbqjt/Kso4oig/m6rSfHadaHe7KwlR8clAuvsDvFPXzGsfWoLT2NmHCFxiz99cv2gaC6YXaJYBRRnHLK8oYusRx4n/046KT4ZwtKmLNw+cZdP+Rn53qJn+QIiCHB/XXFTGdQunce2CMkona+JORZkoLCsv4sUdDTR29DJtytjbVlvFZ4LS1Rdg86Fm3jxwlrcOnuVYczcAVaX53LNmDtcvnMaqqhKyfboXjqJMRMIZrnfWt3PjYhUfJUUYY9h3upM3D5zlzf1nqT3Wgj9oyM/2csXcUu5bW801NWVUTdUs0YqSCSyeWYjXI+w80caNi6e73Z1BqPiMY1q7+nm7rsnybg6cpbGzD4CFMwr4/FXVfOSiMlZWFY/p5IKKoqSGvGwvC6YXsHOMZrhW8RlH9PqDvH+slXcONfFOXTO76tsIGSjMy2JtzVQ+clEZ19SUMaNw7LnYiqKkn+UVhfxy1ymMMWMugEjFZwwTDBl2n2znnbomfneoidqjrfQFQng9wvLyQv7bdTV8ZEEZy8uL8GoqG0VRYlheXsQzW09wtLmb6jE25K7iM4YwxnDo7DneqWvmnbomNh9uprM3AFhDaX9w+Ryuml/K6uoSCsZoygxFUcYOyyvOZ7hW8VEG0NDWY3s2luCE520qSvK47eKZXDl/KlfOK2WqhkErijJMaqZNJjfLw44TbaxbMdvt7gxAxSfNNJ3r493DzfzuUDObDzVzpMlK0Dl1cjZXzJvKVfNKuWr+VCpK8l3uqaIo4x2f18PFswvHZIZrFZ8U097jZ+uRFn53qInNh5rZd7oTgMk5Pi6vLuEPLq9kbc1UzSigKEpKWFZexE/ePYY/GCLLO3bW9an4jDLRudI27T/LbjsiLcfn4bKqEv7i5llcOa+Ui2cX4htDvwiKokxMllcU8cTbR9h/upOlswvd7k4EFZ9RwBjDjhNt/Pz9k7yy5zRn7Vxpy8uLeOD35nPl/KlcUlmk620URUk7y+0M17vq21V8JgrBkOGlXQ08uukQ+890kuPzcP2iaVy/cLrmSlMUZUxQWZJPUX4WO0+08ZnLK93uTgQVnxFypKmLLz+3g+3H26iZNplvfeJiPrps5pjdNVBRlMxERFhWXjTmMh2o+IyAPQ3t3P34Fgzwnf+ynDtWzMajizwVRRmjrCgv5HubztLdHyA/e2y89nXGe5h09vr5o6e3kZvl5Rd/chWfuLRchUdRlDHNisoiQga2Hx873o+KzzD58TtHqW/t4XufuUQzRCuKMi5YXV2KzyO8dfCs212J4Eh8ROQWEdkvInUi8tU453NE5Fn7/BYRqYo696Bdvl9Ebk5mU0SqbRsHbZvZI20jFfz7+/VcXTOVlXNKUtmMoijKqDE5x8fKOcW8daDJ7a5ESCo+IuIFvg/cCiwG7hKRxTHV7gNajTHzgUeAh+1rFwPrgSXALcCjIuJNYvNh4BFjTA3QatsedhvDfRBOaDrXx9Hmbj5yUVkqzCuKoqSMjywo48NTHZxo6Xa7K4Azz2c1UGeMOWyM6Qc2AOti6qwDnrKPnweuF2u5/jpggzGmzxhzBKiz7cW1aV9znW0D2+YdI2xj1GnssPKulRfnpcK8oihKyrhjxWx8HuFbv9rndlcAZ9Fus4ETUZ/rgcsT1THGBESkHSi1y9+NuTac3S6ezVKgzRgTiFN/JG1EEJH7gfvtj+dEpBkYkQ9668MjuWpMM5URPosJiD4LC30O55lQz+JR4NG7R3TpVGDOaPXDifjEC+UyDuskKo/ncQ1VfyRtDCww5jHgsfBnEak1xqyKc23Goc/iPPosLPQ5nEefhYX9HKpGy56TYbd6oCLqcznQkKiOiPiAQqBliGsTlTcBRbaN2LaG24aiKIoyRnEiPu8BNXYUWjbW5P7GmDobgXvt4zuB140xxi5fb0eqVQM1wNZENu1rNtk2sG2+OMI2FEVRlDFK0mE3e37lAeBVwAs8aYzZIyIPAbXGmI3AE8DTIlKH5Y2st6/dIyLPAXuBAPAlY0wQIJ5Nu8mvABtE5BvAdts2I2kjCY8lr5Ix6LM4jz4LC30O59FnYTGqz0Es50FRFEVR0odmOFAURVHSjoqPoiiKknYyUnySpQuaCIjIkyLSKCIfRJWViMhv7NRFvxGRYrtcROSf7OexS0QujbrmXrv+QRG5N15bYxkRqRCRTSLyoYjsEZE/tcsz6lmISK6IbBWRnfZz+LpdPqbTWaUSO9vKdhF5yf6ckc9CRI6KyG4R2SEitXZZ6v8+jDEZ9YMV4HAImAtkAzuBxW73KwX3eQ1wKfBBVNm3ga/ax18FHraPPwr8CmvN1Bpgi11eAhy2/y22j4vdvrdhPoeZwKX2cQFwACulU0Y9C/t+JtvHWcAW+/6eA9bb5T8A/tg+/hPgB/bxeuBZ+3ix/TeTA1Tbf0tet+9vhM/ky8BPgZfszxn5LICjwNSYspT/fWSi5+MkXdC4xxjzFlZUYDTRKYpiUxf9q7F4F2ut1UzgZuA3xpgWY0wr8Bus/HnjBmPMKWPM+/ZxJ/AhVgaMjHoW9v2csz9m2T+GMZzOKpWISDlwG/C4/XlMp/ZygZT/fWSi+MRLFzQoHc8EZbox5hRYL2Vgml2e6JlMqGdlD5dcgvWtP+OehT3MtANoxHo5HMJhOisgOp3VuH4ONv8H+EsgZH92nNqLifcsDPBrEdkmVhoySMPfx9jY0i69OErHk2FcUOqi8YCITAb+HfgfxpgO64tr/KpxyibEszDW+rcVIlIEvAAsilfN/nfCPgcRuR1oNMZsE5Frw8Vxqk74Z2FzlTGmQUSmAb8RkaEyj47as8hEzyeT0/GcsV1k7H8b7fLhpkEaV4hIFpbw/Jsx5ud2cUY+CwBjTBvwBtaYfSams7oK+LiIHMUadr8OyxPKxGeBMabB/rcR60vJatLw95GJ4uMkXdBEJTpFUWzqos/akSxrgHbb1X4VuElEiu1ol5vssnGDPTb/BPChMeY7Uacy6lmISJnt8SAiecANWPNfGZfOyhjzoDGm3FhJMtdj3dsfkIHPQkQmiUhB+Bjr9/oD0vH34XakhRs/WBEbB7DGvL/mdn9SdI/PAKcAP9a3kvuwxqlfAw7a/5bYdQVrc79DwG5gVZSdz2NNpNYBf+j2fY3gOazFcv93ATvsn49m2rMAlmGlq9plv1z+xi6fi/XCrAN+BuTY5bn25zr7/NwoW1+zn89+4Fa37+0Cn8u1nI92y7hnYd/zTvtnT/h9mI6/D02voyiKoqSdTBx2UxRFUVxGxUdRFEVJOyo+iqIoStpR8VEURVHSjoqPoiiKknZUfBRFUZS0o+KjKIqipJ3/ByGtpoS8GSuHAAAAAElFTkSuQmCC\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": 14,
"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": 12,
"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": 13,
"metadata": {},
"outputs": [
{
"ename": "KeyboardInterrupt",
"evalue": "",
"output_type": "error",
"traceback": [
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[1;31mKeyboardInterrupt\u001b[0m Traceback (most recent call last)",
"\u001b[1;32m<ipython-input-13-95555ba396d5>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m 1\u001b[0m \u001b[0mtotal_f\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mupdate_integration_options\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdraws_per_dim\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m20000000\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mmc_sampler\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mNone\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 2\u001b[0m \u001b[0minte\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mtotal_f\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mintegrate\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mlimits\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m(\u001b[0m\u001b[1;36m4250\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m4600\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnorm_range\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mFalse\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 3\u001b[1;33m \u001b[0minte_fl\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mzfit\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mrun\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0minte\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 4\u001b[0m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0minte_fl\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 5\u001b[0m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mpdg\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m\"jpsi_BR\"\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m/\u001b[0m\u001b[0mpdg\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m\"NR_BR\"\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0minte_fl\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mpdg\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m\"psi2s_auc\"\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m/\u001b[0m\u001b[0mpdg\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m\"NR_auc\"\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32m~\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\util\\execution.py\u001b[0m in \u001b[0;36m__call__\u001b[1;34m(self, *args, **kwargs)\u001b[0m\n\u001b[0;32m 79\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 80\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m__call__\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 81\u001b[1;33m \u001b[1;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msess\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mrun\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 82\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 83\u001b[0m \u001b[1;31m# def close(self):\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32m~\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\client\\session.py\u001b[0m in \u001b[0;36mrun\u001b[1;34m(self, fetches, feed_dict, options, run_metadata)\u001b[0m\n\u001b[0;32m 927\u001b[0m \u001b[1;32mtry\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 928\u001b[0m result = self._run(None, fetches, feed_dict, options_ptr,\n\u001b[1;32m--> 929\u001b[1;33m run_metadata_ptr)\n\u001b[0m\u001b[0;32m 930\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0mrun_metadata\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 931\u001b[0m \u001b[0mproto_data\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mtf_session\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mTF_GetBuffer\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mrun_metadata_ptr\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32m~\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\client\\session.py\u001b[0m in \u001b[0;36m_run\u001b[1;34m(self, handle, fetches, feed_dict, options, run_metadata)\u001b[0m\n\u001b[0;32m 1150\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0mfinal_fetches\u001b[0m \u001b[1;32mor\u001b[0m \u001b[0mfinal_targets\u001b[0m \u001b[1;32mor\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mhandle\u001b[0m \u001b[1;32mand\u001b[0m \u001b[0mfeed_dict_tensor\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 1151\u001b[0m results = self._do_run(handle, final_targets, final_fetches,\n\u001b[1;32m-> 1152\u001b[1;33m feed_dict_tensor, options, run_metadata)\n\u001b[0m\u001b[0;32m 1153\u001b[0m \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 1154\u001b[0m \u001b[0mresults\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32m~\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\client\\session.py\u001b[0m in \u001b[0;36m_do_run\u001b[1;34m(self, handle, target_list, fetch_list, feed_dict, options, run_metadata)\u001b[0m\n\u001b[0;32m 1326\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0mhandle\u001b[0m \u001b[1;32mis\u001b[0m \u001b[1;32mNone\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 1327\u001b[0m return self._do_call(_run_fn, feeds, fetches, targets, options,\n\u001b[1;32m-> 1328\u001b[1;33m run_metadata)\n\u001b[0m\u001b[0;32m 1329\u001b[0m \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 1330\u001b[0m \u001b[1;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_do_call\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0m_prun_fn\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mhandle\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfeeds\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfetches\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32m~\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\client\\session.py\u001b[0m in \u001b[0;36m_do_call\u001b[1;34m(self, fn, *args)\u001b[0m\n\u001b[0;32m 1332\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m_do_call\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfn\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 1333\u001b[0m \u001b[1;32mtry\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1334\u001b[1;33m \u001b[1;32mreturn\u001b[0m \u001b[0mfn\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 1335\u001b[0m \u001b[1;32mexcept\u001b[0m \u001b[0merrors\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mOpError\u001b[0m \u001b[1;32mas\u001b[0m \u001b[0me\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 1336\u001b[0m \u001b[0mmessage\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mcompat\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mas_text\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0me\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmessage\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32m~\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\client\\session.py\u001b[0m in \u001b[0;36m_run_fn\u001b[1;34m(feed_dict, fetch_list, target_list, options, run_metadata)\u001b[0m\n\u001b[0;32m 1317\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_extend_graph\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 1318\u001b[0m return self._call_tf_sessionrun(\n\u001b[1;32m-> 1319\u001b[1;33m options, feed_dict, fetch_list, target_list, run_metadata)\n\u001b[0m\u001b[0;32m 1320\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 1321\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m_prun_fn\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mhandle\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfeed_dict\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfetch_list\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32m~\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\client\\session.py\u001b[0m in \u001b[0;36m_call_tf_sessionrun\u001b[1;34m(self, options, feed_dict, fetch_list, target_list, run_metadata)\u001b[0m\n\u001b[0;32m 1405\u001b[0m return tf_session.TF_SessionRun_wrapper(\n\u001b[0;32m 1406\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_session\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0moptions\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfeed_dict\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfetch_list\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtarget_list\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1407\u001b[1;33m run_metadata)\n\u001b[0m\u001b[0;32m 1408\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 1409\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m_call_tf_sessionprun\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mhandle\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfeed_dict\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfetch_list\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;31mKeyboardInterrupt\u001b[0m: "
]
}
],
"source": [
"# total_f.update_integration_options(draws_per_dim=20000000, mc_sampler=None)\n",
"# inte = total_f.integrate(limits = (4250, 4600), 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 = \"p4415\"\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": [
"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": [
"# 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
}
