{
"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 = 10000\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(tf.math.real(-q))))\n",
" \n",
" big_bracket = tf.where(tf.math.real(y) > ztf.constant(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) - ztf.to_complex(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) - ztf.to_complex(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": [
"# 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',\n",
" 'Dbar_mass', 'Dbar_scale', 'Dbar_phase',\n",
" 'DDstar_mass', 'DDstar_scale', 'DDstar_phase',\n",
" 'tau_mass', 'C_tt']\n",
"# 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",
" def P2_D(q):\n",
" Dbar_contrib = ztf.to_complex(self.params['Dbar_scale'])*tf.exp(tf.complex(ztf.constant(0.0), self.params['Dbar_phase']))*ztf.to_complex(h_S(self.params['Dbar_mass'], q))\n",
" DDstar_contrib = ztf.to_complex(self.params['DDstar_scale'])*tf.exp(tf.complex(ztf.constant(0.0), self.params['DDstar_phase']))*ztf.to_complex(h_P(self.params['DDstar_mass'], q))\n",
" return Dbar_contrib + DDstar_contrib\n",
" \n",
" def ttau_cusp(q):\n",
" return ztf.to_complex(self.params['C_tt'])*(ztf.to_complex((h_S(self.params['tau_mass'], q))) - ztf.to_complex(h_P(self.params['tau_mass'], q)))\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) + P2_D(x) + ttau_cusp(x)\n",
"\n",
" vec_f = vec(x, 0.0 * 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) #lower_limit = rho_mass - rho_width,\n",
"# upper_limit = rho_mass + rho_width)\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) #, lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"rho_s = zfit.Parameter(\"rho_s\", ztf.constant(rho_scale), floating = False) #, lower_limit=rho_scale-np.sqrt(rho_scale), upper_limit=rho_scale+np.sqrt(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) #, lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"omega_s = zfit.Parameter(\"omega_s\", ztf.constant(omega_scale), floating = False) #, lower_limit=omega_scale-np.sqrt(omega_scale), upper_limit=omega_scale+np.sqrt(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) #, lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"phi_s = zfit.Parameter(\"phi_s\", ztf.constant(phi_scale), floating = False) #, lower_limit=phi_scale-np.sqrt(phi_scale), upper_limit=phi_scale+np.sqrt(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), lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"jpsi_s = zfit.Parameter(\"jpsi_s\", ztf.constant(jpsi_scale), lower_limit=jpsi_scale-np.sqrt(jpsi_scale), upper_limit=jpsi_scale+np.sqrt(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), lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"psi2s_s = zfit.Parameter(\"psi2s_s\", ztf.constant(psi2s_scale), lower_limit=psi2s_scale-np.sqrt(psi2s_scale), upper_limit=psi2s_scale+np.sqrt(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) #, lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"p3770_s = zfit.Parameter(\"p3770_s\", ztf.constant(p3770_scale), floating = False) #, lower_limit=p3770_scale-np.sqrt(p3770_scale), upper_limit=p3770_scale+np.sqrt(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) #, lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"p4040_s = zfit.Parameter(\"p4040_s\", ztf.constant(p4040_scale), floating = False) #, lower_limit=p4040_scale-np.sqrt(p4040_scale), upper_limit=p4040_scale+np.sqrt(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) #, lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"p4160_s = zfit.Parameter(\"p4160_s\", ztf.constant(p4160_scale), floating = False) #, lower_limit=p4160_scale-np.sqrt(p4160_scale), upper_limit=p4160_scale+np.sqrt(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) #, lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"p4415_s = zfit.Parameter(\"p4415_s\", ztf.constant(p4415_scale), floating = False) #, lower_limit=p4415_scale-np.sqrt(p4415_scale), upper_limit=p4415_scale+np.sqrt(p4415_scale))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Dynamic generation of 2 particle contribution"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<zfit.util.temporary.TemporarilySet at 0x19ee031c8d0>"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"_0 = jpsi_scale*np.cos(jpsi_phase)*jpsi_width/jpsi_mass**3 + psi2s_scale*np.cos(psi2s_phase)*psi2s_width/psi2s_mass**3\n",
"_1 = p3770_scale*np.cos(p3770_phase)*p3770_width/p3770_mass**3 + p4040_scale*np.cos(p4040_phase)*p4040_width/p4040_mass**3\n",
"_2 = p4160_scale*np.cos(p4160_phase)*p4160_width/p4160_mass**3 + p4415_scale*np.cos(p4415_phase)*p4415_width/p4415_mass**3\n",
"\n",
"C_pert = np.random.uniform(0.03, 0.1)\n",
"# c_pert = 0.1\n",
"m_c = 1300\n",
"\n",
"cDDstar_phase = 10\n",
"\n",
"\n",
"DDstar_eta = 0\n",
"Dbar_phase = np.random.uniform(0.0, 2*np.pi)\n",
"DDstar_phase = np.random.uniform(0.0, 2*np.pi)\n",
"DDstar_mass = pdg['D0_M']\n",
"\n",
"if Dbar_phase < np.pi:\n",
" Dbar_phase = 0.0\n",
"else:\n",
" Dbar_phase = np.pi\n",
"\n",
"R = (C_pert/(m_c**2) - ((_0 + _1 + _2)))\n",
"\n",
"Dbar_mass = (pdg['D0_M']+pdg['Dst_M'])/2\n",
"\n",
"Dbar_eta = R/np.cos(Dbar_phase)*(6*Dbar_mass**2)\n",
"\n",
"# print(np.cos(Dbar_phase))\n",
"\n",
"# cDDstar_phase = R_*10*DDstar_mass**2/DDstar_eta\n",
"\n",
"\n",
"# print(Dbar_eta)\n",
"\n",
"\n",
"Dbar_s = zfit.Parameter(\"Dbar_s\", ztf.constant(Dbar_eta), lower_limit=-1.464, upper_limit=1.464)\n",
"Dbar_m = zfit.Parameter(\"Dbar_m\", ztf.constant(Dbar_mass), floating = False)\n",
"Dbar_p = zfit.Parameter(\"Dbar_p\", ztf.constant(Dbar_phase), floating = False)\n",
"DDstar_s = zfit.Parameter(\"DDstar_s\", ztf.constant(DDstar_eta), floating = False)\n",
"DDstar_m = zfit.Parameter(\"DDstar_m\", ztf.constant(DDstar_mass), floating = False)\n",
"DDstar_p = zfit.Parameter(\"DDstar_p\", ztf.constant(DDstar_phase), floating = False)\n",
"\n",
"Dbar_s.set_value(0.0)\n",
"DDstar_s.set_value(0.0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Tau parameters"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"tau_m = zfit.Parameter(\"tau_m\", ztf.constant(pdg['tau_M']), floating = False)\n",
"Ctt = zfit.Parameter(\"Ctt\", ztf.constant(0.0), lower_limit=-0.5, upper_limit=0.5)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup pdf"
]
},
{
"cell_type": "code",
"execution_count": 11,
"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",
" DDstar_mass = DDstar_m, DDstar_scale = DDstar_s, DDstar_phase = DDstar_p,\n",
" Dbar_mass = Dbar_m, Dbar_scale = Dbar_s, Dbar_phase = Dbar_p,\n",
" tau_mass = tau_m, C_tt = Ctt)\n",
" \n",
" \n",
"# print(total_pdf.obs)\n",
"\n",
"# print(calcs_test)\n",
"\n",
"# for param in total_f.get_dependents():\n",
"# print(zfit.run(param))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Test if graphs actually work and compute values"
]
},
{
"cell_type": "code",
"execution_count": 18,
"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": 19,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAD4CAYAAAA+epuFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXxV1bXA8d8iIQNTAkmYEiABwhAUFMIk2mrVgnWgVVoiRVFArIUO9vVZea+tQ/va2jrU1jEKilYJqK1Gq2KRWiemIINMgRCmMCUkEAghw03W++MeMMaEXEKSc+/N+n4++XDvvvvss87Rm5V99j5ni6pijDHGNIc2bgdgjDEmeFmSMcYY02wsyRhjjGk2lmSMMcY0G0syxhhjmk2o2wE0p9jYWE1MTHQ7DGOMCShr1qw5rKpxTdFWUCeZxMREsrKy3A7DGGMCiojsbqq27HKZMcaYZmNJxhhjTLOxJGOMMabZBPWYTF0qKyvJy8ujrKzM7VDOWUREBAkJCbRt29btUIwxpk6tLsnk5eXRsWNHEhMTERG3w2k0VaWwsJC8vDySkpLcDscYY+rU6i6XlZWVERMTE9AJBkBEiImJCYoemTEmeLW6JAMEfII5JViOwxgTvFrd5TJjjDlyooKcghJ2Hj7B0dIKTlZUExbahuh2benTpR1DekYR1c7GOpuCT0lGRCYAjwIhwLOq+odan4cDLwAjgEJgsqrucj6bC8wAqoAfq+oSH9v8K3CrqnZoaB/B6IMPPuDBBx/krbfeory8nKuvvprDhw8zd+5cJk+e7HZ4xgQUT1U1n+4o5L3NB1mRW0ROfkmD2wzp2YkJQ7rz3dRedI+KaIEog1ODSUZEQoDHgSuBPGC1iGSq6uYa1WYAR1S1v4ikAQ8Ak0UkBUgDhgA9gaUiMsDZpt42RSQViK4VSp37aNRRB5i1a9dSWVnJunXr3A7FmICy6/AJFizfxZvr93O4pIJ2YSGMTOzC9cPjGdyjE0kx7YntGE5EaBsqq5Si0gpyC0rYkFfM+1sO8dC/tvHo+9u57oKe/PybA+kZHen2IQUcX3oyo4AcVc0FEJEMYCJQM8lMBO51Xr8KPCbeAYOJQIaqlgM7RSTHaY/62nSS2p+AKcB3GtqHBuDSnrt27WLChAmMHj2atWvXMmDAAF544QU+/PBDfvrTnxIbG8vw4cMByM/PZ+rUqRQUFHDBBRfw2muv0a9fP5ePwBj/tmZ3EU9+kMv7Ww/Rtk0brkzpxrXDenLpwDgi2obUuU1oCMSHRRIfHcklyXHMvqw/e4tKmf/JTl5euYd3Pj/IT69I5rZL+tKmjY2H+sqXJBMP7K3xPg8YXV8dVfWISDEQ45SvqLVtvPO6vjbnAJmqeqDWwHZ9+zhcs5KIzAJmAfTu3fuMB3bfm5vYvP/YGeucrZSenbjn2iEN1svOzmbevHmMGzeO6dOn8/DDD/P000+zbNky+vfvf/qSWNeuXXn22WdPXzozxtRv68FjPLgkm6Vb8unSPowfXdafqWP70LVj4y539erSjnuuHcL0cUnc/9Zmfv/OVj7OOcwjky8gtkN4E0cfnHyZXVZXyq7de6ivzlmVi0hP4LvAXxsZB6qarqqpqpoaF9ckDxFtFr169WLcuHEATJ06laysLJKSkkhOTkZEmDp1qssRGhM4ik5UcNer67nq0Y9YubOIuyYM5ONfXMbPvjmw0Qmmpl5d2pF+0wh+f/35rNpZxA1PfsruwhNNEHnw86Unkwf0qvE+AdhfT508EQkFooCiBratq/xCoD+Q4/Ri2olIjqr2P8M+Gs2XHkdzqT39uLi42KYkG3OWqquVV9bs5ffvbKWkzMNtl/Tlh5f2I7pdWJPvS0S4cVRvBnXvyPTnV3PDk8tZeNtokrt1bPJ9BRNfejKrgWQRSRKRMLwD+Zm16mQC05zXk4BlzlhJJpAmIuEikgQkA6vqa1NV/6mq3VU1UVUTgVInwZxpHwFpz549LF++HICFCxdyxRVXsHPnTnbs2HG6zBhTv92FJ5icvpxfvPY5A7p25O2fXML/fGtwsySYmi7s3ZlXfnARInDz/FXsO3qyWfcX6BpMMqrqwTtOsgTYAixW1U0icr+IXOdUmwfEOAP7PwPudrbdBCzGO0ngXWC2qlbV12YDodS5j0A1ePBgFixYwNChQykqKuLOO+8kPT2dq6++mosvvpg+ffq4HaIxfklVWbhqD1c9+hHZB4/zx0lDWXT7GAa0YI+if9cOLLh1FCXlHm6et5LjZZUttu9AIwHcGWhQamqq1l60bMuWLQwePNiliLx27drFNddcw8aNG8+5LX84HmNaSsHxcub+fQNLt+RzUb8YHvzuMFenFS/fUcjUeSu5YnBXnpo6ImgueYvIGlVNbYq2WuVjZYwxgWdFbiFXPfoRH24/zK+uSeFvM0a7ft/K2H4xzL1qEEs2HeKZj3JdjcVf2WNlXJCYmNgkvRhjWgNV5ekPc/nTkmz6dGnH32aOYlD3Tm6HddqMi5NYvauIB5ds4+sDujKwu00EqKlV9mSC5RJhsByHMfUpPlnJrBfX8Id3tjJ+SDfemDPOrxIMeGed/d93zqdDRCj//ep6PFXVbofkV1pdkomIiKCwsDDgf0GfWk8mIsKeqWSCU07+ca577GP+vTWfX1+TwuNThtMxwj8fWhnbIZz7Jw5hQ14xz3+6y+1w/Eqru1yWkJBAXl4eBQUFbodyzk6tjGlMsPnPtgLmvPQZ4W3bkDFrDKmJXdwOqUHXDO3Jq2vyeHTpdr59Ybw9EcDR6pJM27ZtbSVJY/yUqrLg013c/9ZmBnbvxLPTUokPoIdS/vLqFCb8+UMeei+b318/1O1w/EKru1xmjPFPlVXV/PL1jdz75mYuH9yNV38wNqASDHjvn5l2USIZq/eSffC42+H4BUsyxhjXHS+r5NbnVvPSyj3ccWk/np46gvbhgXmhZc5l/WkfFspflm13OxS/YEnGGOOq/GNlfO/pFazILeRPk4byiwmDAvpR+p3bh3HLRYm8/fkB681gScYY46IdBSV85wnvE43n3zKS76b2anijADDzkiRvb+Z9681YkjHGuOKzPUeY9OSnlHuqyJg1hq8N8N+lOc5WdLswbh7bh7c3HmBPYanb4bjKkowxpsW9v+UQU55ZQafItrx2x0UMTai92nrgm3ZRIiEiPPfpTrdDcZUlGWNMi1q8ei+zXlzDgG4dee2Oi+gT097tkJpFt04RXDO0B69k5bXqpzRbkjHGtJhnP8rlrtc2cFG/GBbeNibob1icfnESJeUeFq3e23DlIGVJxhjT7FSVR5du57f/3MJV53Vn3rSRATtF+WwMTYhmeO9oXl61J+AfZdVYlmSMMc1KVfnd21t4ZOk2Jo1I4K83XkhYaOv51ZM2sje5BSdYs/uI26G4wqf/0iIyQUSyRSRHRL6yIqWzvPIi5/OVIpJY47O5Tnm2iIxvqE0RmSci60Vkg4i8KiIdnPJbRKRARNY5PzPP5cCNMc2vqlr5n398zjMf7eSWixL54w1DCQ1pPQkG4OqhPWgXFsLirNZ5yazB/9oiEgI8DlwFpAA3ikhKrWozgCOq2h94BHjA2TYFSAOGABOAJ0QkpIE271TVYao6FNiDd5nmUxap6gXOz7ONO2RjTEuorKrmzkXrWLhqL7Mv68c916YE9E2WjdU+PJRrhvbgrQ0HKCn3uB1Oi/PlT4pRQI6q5qpqBZABTKxVZyKwwHn9KnC5eNchnQhkqGq5qu4Ecpz26m1TVY8BONtHAq3zQqYxAayssoo7/raGzPX7+cWEQfz3+EFBszRxY0we2YvSiire3nDA7VBanC9JJh6o2c/Lc8rqrKOqHqAYiDnDtmdsU0SeAw4Cg4C/1qh3Q43LaHXeGiwis0QkS0SyguFx/sYEmtIKDzMWrGbplnx+M3EId1zaz+2QXDe8d2f6xLTjzQ373Q6lxfmSZOr686N276K+Omdb7n2heivQE9gCTHaK3wQSnctoS/mi5/TlRlTTVTVVVVPj4oLnDmJjAsGJcg+3zF/N8h2FPPy9Ydw0NtHtkPyCiHD1+T34dEchhSXlbofTonxJMnlAzV5DAlA7HZ+uIyKhQBRQdIZtG2xTVauARcANzvtCVT31X+cZYIQPsRtjWkhJuYdp81exZs8RHk27kOuH24J6NV0ztCdV1co7Gw+6HUqL8iXJrAaSRSRJRMLwDuRn1qqTCUxzXk8Clql3UngmkObMPksCkoFV9bUpXv3h9JjMtcBW532PGvu7Dm8vxxjjB46XVXLzvJWs3XuUv6RdyLXDerodkt8Z3KMjfePa81Yru2TW4N1QquoRkTnAEiAEmK+qm0TkfiBLVTOBecCLIpKDtweT5my7SUQWA5sBDzDb6aFQT5ttgAUi0gnvJbX1wB1OKD8WkeucdoqAW5rkDBhjzsmxskqmzV/F53nFPD7lQiac16PhjVohEeGaoT3567Lt5B8ro2unCLdDahESzHehpqamalZWltthGBO0ik9WcvP8VWzaV8zj3x/O+CHd3Q7Jr2UfPM74P3/I775zPlNG93Y7nHqJyBpVTW2KtlrXXVHGmCZTXFrJTfNWsnl/MU9OHWEJxgcDunUgoXMky7YecjuUFmNJxhhz1o6WVvD9eSvYeuA4T00dwZUp3dwOKSCICFcM7sbHOYcpq6xyO5wWYUnGGHNWjpyoYMozK9l2qISnbxrB5YMtwZyNywd3payymk9yDrsdSouwJGOM8VnRiQqmPLuSnIISnrk5lcsGdXU7pIAzKqkL7cNCWLol3+1QWoQlGWOMTwpLypnyzApyC0qYNy2VrwfRcsktKTw0hK8NiGPZ1kOt4vH/lmSMMQ06XFLOjc+sYFfhCebfMpJLki3BnIvLBnXl0LFysg8ddzuUZmdJxhhzRgXHy7kxfQV7ikqZP20k4/rHuh1SwLvYOYcfbw/+cRlLMsaYeuUfKyMtfTl5R07y/K2juMgSTJPoGR1J37j2rWLw35KMMaZOh46VkZa+ggPFZSyYPooxfWPcDimoXNw/lpU7i6jwVLsdSrOyJGOM+YqDxd4Ec+hYGS9MH8WopC5uhxR0xvWPpbSiinV7j7odSrOyJGOM+ZL9R08yOX05BcfLeWHGKFITLcE0hzF9Y2gj8HGQXzKzJGOMOW3f0ZOkpa+gqKSCF2aMYkQfSzDNJSqyLUMTovl4e3AvrmhJxhgDQN6RUtLSl3OktIIXZ45meO/ObocU9Mb2i2FDXjGlFR63Q2k2lmSMMewtKmXy0ysoLq3kbzNGc0GvaLdDahVGJXbBU61BPS5jScaYVm5PYSlp6SsoKffw8m1jGGYJpsUM79MZEcjadcTtUJqNT0lGRCaISLaI5IjI3XV8Hi4ii5zPV4pIYo3P5jrl2SIyvqE2RWSeiKwXkQ0i8qqIdGhoH8aYxtldeIK09OWcqPDw0szRnBcf5XZIrUpUZFsGduvI6l1FbofSbBpMMiISAjwOXAWkADeKSEqtajOAI6raH3gEeMDZNgXvKplDgAnAEyIS0kCbd6rqMFUdCuwB5pxpH8aYxtl1+ARp6Ss4WVnFyzPHWIJxycjELny2+wiequC8X8aXnswoIEdVc1W1AsgAJtaqMxFY4Lx+FbhcRMQpz1DVclXdCeQ47dXbpqoeA3C2jwS0gX0YY87STifBlHuqefm2MaT07OR2SK3WyKQunKioYuvB4HyOmS9JJh7YW+N9nlNWZx1V9QDFQMwZtj1jmyLyHHAQGAT8tYF9fImIzBKRLBHJKigI7qmBxjTGjoISJj+9nMqqal6+bTSDe1iCcdPIRO8svlU7g/OSmS9Jpq7eQu3nU9dX52zLvS9UbwV6AluAyWcRB6qarqqpqpoaF2dPijWmppz8EtLSV1CtysJZYxjU3RKM23pERZLQOTJox2V8STJ5QK8a7xOA/fXVEZFQIAooOsO2DbapqlXAIuCGBvZhjPFBTv5x0tJXoAoLbxvDgG4d3Q7JOIb37sz6IJ3G7EuSWQ0ki0iSiIThHcjPrFUnE5jmvJ4ELFPvajyZQJozMywJSAZW1demePWH02My1wJbG9iHMaYB2w55E4wIZMwaQ7IlGL8yrFc0+4vLyD9W5nYoTS60oQqq6hGROcASIASYr6qbROR+IEtVM4F5wIsikoO3d5HmbLtJRBYDmwEPMNvpoVBPm22ABSLSCe/lsfXAHU4ode7DGHNm2QePM+WZFYS0ERbOGkO/uA5uh2RquaCXd2bf+rxirkyJcDmapiXB3BlITU3VrKwst8MwxjVbDx5jyjMraRsiLLxtDH0twfilkxVVnHfvEu74ej9+Pn6g2+EgImtUNbUp2rI7/o0JUpv3H+PG9BWEhbQhY9ZYSzB+LDIshIHdOrI+L/jGZSzJGBOENu0vZsqzK4hoG0LGrDEkxbZ3OyTTgGG9olm/9yjBdnXJkowxQWZD3lGmPLOS9mGhLJo1lkRLMAFhWEIUx8o87CosdTuUJmVJxpggsmZ3Ed9/ZiUdI0LJmDWG3jHt3A7J+OjUg0mDbSqzJRljgsTyHYXcNG8VsR3DWXz7WHp1sQQTSJK7diCybUjQPfbfkowxQeDDbQXc8twq4qMjWTRrDD2jI90OyZyl0JA2pPTsxOb9x9wOpUlZkjEmwL2/5RAzF2TRN64DGbPG0LVTcN1n0ZoM7tGRLQeOBdXgvyUZYwLYO58f4PYX1zC4R0cW3jaamA7hbodkzkFKjyiOl3vIO3LS7VCajCUZYwLUG+v2MWfhWob1iubFmaOJbhfmdkjmHJ1acmFTEF0ysyRjTABavHovP120jpGJnXlh+ig6RbR1OyTTBAZ260gbgc0HgifJNPjsMmOMf3lx+S5+9cYmLkmOJf2mVCLDQtwOyTSRyLAQ+sZ1CKrBf+vJGBNAnv0ol1+9sYkrBnfl2WmWYIJRSo9ObAminowlGWMCgKry2LLt/PafW/jW+d154vsjCA+1BBOMUnp2Yt/RkxwtrXA7lCZhScYYP6eq/OGdrTz43ja+c2E8f0m7kLBQ++oGqxRnOexgGZex/1ON8WNV1cr//ONznv4wl5vH9uGh7w4jNMS+tsFs8KkkEyTjMj793yoiE0QkW0RyROTuOj4PF5FFzucrRSSxxmdznfJsERnfUJsi8pJTvlFE5otIW6f8UhEpFpF1zs+vz+XAjfF3FZ5qfpKxloWr9jLnsv7cd90Q2rQRt8MyzSyuYzixHcLYdui426E0iQaTjIiEAI8DVwEpwI0iklKr2gzgiKr2Bx4BHnC2TcG7guUQYALwhIiENNDmS8Ag4HwgEphZYz8fqeoFzs/9jTlgYwLByYoqbn8xi7c2HGDuVYP4+fiBeFckN61BcteObDtU4nYYTcKXnswoIEdVc1W1AsgAJtaqMxFY4Lx+FbhcvN+IiUCGqpar6k4gx2mv3jZV9W11AKuAhHM7RGMCy7GySqbNX8UH2wr4/fXnc/vX+7kdkmlhA7p1ICe/JCgeL+NLkokH9tZ4n+eU1VlHVT1AMRBzhm0bbNO5THYT8G6N4rEisl5E3hGRIXUFKyKzRCRLRLIKCgp8ODxj/EdhSTlTnlnBZ3uO8Je0C7lxVG+3QzIuSO7WkZJyD/uLy9wO5Zz5kmTq6qPXTq/11Tnb8pqeAD5U1Y+c958BfVR1GPBX4PW6glXVdFVNVdXUuLi4uqoY45cOFJ/ke08vZ/uhEp65OZVrh/V0OyTjkgHdOgIExbiML0kmD+hV430CsL++OiISCkQBRWfY9oxtisg9QBzws1NlqnpMVUuc128DbUUk1of4jfF7uw6fYNKTyzl0rJwXpo/iskFd3Q7JuCi5awcAcoJgXMaXJLMaSBaRJBEJwzuQn1mrTiYwzXk9CVjmjKlkAmnO7LMkIBnvOEu9bYrITGA8cKOqVp/agYh0d8Z5EJFRTuyFjTloY/zJlgPHmPTUck5WVrHwtjGM7hvjdkjGZZ3bhxHbITwoejINPrtMVT0iMgdYAoQA81V1k4jcD2SpaiYwD3hRRHLw9mDSnG03ichiYDPgAWarahVAXW06u3wK2A0sd3LK352ZZJOAO0TEA5wE0jQYRsVMq7Yit5DbFmTRPjyUv80cTf+uHd0OyfiJAd06sC0/8HsyEsy/p1NTUzUrK8vtMIyp07sbD/DjjHX07tKOF6aPstUszZfcm7mJV7L2svG+8S0+fV1E1qhqalO0ZbcOG+OCl1bu5ocvfcaQnp145faxlmDMVyR368CJiir2HQ3sBcwsyRjTglSVR5du53//sZFLB3bl5Zlj6NzeFhszX3Vqhtn2AL9kZknGmBZSVa388vWNPLJ0GzcMT+Dpm0bYo/pNvU7NMNse4IP/tmiZMS2grLKKn2as491NB/nB1/vxiwn2mBhzZtHtwohpH8bOwyfcDuWcWJIxppkdK6vktgVZrNxZxK+uSWHGxUluh2QCRFJse3ILLMkYY+qRf6yMac+tZvuh4zyadgETL6j9RCZj6pcU254PtgX247FsTMaYZrL90HG+88Sn7C48wfxbRlqCMWetb1wHCo6Xc7ys0u1QGs2SjDHNYEVuITc8+SnlnmoyZo3hawPsOXrm7CXFtgcI6HEZSzLGNLE31u3jpnkr6dopgn/88CKGJkS7HZIJUH3jAj/J2JiMMU1EVXnigx38aUk2o5O6kH5TKlHt2rodlglgvbu0Q4SAHvy3JGNME/BUVfOrNzaycNVeJl7Qkz9OGkp4qN0DY85NRNsQEjpHkms9GWNar5JyD3Ne/owPsguYfVk//uvKgbRpY/fAmKaRFNuBnYcD965/SzLGnINDx8qY/vxqth48zu++cz5TRttKlqZp9Y1tz5pdRahqQN7Aa0nGmEbadug4tz63mqOlFTw7LZXLBtpCY6bp9Y1rz4mKKvKPl9OtU4Tb4Zw1m11mTCN8kJ3PDU98SkVVNYtuH2sJxjSbU9OYA3Xw35KMMWdpwae7mP78anp1accbs8dxXnyU2yGZIHY6yQTouIxPSUZEJohItojkiMjddXweLiKLnM9Xikhijc/mOuXZIjK+oTZF5CWnfKOIzBeRtk65iMhfnPobRGT4uRy4MWersqqaX72+kXsyN3H54G688gNbB8Y0vx5RkYSFtGFPUanboTRKg0lGREKAx4GrgBTgRhFJqVVtBnBEVfsDjwAPONum4F2KeQgwAXhCREIaaPMlYBBwPhAJzHTKrwKSnZ9ZwJONOWBjGqP4ZCXTn1/Niyt2c/vX+/L01BG0D7chTdP8QtoICZ0j2VMYmEnGl2/JKCBHVXMBRCQDmAhsrlFnInCv8/pV4DHxToOYCGSoajmwU0RynPaor01VfftUoyKyCkiosY8X1Lte9AoRiRaRHqp64GwP2pizsbvwBNOfX82eolL+OGko30vt5XZIppXpHdMueHsyQDywt8b7PKeszjqq6gGKgZgzbNtgm85lspuAd88iDkRklohkiUhWQUFgP73UuG9lbiHffvwTCk9U8OKM0ZZgjCt6d2nHnsJSvH9jBxZfkkxdE7NrH2l9dc62vKYngA9V9aOziANVTVfVVFVNjYuzhxKaxluctZep81bSpX0Yb8wex5i+MW6HZFqp3l3acbzcw9HSwHsasy+Xy/KAmn++JQD766mTJyKhQBRQ1MC29bYpIvcAccDtZxmHMeesqlp54N2tpH+YyyXJsTw2ZThRkfYMMuOe3l3aAbCnqJTO7cNcjubs+NKTWQ0ki0iSiIThHcjPrFUnE5jmvJ4ELHPGTjKBNGf2WRLeQftVZ2pTRGYC44EbVbW61j5udmaZjQGKbTzGNLXi0kpufX416R/mcvPYPsy/ZaQlGOO6PjHeacy7A3BcpsGejKp6RGQOsAQIAear6iYRuR/IUtVMYB7wojOwX4Q3aeDUW4x3koAHmK2qVQB1tens8ilgN7DceYTC31X1fuBt4FtADlAK3NoUJ8CYU7YdOs6sF7LYd/Qkf7j+fNJG2SNijH/o1cU7VX5vMCYZAGfG19u1yn5d43UZ8N16tv0/4P98adMprzMmp2c025d4jTlb7206yJ2L1hEZFkrGrDGM6NPF7ZCMOa1dWCixHcLZXRh4d/3bRH/TqlVXK39dlsMjS7cxLCGKp24aQY8ou8HS+J8+ATqN2ZKMabVKyj381+J1LNl0iOuHx/O775xPRFtbA8b4p95d2rEyt9DtMM6aJRnTKu0uPMFtL2SRk1/Cr65JYfq4xIB8jLppPXp3acfr6/ZR7qkKqAXxLMmYVufDbQX8aOFaROCF6aO5ODnW7ZCMaVDvLu1QhX1HTtI3roPb4fjMnsJsWo3qauUv729n2nOr6N4pgszZF1uCMQGjT4z3XplAm8ZsPRnTKhSXVnLn4nUs25rPty/oye+uP592Yfa/vwkcp2/IDLAHZdq3zAS9TfuLueNvn3Gg+CT3TxzCTWP62PiLCThxHcMJD21D3hFLMsb4jVey9vLL1zfSuV0YGbPGMqJPZ7dDMqZRRIT46Ej2HT3pdihnxZKMCUpllVXc9+ZmFq7aw9i+Mfx1yoXEdgh3Oyxjzkl850j2HbEkY4yr8o6U8sOXPmNDXjE/+Ho/fv7NAYSG2BwXE/jioyPZcuCY22GcFUsyJqj8Z1sBP8lYS1WV8vRNIxg/pLvbIRnTZOKjIzlcUkFZZVXA3DhsScYEBU9VNQ//axtPfLCDgd068uTU4QF1L4Exvojv7H3k0b6jJ+kXIP9/W5IxAe9A8Ul+vHAtq3cdYXJqL+69bgiRYYHxV54xZyM+2kkyRyzJGNMi/p2dz88WraPcU82fJ1/Aty/8yorcxgSNmj2ZQGFJxgSkyqpqHnwvm6f/k8ug7h15/PvDA+YvO2Maq3unCELaSEDNMPNpyo2ITBCRbBHJEZG76/g8XEQWOZ+vFJHEGp/NdcqzRWR8Q22KyBynTEUktkb5pSJSLCLrnJ/T69mY1mX/0ZOkpa/g6f/kMmV0b16fPc4SjGkVQkPa0L1TRHD1ZEQkBHgcuBLIA1aLSKaqbq5RbQZwRFX7i0ga8AAwWURS8K6SOQToCSwVkQHONvW1+QnwFvBBHeF8pKrXNOI4TZB4f8sh/uuV9VR6qvnLjRdy3bCebodkTIuKj2EWKMgAABapSURBVA6se2V8uVw2CshR1VwAEckAJuJdUvmUicC9zutXgcfE+9yOiUCGqpYDO53lmUc59epsU1XXOmXnclwmyJR7qvjTu9k8+/FOUnp04vHvDycptr3bYRnT4uI7RwbUujK+XC6LB/bWeJ/nlNVZR1U9QDEQc4ZtfWmzLmNFZL2IvCMiQ+qqICKzRCRLRLIKCgp8aNL4u5z8Eq5/4lOe/XgnN43pw99/eJElGNNqJXSO5OCxMiqrqt0OxSe+9GTq6lKoj3XqK68rudVus7bPgD6qWiIi3wJeB5K/0ohqOpAOkJqa2lCbxo+pKotW7+W+NzcT0bYNz9ycypUp3dwOyxhXxUdHUq1wsLiMXs6Tmf2ZL0kmD+hV430CsL+eOnkiEgpEAUUNbNtQm1+iqsdqvH5bRJ4QkVhVPezDMZgAc7S0grtf+5x3Nx1kXP8YHv7eBXTrFOF2WMa4ruY05kBIMr5cLlsNJItIkoiE4R3Iz6xVJxOY5ryeBCxTVXXK05zZZ0l4ex6rfGzzS0SkuzPOg4iMcmIPnAuTxmfLdxQy4c8f8f7WQ8y9ahAvTh9tCcYYR80bMgNBgz0ZVfWIyBxgCRACzFfVTSJyP5ClqpnAPOBFZ2C/CG/SwKm3GO8kAQ8wW1WrwDtVuXabTvmPgbuA7sAGEXlbVWfiTV53iIgHOAmkOYnMBInKqmoe+dc2nvzPDpJi2vPMzeM4PyHK7bCM8Ss9owPrhkwJ5t/TqampmpWV5XYYxge7Dp/gJxlrWZ9XTNrIXvz62hRbudKYeqT+dimXD+rKA5OGNkv7IrJGVVOboi37FhtXnRrc/81bmwkNacOT3x/OVef3cDssY/xaz+gIDhwrczsMn1iSMa7JP17G3a99zrKt+VzUL4YHvzvs9KUAY0z9ekRFkFtwwu0wfGJJxrji7c8P8L//+JzSiiruuTaFaWMTadPGbsA1xhc9oiL5NCcw5j1ZkjEtqri0knsyN/L6uv0MTYji4e8No3/Xjm6HZUxA6R4VwfFyD8fLKukY0dbtcM7IkoxpMR9tL+C/X9lAQUk5P70imdmX9aetLYtszFnrEeWd0n+wuMySjDEnK6r4wztbWLB8N/3i2pN+80UMTYh2OyxjAlaPKO/Y5YHiMpK7+feVAEsypll9tucIP1+8ntzDJ5g+Lom7JgwMmLXJjfFXp3oyB4r9/14ZSzKmWZRVVvHQe9nM+3gn3TtF8PLM0VzUP7bhDY0xDerWKQIRb0/G31mSMU0ua1cRd726gdzDJ7hxVG/+51uD/P66sTGBJCy0DbEdwjloSca0JqUVHv60JJvnP91Fz6hIXpo5mnHWezGmWfSIimC/JRnTWqzILeQXr21gd2EpN4/tw10TBtEh3P73Mqa5dO8Uwa5C/78h034LmHNyotzDA+9u5YXlu+kT046MWWMY0zfG7bCMCXo9oyNZHgArZFqSMY32Sc5hfvHaBvYdPcn0cUn8fPwAe6ilMS2ke1QEx8s8lJR7/Pqqgf9GZvxWcWklv39nCxmr99I3tj2v3D6W1MQubodlTKvyxQ2ZJ/36qRmWZIzPVJW3Nhzgvjc3caS0ktu/3pc7rxhg970Y44KaN2T6c5Lx6ZkeIjJBRLJFJEdE7q7j83ARWeR8vlJEEmt8NtcpzxaR8Q21KSJznDIVkdga5SIif3E+2yAiwxt70Obs5R0pZfrzq/nRwrX0jI4kc8445l412BKMMS45fUPmUf+eYdZgT0ZEQoDHgSuBPGC1iGSq6uYa1WYAR1S1v4ikAQ8Ak0UkBe8qmUOAnsBSERngbFNfm58AbwEf1ArlKrzLNycDo4EnnX9NM/JUVfP8p7t46L1tiMCvr0lh2kWJhNgTk41x1aklyf39hkxfLpeNAnJUNRdARDKAiXiXVD5lInCv8/pV4DEREac8Q1XLgZ3O8syjnHp1tqmqa52y2nFMBF5wllxeISLRItJDVQ+czQEb323cV8zcv3/O5/uKuWxgHL/59nkkdG7ndljGGL64IdPfHy3jS5KJB/bWeJ/HV3sQp+uoqkdEioEYp3xFrW3jndcNtelLHPHAl5KMiMwCZgH07t27gSZNXU5WVPHnpdt49uOddG4XxmNTLuTq83vUlfiNMS7qGR0RFD2Zun6zqI916iuvayyodpuNiQNVTQfSAVJTUxtq09Ty7+x8fv3GRvYWnSRtZC/mXjWYqHb2SBhj/FH3ThHsLix1O4wz8iXJ5AG9arxPAPbXUydPREKBKKCogW0barMxcZhG2n/0JPe/uZl3Nx2kb1x7u6nSmAAQCDdk+pJkVgPJIpIE7MM7kD+lVp1MYBqwHJgELFNVFZFM4GUReRjvwH8ysApvr6ShNmvLBOY44zejgWIbjzl3FZ5q5n+yk0eXbkdR/nv8QGZekkR4qM0aM8bfBcINmQ1G5YyxzAGWACHAfFXdJCL3A1mqmgnMA150BvaL8CYNnHqL8U4S8ACzVbUKvFOVa7fplP8YuAvoDmwQkbdVdSbwNvAtIAcoBW5tqpPQWi3fUciv3thITn4JVwzuxj3XptCriw3sGxMoaq6Q2b9rB5ejqZt4J2sFp9TUVM3KynI7DL+Tf7yM3/1zC6+v209C50juvXYIV6R0czssY8xZWrWziO89vZy/zRjNxclN98RzEVmjqqlN0ZZ/9q9Ms6iqVv62YjcPLsmm3FPNj77Rnx9e2p/IMLs0ZkwgOtWT2e/H05gtybQSn+05wq9e38im/ce4JDmW+64bQt84/+xeG2N807VTOIBfL15mSSbIHTpWxh/e2co/1u6jW6dwHp8ynG+d393ueTEmCISHhhDbIcyv75WxJBOkyj1VzPt4J48ty8FTpcy+rB8/vLQ/7f10BooxpnG6R0Vw0C6XmZaiqizdks9v/7mZ3YWlXJnSjV9ePZg+Me3dDs0Y0wy6d4ok74j/3pBpSSaI5OQf5743N/PR9sP079qBF2eM4pLkOLfDMsY0o57REazeVeR2GPWyJBMEik9W8uel23hh+W7ahYXw62tSuGlsH9qG+LSSgzEmgHWPiqD4ZCWlFR6/XJnW/yIyPvNUVbMoay8PvbeNI6UVpI3szc+/OYCYDuFuh2aMaSE1b8j0xxmjlmQCkKryQXYBv3t7C9vzSxiV2IVfX5vCefFRbodmjGlh3Tt5V8i0JGOaxKb9xfzu7S18klNIYkw7npo6gvFDutmUZGNaqdMrZPrpNGZLMgHiYHEZD76XzWuf5REV2ZZ7rk3h+6P7EBZq4y7GtGbdTycZ/5zGbEnGz50o9/D0f3aQ/lEu1dVw2yV9mX1pf1vjxRgDQETbELq0998bMi3J+KmqamWxM6h/uKSca4b24K7xg+gdY09JNsZ8WfdOEX77aBlLMn6mulp5d9NBHnovmx0FJxjRpzPpN49geO/ObodmjPFTPaL8dxlmSzJ+QlX5z7YCHnwvm437jtG/aweemjqc8UPsOWPGmDPrHhXB2r1H3Q6jTj6NGovIBBHJFpEcEbm7js/DRWSR8/lKEUms8dlcpzxbRMY31KaIJDltbHfaDHPKbxGRAhFZ5/zMPJcD9yerdxUx+ekV3PLcao6WVvLQd4ex5KdfY8J5PSzBGGMa1CMqgqITFZRVVrkdylc02JMRkRDgceBKIA9YLSKZqrq5RrUZwBFV7S8iacADwGQRScG7SuYQvMsvLxWRAc429bX5APCIqmaIyFNO20862yxS1TnneMx+Y+O+Yh56L5t/ZxcQ1zGc30wcwuSRvW3GmDHmrPSI+uJemcRY/3pOoS+Xy0YBOaqaCyAiGcBEvEsqnzIRuNd5/SrwmHj/BJ8IZKhqObDTWZ55lFPvK22KyBbgG8AUp84Cp91TSSYo7Cgo4eF/beOfGw4QFdmWX0wYxC0XJdriYcaYRql5r0wgJpl4YG+N93nA6PrqqKpHRIqBGKd8Ra1t453XdbUZAxxVVU8d9QFuEJGvAduAO1W1Zht+b9/Rkzy6dBuvrskjom0IP/pGf2Ze0peoSJuObIxpvFP3yhw85n/3yviSZOoaFFAf69RXXtf1oDPVB3gTWKiq5SLyA7y9nG98JViRWcAsgN69e9fRXMsrOF7OEx/k8NKKPQDcclESP7ysH7H2jDFjTBPo7sd3/fuSZPKAXjXeJwD766mTJyKhQBRQ1MC2dZUfBqJFJNTpzZyur6qFNeo/g3fs5itUNR1IB0hNTa2dDFtU8clKnvkwl/mf7KTcU82k4Qn8+Ipk4qMj3QzLGBNk2oWFEhXZ1i/vlfElyawGkkUkCdiHdyB/Sq06mcA0YDkwCVimqioimcDLIvIw3oH/ZGAV3h7LV9p0tvm300aG0+YbACLSQ1UPOPu7DtjSyGNudqUVHp7/dBdPfbCDY2UerhnagzuvHEA/P3x4nTEmOPSIimD/0QBMMs4YyxxgCRACzFfVTSJyP5ClqpnAPOBFZ2C/CG/SwKm3GO8kAQ8wW1WrAOpq09nlL4AMEfktsNZpG+DHInKd004RcMs5H30Tq/BUs2j1Hh59P4fDJeV8Y1BX/uubAxjS056ObIxpXj2iIvxyTEZUXb2i1KxSU1M1Kyur2fdTXa28uWE/D723jT1FpYxK7MJdEwaSmtil2fdtjDEAc//+Of/afJCsX155zm2JyBpVTW2CsOyO/3P14bYCfv/OVrYcOMag7h157paRXDowzm6iNMa0qB5RERwuqaDcU0V4qP/cDmFJppH2FpXym7c2897mQ/TqEsmfJ1/AdcN60qaNJRdjTMs7NcMs/1g5vbr4z4N0LcmcpepqZf4nO/nTkmzaiPDf4wcy85Ikv/rLwRjT+py6IXP/0ZOWZAJVwfFy7ly0jo9zDnPF4K7cN/E8m45sjPELPU7fkOlfM8wsyfgoJ7+EW55bxeGScn5//fmkjexl4y7GGL/R3Xl+mb/dkGlJxge7C0+Qlr4cEBbfPpahCdFuh2SMMV/SITyUjhGhfndDpiWZBpSUe5g2fxVV1corPxhL/652Q6Uxxj95Fy/zr3tlLMk04P/+uYXdRaUsmmUJxhjj37pHRfrd5TJbuOQMth48xsJVe5gxLolRSXZjpTHGv/Xo5H/LMFuSOYPH/72DjuGhzPlGf7dDMcaYBnWPiuBwSTkVnmq3QznNkkw9ik9WsmTTQW4YkUB0uzC3wzHGmAb1jI5AFfKP+09vxpJMPZZuPkSFp5pvXxjfcGVjjPED/jiN2ZJMPVbuLCQqsi1D4+0JysaYwJDQ2ZtkdheWuhzJFyzJ1CNr9xFGJna2Z5EZYwJG7y7tCG0j5BaUuB3KaZZk6uCpqmZPYSkDunV0OxRjjPFZ25A29O7SjtyCE26HcpolmTrsP1qGp1rpE+M/D5kzxhhf9I3rQO7hAOvJiMgEEckWkRwRubuOz8NFZJHz+UoRSazx2VynPFtExjfUpogkOW1sd9oMa2gfTa2gxDtodmoQzRhjAkW/uPbsKiylqto/FqRsMMmISAjwOHAVkALcKCIptarNAI6oan/gEeABZ9sUvEsxDwEmAE+ISEgDbT4APKKqycARp+1699EcTpRXAdA+zB7fb4wJLMndOlLhqSYn3z96M770ZEYBOaqaq6oVQAYwsVadicAC5/WrwOXifUTxRCBDVctVdSeQ47RXZ5vONt9w2sBp89sN7KPJlVZ4k0ykJRljTIAZ2y8GgHc3HnQ5Ei9fkkw8sLfG+zynrM46quoBioGYM2xbX3kMcNRpo/a+6tvHl4jILBHJEpGsgoICHw7vq+I6hnHVed2JaR/eqO2NMcYt8dGRXH9hPCk9O7kdCuDbAzLr6i3UvthXX536yutKbmeq72scqGo6kA6QmpraqIuSI/p0YUQfe1aZMSYwPTz5ArdDOM2Xnkwe0KvG+wRgf311RCQUiAKKzrBtfeWHgWinjdr7qm8fxhhj/JQvSWY1kOzM+grDO5CfWatOJjDNeT0JWKaq6pSnOTPDkoBkYFV9bTrb/NtpA6fNNxrYhzHGGD/V4OUyVfWIyBxgCRACzFfVTSJyP5ClqpnAPOBFEcnB27tIc7bdJCKLgc2AB5itqlUAdbXp7PIXQIaI/BZY67RNffswxhjjvySYOwOpqamalZXldhjGGBNQRGSNqqY2RVt2x78xxphmY0nGGGNMs7EkY4wxptlYkjHGGNNsgnrgX0QKgBN4778xEIudi1PsXHjZefiCnQuvWKC9qsY1RWNBnWQARCSrqWZJBDo7F1+wc+Fl5+ELdi68mvo82OUyY4wxzcaSjDHGmGbTGpJMutsB+BE7F1+wc+Fl5+ELdi68mvQ8BP2YjDHGGPe0hp6MMcYYl1iSMcYY02yCOsmIyAQRyRaRHBG52+14mpqIzBeRfBHZWKOsi4j8S0S2O/92dspFRP7inIsNIjK8xjbTnPrbRWRaXfvydyLSS0T+LSJbRGSTiPzEKW9V50NEIkRklYisd87DfU55koisdI5pkbPEBs4yHIuc87BSRBJrtDXXKc8WkfHuHNG5E5EQEVkrIm8571vluRCRXSLyuYisE5Esp6z5vx+qGpQ/eJcQ2AH0BcKA9UCK23E18TF+DRgObKxR9kfgbuf13cADzutvAe/gXWF0DLDSKe8C5Dr/dnZed3b72BpxLnoAw53XHYFtQEprOx/O8XRwXrcFVjrHtxhIc8qfAu5wXv8QeMp5nQYscl6nON+ZcCDJ+S6FuH18jTwnPwNeBt5y3rfKcwHsAmJrlTX79yOYezKjgBxVzVXVCiADmOhyTE1KVT/kq6uDTgQWOK8XAN+uUf6Ceq3AuwJpD2A88C9VLVLVI8C/gAnNH33TUtUDqvqZ8/o4sAWIp5WdD+d4Spy3bZ0fBb4BvOqU1z4Pp87Pq8DlIiJOeYaqlqvqTiAH73cqoIhIAnA18KzzXmil56Iezf79COYkEw/srfE+zykLdt1U9QB4f/ECXZ3y+s5H0J0n5zLHhXj/im9158O5PLQOyMf7S2AHcFRVPU6Vmsd0+nidz4uBGILgPDj+DNwFVDvvY2i950KB90RkjYjMcsqa/fvR4MqYAUzqKGvN87XrOx9BdZ5EpAPwGvBTVT3m/UO07qp1lAXF+VDv6rMXiEg08A9gcF3VnH+D9jyIyDVAvqquEZFLTxXXUTXoz4VjnKruF5GuwL9EZOsZ6jbZuQjmnkwe0KvG+wRgv0uxtKRDTrcW5998p7y+8xE050lE2uJNMC+p6t+d4lZ7PlT1KPAB3mvq0SJy6o/Kmsd0+nidz6PwXoINhvMwDrhORHbhvVz+Dbw9m9Z4LlDV/c6/+Xj/+BhFC3w/gjnJrAaSnZkkYXgH8jJdjqklZAKnZnxMA96oUX6zM2tkDFDsdI+XAN8Ukc7OzJJvOmUBxbl2Pg/YoqoP1/ioVZ0PEYlzejCISCRwBd7xqX8Dk5xqtc/DqfMzCVim3hHeTCDNmXGVBCQDq1rmKJqGqs5V1QRVTcT7/V+mqt+nFZ4LEWkvIh1Pvcb7//VGWuL74faMh+b8wTtDYhvea9L/63Y8zXB8C4EDQCXevzBm4L2G/D6w3fm3i1NXgMedc/E5kFqjnel4BzNzgFvdPq5GnouL8XbbNwDrnJ9vtbbzAQwF1jrnYSPwa6e8L95fjDnAK0C4Ux7hvM9xPu9bo63/dc5PNnCV28d2juflUr6YXdbqzoVzzOudn02nfh+2xPfDHitjjDGm2QTz5TJjjDEusyRjjDGm2ViSMcYY02wsyRhjjGk2lmSMMcY0G0syxhhjmo0lGWOMMc3m/wFq2ptM7JrRDQAAAABJRU5ErkJggg==\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": 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": 15,
"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": null,
"metadata": {},
"outputs": [],
"source": [
"total_f.update_integration_options(draws_per_dim=200000, mc_sampler=None)\n",
"# inte = total_f.integrate(limits = (2000, x_max), 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() - _))))\n",
"\n",
"# print(pdg['NR_BR']/pdg['NR_auc']*inte_fl)\n",
"# print(0.25**2*4.2/1000)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Sampling\n",
"## 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.04, 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(1000, 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(1040, 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",
"# for param in total_f.get_dependents():\n",
"# param.randomize()\n",
" \n",
"# # for param in total_f.get_dependents():\n",
"# # print(zfit.run(param))\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)\n",
"# # print(\"Results:\", result.params)\n",
"# print(\"Hesse errors:\", result.hesse())"
]
},
{
"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, 10e-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": [
"# def jpsi_res(q):\n",
"# return resonance(q, _mass = jpsi_mass, scale = jpsi_scale,\n",
"# phase = jpsi_phase, width = jpsi_width)\n",
"\n",
"# def psi2s_res(q):\n",
"# return resonance(q, _mass = psi2s_mass, scale = psi2s_scale,\n",
"# phase = psi2s_phase, width = psi2s_width)\n",
" \n",
"# def p3770_res(q):\n",
"# return resonance(q, _mass = p3770_mass, scale = p3770_scale,\n",
"# phase = p3770_phase, width = p3770_width)\n",
" \n",
"# def p4040_res(q):\n",
"# return resonance(q, _mass = p4040_mass, scale = p4040_scale,\n",
"# phase = p4040_phase, width = p4040_width)\n",
" \n",
"# def p4160_res(q):\n",
"# return resonance(q, _mass = p4160_mass, scale = p4160_scale,\n",
"# phase = p4160_phase, width = p4160_width)\n",
" \n",
"# def p4415_res(q):\n",
"# return resonance(q, _mass = p4415_mass, scale = p4415_scale,\n",
"# phase = p4415_phase, width = p4415_width)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# 0.15**2*4.2/1000\n",
"# result.hesse()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Analysis"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"# zfit.run.numeric_checks = False \n",
"\n",
"Ctt_list = []\n",
"Ctt_error_list = []\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",
" ### Generate data\n",
" \n",
" print(\"Toy {}: Generating data...\".format(toy))\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",
" with open(\"data/zfit_toys/toy_{0}/{1}.pkl\".format(toy, call), \"wb\") as f:\n",
" pkl.dump(sam, f, pkl.HIGHEST_PROTOCOL)\n",
" \n",
" print(\"Toy {}: Data generation finished\".format(toy))\n",
" \n",
" ### Load data\n",
" \n",
" print(\"Toy {}: Loading data...\".format(toy))\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",
" data = zfit.data.Data.from_numpy(array=total_samp[:int(nevents)], obs=obs)\n",
" \n",
" print(\"Toy {}: Loading data finished\".format(toy))\n",
"\n",
" ### Fit data\n",
" \n",
" print(\"Toy {}: Fitting pdf...\".format(toy))\n",
"\n",
" for param in total_f.get_dependents():\n",
" param.randomize()\n",
"\n",
" nll = zfit.loss.UnbinnedNLL(model=total_f, data=data, fit_range = (jpsi_mass+50.0, psi2s_mass-50.0))\n",
"\n",
" minimizer = zfit.minimize.MinuitMinimizer(verbosity = 5)\n",
" # minimizer._use_tfgrad = False\n",
" result = minimizer.minimize(nll)\n",
" \n",
" print(\"Toy {}: Fitting finished\".format(toy))\n",
"\n",
" print(\"Function minimum:\", result.fmin)\n",
" print(\"Hesse errors:\", result.hesse())\n",
" \n",
" params = result.params\n",
" Ctt_list.append(params[Ctt]['value'])\n",
" Ctt_error_list.append(params[Ctt]['minuit_hesse']['error'])\n",
"\n",
" print(\"Toy {0}/{1}\".format(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"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print('Mean Ctt value = {}'.format(np.mean(Ctt_list)))\n",
"print('Mean Ctt error = {}'.format(np.mean(Ctt_error_list)))"
]
},
{
"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
}
