{
"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": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"def formfactor(q2, subscript, b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2): #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",
"\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",
" b0 = [b0_0, b0_1, b0_2]\n",
"\n",
" for i in range(N):\n",
" _sum += b0[i]*(tf.pow(z,i))\n",
"\n",
" return ztf.to_complex(prefactor * _sum)\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",
" bT = [bT_0, bT_1, bT_2]\n",
" for i in range(N):\n",
" _sum += bT[i] * (tf.pow(z, i) - ((-1)**(i-N)) * (i/N) * tf.pow(z, N))\n",
" else:\n",
" bplus = [bplus_0, bplus_1, bplus_2]\n",
" for i in range(N):\n",
" _sum += bplus[i] * (tf.pow(z, i) - ((-1)**(i-N)) * (i/N) * tf.pow(z, N))\n",
"\n",
" return ztf.to_complex(prefactor * _sum)\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",
"\n",
"def axiv_nonres(q, b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_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",
" 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 = 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. * tf.pow(kabs,2) * tf.pow(beta,2) * tf.pow(tf.abs(ztf.to_complex(C10eff)*formfactor(q2, \"+\", b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2)),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(ztf.to_complex(C10eff) * formfactor(q2, \"0\", b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2)), 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 * q\n",
"\n",
"def vec(q, funcs, b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2):\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 = 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 = tf.pow(kabs,2) * (1. - 1./3. * beta)\n",
"\n",
" abs_bracket = tf.pow(tf.abs(c9eff(q, funcs) * formfactor(q2, \"+\", b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2) + ztf.to_complex(2.0 * C7eff * (mb + ms)/(mB + mK)) * formfactor(q2, \"T\", b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2)),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 * q\n",
"\n",
"def c9eff(q, funcs):\n",
"\n",
" C9eff_nr = ztf.to_complex(ztf.constant(pdg['C9eff']))\n",
"\n",
" c9 = C9eff_nr + 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": [
"## Build pdf"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"class total_pdf(zfit.pdf.ZPDF):\n",
" _N_OBS = 1 # dimension, can be omitted\n",
" _PARAMS = ['b0_0', 'b0_1', 'b0_2', \n",
" 'bplus_0', 'bplus_1', 'bplus_2', \n",
" 'bT_0', 'bT_1', 'bT_2', \n",
" '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",
" 'Dstar_mass', 'DDstar_scale', 'DDstar_phase', 'D_mass',\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",
" b0 = [self.params['b0_0'], self.params['b0_1'], self.params['b0_2']]\n",
" bplus = [self.params['bplus_0'], self.params['bplus_1'], self.params['bplus_2']]\n",
" bT = [self.params['bT_0'], self.params['bT_1'], self.params['bT_2']]\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 ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['jpsi_mass'], 2)) * resonance(q, _mass = self.params['jpsi_mass'], \n",
" scale = self.params['jpsi_scale'],\n",
" phase = self.params['jpsi_phase'], \n",
" width = self.params['jpsi_width'])\n",
" def psi2s_res(q):\n",
" return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['psi2s_mass'], 2)) * resonance(q, _mass = self.params['psi2s_mass'], \n",
" scale = self.params['psi2s_scale'],\n",
" phase = self.params['psi2s_phase'], \n",
" width = self.params['psi2s_width'])\n",
" def p3770_res(q):\n",
" return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['p3770_mass'], 2)) * resonance(q, _mass = self.params['p3770_mass'], \n",
" scale = self.params['p3770_scale'],\n",
" phase = self.params['p3770_phase'], \n",
" width = self.params['p3770_width'])\n",
" \n",
" def p4040_res(q):\n",
" return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['p4040_mass'], 2)) * resonance(q, _mass = self.params['p4040_mass'], \n",
" scale = self.params['p4040_scale'],\n",
" phase = self.params['p4040_phase'], \n",
" width = self.params['p4040_width'])\n",
" \n",
" def p4160_res(q):\n",
" return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['p4160_mass'], 2)) * resonance(q, _mass = self.params['p4160_mass'], \n",
" scale = self.params['p4160_scale'],\n",
" phase = self.params['p4160_phase'], \n",
" width = self.params['p4160_width'])\n",
" \n",
" def p4415_res(q):\n",
" return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['p4415_mass'], 2)) * resonance(q, _mass = self.params['p4415_mass'], \n",
" scale = self.params['p4415_scale'],\n",
" phase = self.params['p4415_phase'], \n",
" 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['Dstar_mass'], q)) + ztf.to_complex(h_P(self.params['D_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, funcs, b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2)\n",
"\n",
" axiv_nr = axiv_nonres(x, b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2)\n",
"\n",
" tot = vec_f + axiv_nr\n",
"\n",
" return tot"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup parameters"
]
},
{
"cell_type": "code",
"execution_count": 6,
"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": [
"# formfactors\n",
"\n",
"b0_0 = zfit.Parameter(\"b0_0\", ztf.constant(0.292), floating = False) #, lower_limit = -2.0, upper_limit= 2.0)\n",
"b0_1 = zfit.Parameter(\"b0_1\", ztf.constant(0.281), floating = False) #, lower_limit = -2.0, upper_limit= 2.0)\n",
"b0_2 = zfit.Parameter(\"b0_2\", ztf.constant(0.150), floating = False) #, lower_limit = -2.0, upper_limit= 2.0)\n",
"\n",
"bplus_0 = zfit.Parameter(\"bplus_0\", ztf.constant(0.466), lower_limit = -2.0, upper_limit= 2.0)\n",
"bplus_1 = zfit.Parameter(\"bplus_1\", ztf.constant(-0.885), lower_limit = -2.0, upper_limit= 2.0)\n",
"bplus_2 = zfit.Parameter(\"bplus_2\", ztf.constant(-0.213), lower_limit = -2.0, upper_limit= 2.0)\n",
"\n",
"bT_0 = zfit.Parameter(\"bT_0\", ztf.constant(0.460), floating = False) #, lower_limit = -2.0, upper_limit= 2.0)\n",
"bT_1 = zfit.Parameter(\"bT_1\", ztf.constant(-1.089), floating = False) #, lower_limit = -2.0, upper_limit= 2.0)\n",
"bT_2 = zfit.Parameter(\"bT_2\", ztf.constant(-1.114), floating = False) #, lower_limit = -2.0, upper_limit= 2.0)\n",
"\n",
"\n",
"#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, 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",
"\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), floating = False) #, 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), floating = False) #, 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), lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"p3770_s = zfit.Parameter(\"p3770_s\", ztf.constant(p3770_scale), 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), lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"p4040_s = zfit.Parameter(\"p4040_s\", ztf.constant(p4040_scale), 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), lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"p4160_s = zfit.Parameter(\"p4160_s\", ztf.constant(p4160_scale), 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), lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"p4415_s = zfit.Parameter(\"p4415_s\", ztf.constant(p4415_scale), 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": 7,
"metadata": {},
"outputs": [],
"source": [
"m_c = 1300\n",
"\n",
"Dbar_phase = 0.0\n",
"DDstar_phase = 0.0\n",
"Dstar_mass = pdg['Dst_M']\n",
"Dbar_mass = pdg['D0_M']\n",
"D_mass = pdg['D0_M']\n",
"\n",
"Dbar_s = zfit.Parameter(\"Dbar_s\", ztf.constant(0.0), lower_limit=-0.3, upper_limit=0.3)\n",
"Dbar_m = zfit.Parameter(\"Dbar_m\", ztf.constant(Dbar_mass), floating = False)\n",
"Dbar_p = zfit.Parameter(\"Dbar_p\", ztf.constant(Dbar_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)#, floating = False)\n",
"DDstar_s = zfit.Parameter(\"DDstar_s\", ztf.constant(0.0), lower_limit=-0.3, upper_limit=0.3)#, floating = False)\n",
"Dstar_m = zfit.Parameter(\"Dstar_m\", ztf.constant(Dstar_mass), floating = False)\n",
"D_m = zfit.Parameter(\"D_m\", ztf.constant(D_mass), floating = False)\n",
"DDstar_p = zfit.Parameter(\"DDstar_p\", ztf.constant(DDstar_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)#, floating = False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Tau parameters"
]
},
{
"cell_type": "code",
"execution_count": 8,
"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": [
"## Load data"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"x_min = 2*pdg['muon_M']\n",
"x_max = (pdg[\"Bplus_M\"]-pdg[\"Ks_M\"]-0.1)\n",
"\n",
"# # Full spectrum\n",
"\n",
"# obs = zfit.Space('q', limits = (x_min, x_max))\n",
"\n",
"# Jpsi and Psi2s cut out\n",
"\n",
"obs1 = zfit.Space('q', limits = (x_min, jpsi_mass - 50.))\n",
"obs2 = zfit.Space('q', limits = (jpsi_mass + 50., psi2s_mass - 50.))\n",
"obs3 = zfit.Space('q', limits = (psi2s_mass + 50., x_max))\n",
"\n",
"obs = obs1 + obs2 + obs3\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 pdf"
]
},
{
"cell_type": "code",
"execution_count": 10,
"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",
" Dstar_mass = Dstar_m, DDstar_scale = DDstar_s, DDstar_phase = DDstar_p, D_mass = D_m,\n",
" Dbar_mass = Dbar_m, Dbar_scale = Dbar_s, Dbar_phase = Dbar_p,\n",
" tau_mass = tau_m, C_tt = Ctt, b0_0 = b0_0, b0_1 = b0_1, b0_2 = b0_2,\n",
" bplus_0 = bplus_0, bplus_1 = bplus_1, bplus_2 = bplus_2,\n",
" bT_0 = bT_0, bT_1 = bT_1, bT_2 = bT_2)\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": 11,
"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, int(2e6))\n",
"\n",
"probs = total_f.pdf(test_q, norm_range=False)\n",
"\n",
"calcs_test = zfit.run(probs)\n",
"res_y = zfit.run(jpsi_res(test_q))\n",
"# b0 = [b0_0, b0_1, b0_2]\n",
"# bplus = [bplus_0, bplus_1, bplus_2]\n",
"# bT = [bT_0, bT_1, bT_2]\n",
"# f0_y = zfit.run(tf.math.real(formfactor(test_q,\"0\", b0, bplus, bT)))\n",
"# fplus_y = zfit.run(tf.math.real(formfactor(test_q,\"+\", b0, bplus, bT)))\n",
"# fT_y = zfit.run(tf.math.real(formfactor(test_q,\"T\", b0, bplus, bT)))"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\ipykernel_launcher.py:12: UserWarning: Creating legend with loc=\"best\" can be slow with large amounts of data.\n",
" if sys.path[0] == '':\n",
"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\IPython\\core\\pylabtools.py:128: UserWarning: Creating legend with loc=\"best\" can be slow with large amounts of data.\n",
" fig.canvas.print_figure(bytes_io, **kw)\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAaUAAAD4CAYAAABMtfkzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO29eXyc1XX//z4zo5Es2ZKszZtsJNvCRgazOWYxAWMImJLEfFtoTEpCElraNLRN0gX4pkm+JeHX0KahSbNSSELSBENIWlxCMAQbEjYbs+NFtrxgCy9aLVn7LOf3x/OMNB7PaGZkSTOSzvv18kvP3Ln33DuPpecz595zzxVVxTAMwzCyAU+mB2AYhmEYEUyUDMMwjKzBRMkwDMPIGkyUDMMwjKzBRMkwDMPIGnyZHkC2UVZWplVVVZkehmFkPbuOHicvx8u8kvwTytt7Ahxo7aamYhp5OZn93ru7sRO/18NppfkEw8qOwx3MLp5CaYE/o+OaiLz66qvNqlp+qnZMlGKoqqpi69atmR6GYWQ9l/3rJs6ZW8w31557Qvlv3j7Mp3/2Gr/4m/dzxqzCDI3OYfW//465Jfn858eX0dbVz7lfeZovfaiWT66ozui4JiIi8u5I2LHpO8MwhkVYFY/ISeUej1MWCmfHHsjICLNtXEZ8TJQMwxgW4TBxRcnnPvzDWbAxXxUiQ/Rm0biMxJgoGYYxLBxP6eTyiEcSzAKPRFHE9ZW8EvGUMjkiIxm2pmQYxrAIhXXA+4gm8vAPZ4MoRXlKHvcr+Fh7SoFAgIaGBnp7e8e039EiLy+PyspKcnJyRsW+iZJhGMMirCBDTN9lw9qNEjV9J5kZV0NDA9OmTaOqqiru/RpPqCotLS00NDRQXT06wSI2fWcYxrAIq+KN8wTJpoACVR0QAm+GxtXb20tpaem4FyRwvoSUlpaOqtdnomQYxrBIFH038PDPgoAC1cHoOxFBJDOBDhNBkCKM9mcxUTIMY1iEwklEKRs8JU58iHpFsmJcRmJMlAzDGBaq8UPCM7V2Ew9VJXqEHo9khQeXjTz77LN88IMfBKCvr48rr7ySc845h4cffnhMx2GBDoZhDAsn+u7k8uzzlAZfe0WyIiow23n99dcJBAK88cYbY953Sp6SiKwWkToRqReRO+K8nysiD7vvbxaRqqj37nTL60Tk6mQ2RaTatbHbtekfqg8RKRWRTSLSKSLfTjD+9SLyTmq3xDCMVEiY0SESEp4FHkn0mhI4gjkZ9ynt37+fxYsXc/PNN7N06VKuv/56uru7efLJJ1m8eDGXXHIJv/rVrwBobGzkpptu4o033uCcc85hz549YzrWpJ6SiHiB7wAfABqAV0Rkvapuj6p2C9CmqgtFZC1wD/AREakF1gJLgNnAb0XkdLdNIpv3APeq6joR+b5r+3uJ+gB6gS8CZ7r/Ysf/h0BnWnfFMIykhFUHIu2i8XmzbPNslHB6MhToEOGf/ncb2w91jKjN2tmFfPlDS5LWq6ur44EHHmDFihV86lOf4hvf+AY/+MEP2LhxIwsXLuQjH/kIABUVFdx///18/etf5/HHHx/RsaZCKp7ScqBeVfeqaj+wDlgTU2cN8KB7/ShwhTi/CWuAdarap6r7gHrXXlybbptVrg1cm9cN1Yeqdqnq8zjidAIiMhX4PPDVFD6nYRhpEFbiZ3TIqjWleJ5S5seVCebOncuKFSsAuOmmm9i6dSvV1dXU1NQgItx0000ZHqFDKmtKc4CDUa8bgAsS1VHVoIi0A6Vu+csxbee41/FslgLHVDUYp36iPpqHGPtXgH8Duof6gCJyK3ArwLx584aqahiGSyisA0EN0WRTjjmN2eDrzXCgQyoezWgRG8rd3t6elaHqqXhK8UYd+7+aqM5Ilac6jsEBiZwDLFTV/05UZ8CI6n2qukxVl5WXn/JxIIYx4VH3wT50RocxHVJcnM2zg6+9HiEUyrxYZoIDBw7w0ksvAfDQQw9x5ZVXsm/fvoE1o4ceeiiTwxsgFVFqAOZGva4EDiWqIyI+oAhoHaJtovJmoNi1EdtXoj4ScRFwvojsB54HTheRZ4f8pIZhpERkvcgXZ/5uMKND5lUp9puuVyZvSPgZZ5zBgw8+yNKlS2ltbeVzn/sc9913H9deey2XXHIJp512WqaHCKQ2ffcKUCMi1cB7OIELH42psx64GXgJuB7YqKoqIuuBn4vIN3ACHWqALTi/JyfZdNtscm2sc20+NlQfiQatqt/DCZDAjdR7XFVXpvB5DcNIQmRdJl6gQzZl445OyArOeCdrSLjH4+H73//+CWWrV69m586dJ9VduXIlK1euHKORnUhSUXLXb24DNgBe4Iequk1E7gK2qup64AHgpyJSj+O9rHXbbhORR4DtQBD4jKqGAOLZdLu8HVgnIl8FXndtk6gP19Z+oBDwi8h1wFUx0YGGYYwgkfWieJ5SVqUZijq6AjK/pmQkJ6XNs6r6BPBETNmXoq57gRsStL0buDsVm275XpzovNjyofqoSjL+/cQJFzcMY3hEpu/iHl0REaUscJViPaXJmmaoqqqKd94ZH1s1Lc2QYRhpE5kCGzLNUBY8+2MzOng8kpGowCFWGsYdo/1ZTJQMw0ibiLcR2SgbjdebXYf8RYc6ZMJTysvLo6WlZUIIU+Q8pby8vFHrw3LfGYaRNqEUPKVsyOgAJx7Z7slAmqHKykoaGhpoamoa245HicjJs6OFiZJhGGkTCRaIt6aUqWPH4xGOXVPyjP24cnJyRu2U1omITd8ZhpE2wdAQgQ5ZlWYoJvpukgY6jCdMlAzDSJuItzFUmqFsmL7LlkAHI3VMlAzDSJuhAh1ExMnGnQ2iFJuQ1TylrMdEyTCMtBkq0AGyZ5Oqk/su6uiKSZwlfLxgomQYRtoMFegQKc8KTynmtVds+i7bMVEyDCNtQkNkdADn4Z8Na0roid7cZD5PabxgomQYRtoMiFKC6btsmSYLxxxd4fFIVmSaMBJjomQYRtoMiFKcQAdwErVmwzRZ7NEVPo8MHKnR0Rvg+d1DnRFqZAITJcMw0iaZp+T1ZMf03UkJWT0ysMfq9kff4qYHNtN0vC9DozPiYaJkGEbaJFtT8ki2BDqcGH2X4x0Uy9cPHAOgvSeQkbEZ8TFRMgwjbZJF3/myyVOKeu3zeAi6ye/ycpzH3/FeE6VswkTJMIy0SRp9582OQIeYJOH4vEIgdOKpuX3BzJ/7ZAxiomQYRtok2zyb4/EQyIJD/lBOyH2X4/EQdAMdIqX9JkpZhYmSYRhpM9Rx6AA53uwQJY05uiLHNxjoEFlrMlHKLkyUDMNIm6GyhIMzTRbMgg1BsUdX+KI8uMjQ+7NAPI1BTJQMw0ibcLJAB6+HQDasKcUcXZETtaYUKTdPKbtISZREZLWI1IlIvYjcEef9XBF52H1/s4hURb13p1teJyJXJ7MpItWujd2uTf9QfYhIqYhsEpFOEfl2lJ18Efm1iOwUkW0i8rX0b49hGPEIJgl0yPHIQJRbJok9usLnjVpTinhKJkpZRVJREhEv8B3gGqAWuFFEamOq3QK0qepC4F7gHrdtLbAWWAKsBr4rIt4kNu8B7lXVGqDNtZ2wD6AX+CLwd3GG/3VVXQycC6wQkWuSfV7DMJKTLNAhW6bvYkPCczyOpxSdPbwvC8TTGCQVT2k5UK+qe1W1H1gHrImpswZ40L1+FLhCnP/xNcA6Ve1T1X1AvWsvrk23zSrXBq7N64bqQ1W7VPV5HHEaQFW7VXWTe90PvAaM3sHyhjGJSCnQIZwlD/so4fR5nUdeKKwWfZelpCJKc4CDUa8b3LK4dVQ1CLQDpUO0TVReChxzbcT2laiPpIhIMfAh4JkE798qIltFZGtTU1MqJg1jUpM00MGTeU9JNbJ2NEiOK0rBsOJxn359wdAYj8wYilREKd5vXexvW6I6I1We6jhOQkR8wEPAt1R1b7w6qnqfqi5T1WXl5eXJTBrGpCfiKXmGCnTI8LRYJB9s9AxjjptANhAKDwQ6BIKZn2Y0BklFlBqAuVGvK4FDieq4IlAEtA7RNlF5M1Ds2ojtK1EfybgP2K2q/55CXcMwUiCiN4mm7/zZIEruz+h1r8h4gyEdWBfrD5mnlE2kIkqvADVuVJwfJ3BhfUyd9cDN7vX1wEZ1fOf1wFo3cq4aqAG2JLLpttnk2sC1+ViSPhIiIl/FEa/PpvA5DcNIkcjxD0MGOmQ4JDwcZ/ousqYUCIUHRNPWlLILX7IKqhoUkduADYAX+KGqbhORu4CtqroeeAD4qYjU43gva92220TkEWA7EAQ+o6ohgHg23S5vB9a5gvK6a5tEfbi29gOFgF9ErgOuAjqALwA7gdfcSJtvq+r96d8mwzCiiXgZiTwlJ/FppteUnJ9xp+/COiCagSyIEjQGSSpKAKr6BPBETNmXoq57gRsStL0buDsVm275XpzovNjyofqoSjD0+H8xhmGcEpEHeqI1JWeTaqan705MJwSOWAIEQ+EBD8kSsmYXltHBMIy0SZ7RIfPTd/Em93N8kek7tem7LMVEyTCMtEkW6ODLlizhxEzfRQIdwuEB0bTcd9mFiZJhGGmTLNAhJwsyOgysKXHy5tlgSAm4HlLAPKWswkTJMIy0iQQHRAIHYonOMZcpImtK0c6cL2qfUsRDMk8puzBRMgwjbYLhMF6PnBBEEI1znpKSZNfGqBKOF33nGVxTGpi+M08pqzBRMgwjbYIhTbieBNFrN5kTpYggRk8xRjy7/mA4avOsiVI2YaJkGEbaBEI6kEcuHtFrN5li0FM6eU2pJzCYxcE8pezCRMkwjLQJhsMD6zPxGNykmrkH/qCnNFgWGZeJUvZiomQYRtoEQjqwETUe0TnmMkXEU/LE2Tzb0x8cKLPpu+zCRMkwjLQJhsIJI+8gevoucw/8cBxPKTfHGVdX36CnlC37qQwHEyVjSNq7AwMLwoYRIRBKdfouk57SyWmGct2MDp19UZ6STd9lFSZKRkK6+4OcfddTfOXx7ZkeipFlBMJJAh0iodcZfOBrnOm7vBwvMChKuT6PiVKWYaJkJCQyxfH4W7HHZxmTnWAoPLDnJx6RHHOZ3EAbd/rOHdfxXkeUCnJ9JkpZhomSYRhpEwzp0NN3nkjmhOwKdIh4Sl19EVHyWqBDlmGiZBhG2gTCOhDMEI+s2KcUjqwpDZb5PIJHBqfvCvw++kPhjGaeME7ERMkwjLRxpu+SBzpk8qjxeGtKIkJejndAlKb4vahmNvOEcSImSoZhpE2y6btcnzNNlskD9EKRNaWYp1yuz0Nn76CnBBYWnk2YKBmGkTaBcHjI6LvIfqBMBhGE4+S+A0cwI55Svt8RTwt2yB5MlAzDSJtkCVn9rmBl0lOKl5AVIC/HM7imlOt4SiZK2YOJkmEYaeNsnk38+MjLybwoxYu+g/ieUibHaZxISqIkIqtFpE5E6kXkjjjv54rIw+77m0WkKuq9O93yOhG5OplNEal2bex2bfqH6kNESkVkk4h0isi3Y8Z1voi87bb5liQ6/MUwjLQIhnXINEORNaXsmL47sTwvZ3DDbMRTsjWl7CGpKImIF/gOcA1QC9woIrUx1W4B2lR1IXAvcI/bthZYCywBVgPfFRFvEpv3APeqag3Q5tpO2AfQC3wR+Ls4w/8ecCtQ4/5bnezzGoaRnGAoPGRC1sgm1b5g5qLvIvt2Y7+LRgQTotaUTJSyhlQ8peVAvaruVdV+YB2wJqbOGuBB9/pR4ArXK1kDrFPVPlXdB9S79uLadNuscm3g2rxuqD5UtUtVn8cRpwFEZBZQqKovqTO5/JMoW4ZhnAKBJNF3/ogoBbLPU4oEYYAFOmQjqYjSHOBg1OsGtyxuHVUNAu1A6RBtE5WXAsdcG7F9JepjqHE3JBk3ACJyq4hsFZGtTU1NQ5g0DAOc6S7/UNF3vsx7IPH2KcGJntK0vBzApu+yiVREKd7XodidZonqjFR5quNIZUwnF6rep6rLVHVZeXn5ECYNwwBnTWnceEoxT7kp/kFRKnRFaSwCHX63q4lP/mgLPf2Zm9IcD6QiSg3A3KjXlUBshs6BOiLiA4qA1iHaJipvBopdG7F9JepjqHFXJhm3YRjDIJBkTcnrEXweyeyaUpyjKwCm5flOuh6L6bt/3VDHpromtuwf6rFlpCJKrwA1blScHydwYX1MnfXAze719cBGdx1nPbDWjZyrxgk22JLIpttmk2sD1+ZjSfqIi6oeBo6LyIXuWtXHo2wZhnEKBENDR99B5o+FSBQSHvGOYGxF6XB7DwB7mzpHva/xjC9ZBVUNishtwAbAC/xQVbeJyF3AVlVdDzwA/FRE6nG8l7Vu220i8giwHQgCn1HVEEA8m26XtwPrROSrwOuubRL14draDxQCfhG5DrhKVbcDnwZ+DEwBfuP+MwzjFAmGh96nBJCb482SzbMnlkd7SvkDaYZGN/edqg4cl3HoWM+o9jXeSSpKAKr6BPBETNmXoq57gRsStL0buDsVm275XpzovNjyofqoSlC+FTgz3ntGcnTIJTtjsqKqBEI6ZEJWcDylzE7fOT9P8pSmDHpKkbWv0U4c2xMIDQj0oWO9SWpPbiyjg5EY94/asvob0US8isgDPRH+jE/fnXx0BUBhlKc0VgEZHT2Dx6+3dPWNal/jHRMlIymWB8OIJuL9RIdWx8PxlDIvSrGe0gnTd+6hfz2B0fWU2nsCA9dtXYEhahomSkZCzEEy4hHxfpJ5Srm+TK8pOT9jRamkIHfgOj/XEaXuUQ7TjohS5fQptHb3j2pf4x0TJSMhNm1nxCOyITY3y6fvQuH4gQ4LygsAOK00H7/Xg88jA8ejjxYRUaouK6Ctq99Ouh2ClAIdjMmJBToY8YisvyT3lDId6JBon1IOd61ZwjlzixER8v3eUfeUOlxRqiot4Pe7mzneFzwhNN0YxETJMIy0iHhKqYhS5yh7IEMRcUa8caIEP35R1cB1vt9Hd//YeEpVZY6X1trZb6KUAJu+MxJiMwxGPCKeUvJAB292pBlKEqiTn+ula4zWlKpK8wFo6bJ1pUSYKBkJGdQkC78zBons6UnmKeXleOjNwn1KsRT4faOej66jN8C0XB9lU50gi2MW7JAQEyUjIYOLseYyGYNEIuqGyhIOMMXvG/W1mqFItE8plny/d0wCHQqn5FBS4AegrdvCwhNhomQYRlpERCn6XKJ45Pu9dGd0TSn+PqVYCnJHXzw7XFGaHhElm75LiImSkZDBNSWbvjMG6U/RU8r3e+kOhDIW/pzq9N0Uv5euUQ506OgJUjTFR4HfS45XbK/SEJgoGYaRFhFRSrZPKd/vQ3VsziqKR6qBDgV+76ivKbX3BCjMy0FEmJ7vN09pCEyUDMNIi4HpuyTRd5Gjxkd7vSYREU8pdp9SLPl+36iHrh/vDQwkgi0p8NNqopQQEyUjIYOzLhboYAySapqhyAmvmQp2SHR0RSxFU3Lo7AsOZIAYDTp6BzfLTs/302bTdwkxUTISokOeRG9MVvqDqYWERzyl0U52moiIyCTzlIrzc1AdzLowGuPo7AtSOMXJVWCe0tCYKBmGkRZ9Ka8pZXb6LiJKviSuUnG+48EcGyVR6nQP95sW8ZQKciwkfAhMlIyEWEYHIx6pTt9FTnUd7SCCREQCHeKlGYqmeEpk79DoeC8dvY4ARc5xmp7v51h3P+FRnC4cz5goGQmxPxkjHn3BMB5J7oHkZ3hNKZimp9Q+St5LRJSmRa0phXWw3DgREyXDMNKiPxTG7/OkENXmilKG15Q8SUXJ8ZSO9YySp+SeOhu9pgTYulICTJSMhNiZL0Y8evpDTMkZOhwcnDRDQMayOqS8puSGao/WibCD03eRNaXRnS4c76QkSiKyWkTqRKReRO6I836uiDzsvr9ZRKqi3rvTLa8TkauT2RSRatfGbtem/xT6+JyIbBORd0TkIRHJS+/2TG4GYu8s+M6Iors/NLBeNBQFGZ6+i4hSsjWlyP6h0Qp0OO4GOkREqSQ/4inZ9F08koqSiHiB7wDXALXAjSJSG1PtFqBNVRcC9wL3uG1rgbXAEmA18F0R8SaxeQ9wr6rWAG2u7eH0MQf4a2CZqp4JeN16RpqYw2RE0xsIkZck7x0M7lPKVEh4MEVR8nqE6fk5NHf2jco4IqHmkem7yBqWZXWITyqe0nKgXlX3qmo/sA5YE1NnDfCge/0ocIU4E85rgHWq2qeq+4B6115cm26bVa4NXJvXDbMPcA4xnCIiPiAfOJTC5zVcTIyMeHT3B1PylPxeDzleGfAUxppUPSWAGYV5NHb0jso4Ip9/am7MmpJN38UlFVGaAxyMet3glsWto6pBoB0oHaJtovJS4JhrI7avtPpQ1feArwMHgMNAu6o+Fe8DisitIrJVRLY2NTUlvBGTj9RS/xuTi55AaMALGgoRoTAvJ2NRZoNrSskfczOL8jgySqLU3hOgwO/F5yawzfd78fs85iklIBVRivdIiv0OnajOSJWn3YeITMfxoqqB2UCBiNwUpy6qep+qLlPVZeXl5fGqTGrMYzKiSTXQAZz1mtHKlJCMyPRdCo4SMwvzONI+OtN3bd39A8EN4Ih1Sb5ldUhEKqLUAMyNel3JydNgA3XcqbIioHWItonKm4Fi10ZsX+n2cSWwT1WbVDUA/Aq4OIXPa7iYGBnx6AmEBsK9k1E4JYeODE3fhcOK1yNJQ9fBmb5r6eojEBr5jOYtXf2URokSOFN4diR6fFIRpVeAGjcqzo8TLLA+ps564Gb3+npgozrxxOuBtW7kXDVQA2xJZNNts8m1gWvzsWH2cQC4UETy3bWnK4Adqd0WAyz6zohPdzqeUp4vo55SKutJ4EzfqULj8ZH3llq7+gbWkSLMKsrjSPvoTBeOd5KKkrt+cxuwAeeh/oiqbhORu0Tkw261B4BSEakHPg/c4bbdBjwCbAeeBD6jqqFENl1btwOfd22VuraH08dmnICI14C33c963zDukWEYUfSmuKYEEU8pU2tKYbwpfqOaVeTsFjl0rGfEx9HWFaCkIPeEstFcwxrvJA+hAVT1CeCJmLIvRV33AjckaHs3cHcqNt3yvQxGz0WXD6ePLwNfjtfGSI5N3xnxSM9TytyaUiicfONshOqyAgD2NXfxvqqSER1HS1cfJQU5J5TNKsqjtavfDa9P7V5OFiyjg5EQtex3RgyqmtaaUtGUHDp6ghnJDhIKh/F6UxOlOcVTyPEK+5q7htVXouSq3f1BegPhOJ7SFACOmrd0EiZKRkLMUzJi6QuGUR1MIZSMwik++kPhjByJHgxrytN3Pq+HeSX57G3qTLuf3UePc85dT/HPvzl5ybr5uBPMUDr15DUlgMO2rnQSJkqGYaRMJGXQlBQyOsBgap1MTOGFNfVAB4D55VOH5Sn99OV36egN8oPn9tISkxXiULuzRjXb9YwizHRFyYIdTsZEyUiIeUpGLJGUQekEOkBmjmkIhjTlNSWAhRWOKPUF00uLtO1QB9Pcs5Ke3HbkhPcigROzi09Muzmz0A2saB/5wIrxjomSkZDImpJFhBsROgdS5uQkqekwkIE7AyethsKa9NiKaM6aU0QgpNQdOZ5yG1WlvrGTD509m6rSfJ58J5EonegpFeT6KCnwc7C1O+W+JgsmSkZSzGEyIhwfOLAutTWlyFpK7LTWWBDS9Dyls+YUAfBmQ3vKbZo7+2nvCbCwfCpXL5nJS3taaI+aqnzvWA8lBf64EXbVZQXDDqyYyJgoGQmx6Tsjlkhy0VRFqWyqE3XW3Dn22QvS2TwLUDl9CiUFft5uOJZym/pGJzBiYcVUrloyk2BYebauceD9PY1dLCgviNu2qtREKR4mSkZSbPrOiBB7tHcyIpkMRutYiKEIhdITJRHhnLnFbNnXmnKb+qZBUTp3bjHl03LZ4K4rqSp1R49z+oxpcdvOLy/gaEcf3f2ZScOUrZgoGUkxh8mI0NmXnqeU4/VQnJ9DSwY8pZAq3hQyhEdzycIy9rd0p7zWs6exkwK/l1lFeXg8wgdqZ/BsXRO9gRANbT209wRYPDO+KFWVOh7U/mZbV4rGRMlIiE3fGbGkO30HzhReJjylYCic1poSwKWnlwHw+93NKdWvb+xkQcXUgaSvVy+ZSXd/iN/vbubFPY6NC+aXxm0bySKxZxh7oyYyJkpGQiz6zojleG8Ar0dSTjMEUFrgz4inFAgpfl96j7gF5VOZXZTHMzuOplS/vrGTheVTB15fNL+UGYW5fOuZ3fxs8wHmFE+hpmJq3LYLKgrweYQdhzvSGuNEx0TJMIyUOd4bZFqeL6XjICKUTcuMp9QfDJOTYpqhCCLCtUtn8dyupqSH8B3vDXCko5cFUaLj93n4wrW1vP1eO281tPMXKxckvFe5Pi81M6ax7ZCJUjSp++DGpMOm74xYjvcGB471TpWyAj9NmRClUJhpOek/4q47dw7/+ft9/Prtw9x04WkJ6+1pciLnFsZ4Qh8+ezZlBX6O9wW5qnbGkH2dObuQjTsbUdW0hH4iY56SkRDTJCMWx1NKLfIuwsyiKRzvDQ4ESYwVgVCY3DSn7wBqZxWyeOY0/uvld4dMJLvL3WQbL7ru4oVlXL1kZlKhWTK7kJaufo52DC3aXX1BvvFU3aRIS2SiZBhGynT0BtIKcgCYM93JZvBe29im1AmEwuR403/EiQi3XFLNziPHhwx42HX0OLk+J5HrcDmrshiANw62DVnvl6818K2N9Xz119uH3dd4wUTJSEgmjhswspu2rn5K8v3JK0Yxx02x896xsQ19dtaUhveI+/A5s5lRmMs3n9md8O9gV2MnC8qnprUXKpaz5hQxJcfLy3uH3hv1+gFnQ+8776WebWK8YqJkJMQkyYilrbufkqnpiVJlxjwlHbYo5fq8fO7K03n13TYef+vwSe+Hw8rbDcdYMrvwlMbo93lYVjWdl/e2DFlvr5v54UBrN72B9BLGjjdMlAzDSIlwWGnrDqTtKZVPzcXv9dAwCkeND0V/KJx2SHg0NyybS+2sQr7y+HZaYyLx9jZ30tYd4H3Vp35K7YXzS9l55PhJfURQVfY1dVI0JYewwv6WiZ2ayETJSIjN3hnRdPQGCIWV6QXpiZLHI8wqzhBeevcAACAASURBVBtzT6k/GMafZkh4NF6P8C/XL+VYd4C//8WbhKJOl920swmAC6vjb4xNh4sXODZ+t6sp7vutXf109Ab5gBvJt6fRRMmYtJgqGYNEvsmXFKQXfQcwrySfd1vGdk1puIEO0Zw5p4gvXHsGz+xs5B//5x1CYSUUVn7x6kHOmlPEvNLhBzlEOLuymIqonHmxRJK2Xr6oAjBPCQARWS0idSJSLyJ3xHk/V0Qedt/fLCJVUe/d6ZbXicjVyWyKSLVrY7dr038KfRSLyKMislNEdojIRendHsMwIrR1O6I0Pc3pO3D28uxp6iQcHrsvOoFQmJxTmL6LcPPFVfzlygU8tOUAH/nBS/zNutfZdbSTP7t0/giMkoGcec/taoq7XrTfFfPa2YXMKMxlb9MkFyUR8QLfAa4BaoEbRaQ2ptotQJuqLgTuBe5x29YCa4ElwGrguyLiTWLzHuBeVa0B2lzbaffhtvkm8KSqLgbOBnakemMMm74zTiSSKqi0IDfttjUV0+juD/HeGK0rqaqTZugUPaUIf3/1Iv71+qU0tPXw5DtH+PNL5/OhpbNGxDYM5sx7tu7kKbx3W7rweoQ5xVOoKi0wTwlYDtSr6l5V7QfWAWti6qwBHnSvHwWuEGfX2Bpgnar2qeo+oN61F9em22aVawPX5nXD6UNECoFLgQcAVLVfVVM/KMWwyTvjBCKeUnF++tN3NTOcrAeR84dGm/5QGOCUAh2iERFuWDaXl+5cxc6vrObOPzhjRDMwXLyglIppufxi68GT3nu3pZvZxXn4fR6qywrYP8HPYErlf2wOEH2nGtyyuHVUNQi0A6VDtE1UXgocc23E9pVuH/OBJuBHIvK6iNwvInFP2xKRW0Vkq4hsbWqKv9hoGJOdSNaB8mnpe0qRpKW7G1M/avxUCIScr1Tp5r5LhojgGyHvKxqf18P151eyqa7xpKwN77Z0cVqJ8+iqLiugpav/hNNtJxqp3N14/6uxX6IT1Rmp8uH04QPOA76nqucCXcBJ62EAqnqfqi5T1WXl5eXxqkxKbPrOiOZwe2/Co72TMb3Az4zC3DFLPtrT76zNpJPNPNN85H1zCSs8tOXAQJmqsre5i9PcgIqqssgZTBPXW0pFlBqAuVGvK4FDieqIiA8oAlqHaJuovBkodm3E9jWcPhpUdbNb/iiOSBkpEtnJbnkiDYCjHb3MLMwbdvtz5hbzxsGxmUEfECX/+Mk5fVppAVeeMYMHX9o/kCewoa2H471BlswuAgbPYJrI60qpiNIrQI0bFefHCSpYH1NnPXCze309sFGdJ9p6YK0bOVcN1ABbEtl022xybeDafGw4fajqEeCgiCxy21wBTPzEUaOAeUwGwJH2XmYWDV+Uzps3nXdbumkZg4zhPYHx5ykB3LZqIce6A/zkpf0AvNXgpBWKZI6YV5KPyGCY+EQkqSi56ze3ARtwotceUdVtInKXiHzYrfYAUCoi9cDncafJVHUb8AiOGDwJfEZVQ4lsurZuBz7v2ip1bafdh9vmr4CfichbwDnA/5fuDZrMmBYZ0RzpODVROnfedGAwj9to0t3veBr5/vElSufMLeaKxRV8e2M9B1u7eW5XI9PyfNS6opSX42V20ZQJLUop+baq+gTwREzZl6Kue4EbErS9G7g7FZtu+V6c6LzY8uH08QawLF4bIzkRD8mm74zeQIjWrv5Tmr5bWlmE3+fhpb0tXJnknKFTJeIpDWf9K9P805olXH3v7/jYA5s53N7Lh86efcIm4IkegWcZHYyk2PSdcdiNCJt1Cp5SXo6Xi+aXsmln40gNKyGRNaXx5ikBVE7P5wcfW0ZHb5AZhXl89sqaE96vLitgX3PXhM3iP35WAY0xR20Cz3CJfDOPLLQPl8sXlfP//nc7+5u7BiLJRoOIpzQeRQngkpoyXv3HKwFO2g9VVVZAR2+Q1q5+SqemH56f7ZinZCTGpu8Ml8gaxqkKyRVnONN2v3775OMgRpLu/vE7fRdBROJu0K0uc8LDJ2oEnomSYRhJ2d/SxbRcH6VpZgiPZW5JPsurSnj01YZRnX4az9N3yaguczYi72se2wS3Y4WJkpEQm7wzIuxzp9tGIrXO9csq2dfclfS01VPheK+T8WBaXvopkbKdyulTyPHKmGXHGGtMlIyETNB1VGMY7G3qOuX1pAgfWjqbsql+vvts/YjYi8ex7gBTcrwjlvsum8jxejhjViFvHZyYR6NPvP8xwzBGlPbuAO8d6+GMWad29HeEKX4vf/b++fx+dzObkxwDPlzaewLDShw7XlhaWcQ777WP6VEgY4WJkpEQi74zALYddr6RRzZwjgQfu+g05hRP4f/+99v0BU8+Q+hUOdYToGjKRBalYo73Bdk7AfcrmSgZCbHpOwNgu5tEtXaEPCWAfL+Pu//Pmexp6uKrj4/8MWftPQEKJ7AonV1ZDMCbY5RLcCwxUTIMY0jeea+dimm5wzqyYihWLqrg1kvn89OX3+X+3+8dUdvHuvspnsCitLBiKgV+L68fbMv0UEYc2zxrJMQcJUNV2bKvlWVV00fF/u2rF3OwtZuv/noHHb1BPntFDR7PqUf4He3o48L5pSMwwuzE6xGWVZWMagRjpjBPyUjIRE1jYqROQ1sPh9p7uaB6dB7wXo/wHzeeyx+dV8m3ntnNzT/aQkPbqe2/6e4P0t4TOKXkseOBC+eXUt/YSePx3uSVxxEmSoZhJORlNzrugvklo9aHz+vh6zcs5Z//8Cxe2d/KFf/2HF/fUEdbV/+w7B0ZgTx944GLFjhfFDZPMG/JRMlIiPlJxjM7GqmYlsvpFdNGtR8R4cbl89j4tyu5eslMvr2pnou+9gxffuwd9jZ1pmXr3VbH06qcnj8aQ80azpxdyNRcHy+NUlh9pjBRMhJjqjSp6Q2EeG5XE1ctmTEi6zypMLt4Ct+68Vye/tylfGjpbH6+5QCr/u05/uT+l/nN24cJhMJJbdQdcTIdjLaQZhqf18Py6hJerG/O9FBGFBMlIyGRfUqCZWSdjDxb10hPIMTVS2aOed81M6bxrzeczQt3rOLvrjqd/c3dfPpnr7Hiaxv5xlN1HG7vSdj2nffamV2UR9EE3jwb4dKaMva3dE+o85VMlIyEhN0vpWP0JdnIMta9cpAZhblclMEotoppedy2qobf/cPl3P/xZSyZXch/bKpnxdc28mc/2cpzu5pOyGoQDIV5ob6ZCxdM3Mi7aFYuqgCcLxATBQsJNxISdqPvRiIJpzG+cI7ibuKvLl+Iz5v5765ej3Bl7QyurJ3BwdZufr7lAI+8cpCntx9lXkk+H71gHn+8bC4bth2hrTvAtWfNyvSQx4SqsgLmlxWwqa6JT6yozvRwRgQTJSMhEzCtlpEiP/jdHnweYe3yeZkeyknMLcnn9tWL+eyVNWzYdpT/evldvvabnfzLkzsJKyyvKuFy14OYDFy+uIKfvvwuPf0hpkyAozpMlIyERPYpeTL/RdkYQxraunn4lYP88bK5zC6ekunhJCTX5+XDZ8/mw2fPZtfR4/zqtfeYlufj5ourxiwwIxu4fFEFDzy/j5f2NrNq8YxMD+eUSelxIyKrRaROROpF5I447+eKyMPu+5tFpCrqvTvd8joRuTqZTRGpdm3sdm36h9uH+55XRF4XkcdTvy0GDHpKFugwubjrf7fjEeEzly/M9FBS5vQZ07jjmsV85vKFTM2dXN+131c9nXy/l407J8a6UlJREhEv8B3gGqAWuFFEamOq3QK0qepC4F7gHrdtLbAWWAKsBr7risRQNu8B7lXVGqDNtZ12H1Fj+xtg5DM+TgIG15QyPBBjzHji7cM8tf0on/vA6VntJRmD5Pq8rFhYxqadTRMiC0sqntJyoF5V96pqP7AOWBNTZw3woHv9KHCFOKvja4B1qtqnqvuAetdeXJtum1WuDVyb1w2zD0SkErgWuD+122FEExElj6nSpGBPUye3P/oWSyuL+NNLJsai+WRh1eIK3jvWw66j6W00zkZSEaU5wMGo1w1uWdw6qhoE2oHSIdomKi8Fjrk2YvtKtw+Afwf+ARhyx52I3CoiW0Vka1NT01BVJxU6MH1nTHSajvfxZw9uJcfn4Xs3nZ8VEXdG6kQCOybCFF4qv3nxnkmxPmKiOiNVnnYfIvJBoFFVX43z/omVVe9T1WWquqy8vDxZ9UnDwCF/pkoTmubOPv7k/pc53N7LfR87nzk2bTfumFmUx5LZhWzceTTTQzllUhGlBmBu1OtK4FCiOiLiA4qA1iHaJipvBopdG7F9pdvHCuDDIrIfZ3pwlYj8Vwqf13CJbJ41TZq47DzSwZpvv8CB1m5++In3saxq9BKvGqPLFYsrePXdNo51Dy+RbbaQiii9AtS4UXF+nKCC9TF11gM3u9fXAxvVWXFbD6x1I+eqgRpgSyKbbptNrg1cm48Npw9VvVNVK1W1yrW/UVVvSvG+GNjm2YmMqvKLrQf5o+++SDAc5pE/v2gg67QxPrl8cQVhhed2je8liKSxk6oaFJHbgA2AF/ihqm4TkbuAraq6HngA+KmI1ON4L2vdtttE5BFgOxAEPqOqIYB4Nt0ubwfWichXgddd2wynD+PUiKwpTaItH5OCxo5evvA/7/D09qMsry7hW2vPnfBnD00Gzq4sprTAz8adjaw5J3bZf/yQUkC/qj4BPBFT9qWo617ghgRt7wbuTsWmW74XN3oupjztPqLefxZ4NtH7RnwGPCWbwJsQ9AZCPPD8Pr6zqZ5gWPnHa8/gUyuqJ9VG04mMxyOsXFTBb3ccJRgKj9tglcm1y8xIi4HNs/bMGtcEQmH++7X3+NbG3TS09XBV7Qy+cO0ZnFZakOmhGSPMFWdU8MvXGnj94DHeN07XB02UjITYmtL4pi8Y4tFXG/jes3toaOvhzDmF3PNHS1mxsCzTQzNGiUtqyvB5hGd2NJooGRMPHZi+M8YTzZ19PLT5AP+1+V2OdvRxztxi7lqzhMsXVdgXjAlOYV4Oy6tL2LSzkTuuWZzp4QwLEyUjITZ9N7548+AxHnxxP4+/dZj+UJj315Tx9RvO5pKFZSZGk4hViyv46q930NDWPS6PhDdRMhKilvsu6+nsC/Lrtw7x0JaDvHHwGAV+Lzcun8vHL65iQfnUTA/PyACXu6K0aWcjH7uoKtPDSRsTJSMhliU8O1FVtuxr5ZGtDTzx9mF6AiEWlBfw5Q/Vcv35lUzLm/jHgBuJmV9WQFVpPs+YKBkTjcGErBkeiAHAkfZefvlaA7/YepD9Ld1MzfVx3bmzuWHZXM6dW2xTdAbgBCZdvriCn20+QHd/kHz/+HrMj6/RGmPKQBZ8e9hljOO9AZ7adpT/eeM9XqhvJqxwQXUJf7WqhmvOmjnuHjjG2HDF4hn86IX9vFjfwpW14+vgP/uNNhIStui7jNAfDPNsXSOPvXmI324/Sl8wTOX0KfzlyoVcf34lVWW2v8gYmuXVJRT4vWysazRRMiYOQXdRyWfzd6NOOKxs2d/KY2+8xxNvH6G9J0BJgZ+PvG8ua86ZzXnzptv0nJEyfp+H99eUs2lnI6o6rn53TJSMhARDjih5TZRGBVVl++EO1r9xiPVvHuJwey/5fi9X1c5gzblzuGRhGTnjNFWMkXlWLa7gyW1H2HH4OLWzCzM9nJQxUTISEnTPrrCTZ0eOiBA98fZhnnj7CPuau/B5hMtOL+eOaxbzgdoZtk5kjAiXLXLOhnu+vslEyZgYRKbv9KQzHY10UFW2HYoI0WH2t3Tj9QgXzS/lT99fzTVnzqKkwJ/pYRoTjBmFeSysmMrz9S3ceumCTA8nZUyUjIQEQ46npKZJaRMRol+/fZjfRAnRxQtK+fPLFnBV7QxKp+ZmepjGBOeShWWse+UAfcEQuT5vpoeTEiZKRkICIVOjdFBV3nnPFaJ3DvNulBD9xWULuGrJTPOIjDFlxcIyfvzifl4/cIwL54+PQxxNlIyEhAam74xEhMLK6wfa2LDtCBu2HeVAazc+j3DxwjL+cuUCrqqdyXQTIiNDXDC/BK9HeKG+2UTJGP9EAh2ME+kLhnixvoWnth/h6e1Hae7sx+/1cNGCUm67fCEfqJ1hQmRkBYV5OZxdWcTz9c387VWLMj2clDBRMhIyMH1nrhLHewNsqmviqW1HeLauic6+IFNzfaxcVM7VS2ayclG55ZwzspIVC8v4zqZ6OnoDFI6D31ETpQnO6wfaqJyeT/m09BfVQ5M8+q7xeC+/3d7Ihm1HeHFPM4GQUjbVz4fOnsVVS2Zy8YLScbN4bExeViws4z821vPynhauWjIz08NJionSBOf/fPdFKqblsuULV6bdNhCafNN3+5q7eHq7sz702oE2VGFeST6fuLiKq5fM5Nx5020zsTGuOHdeMVNyvLxQ3zxxRElEVgPfBLzA/ar6tZj3c4GfAOcDLcBHVHW/+96dwC1ACPhrVd0wlE0RqQbWASXAa8DHVLU/3T5EZK5bfyYQBu5T1W+me4PGM2HX02k83jes9n3BiR8S3hcMsWVfK5t2NrGprpF9zV0ALJldyOeuPJ2rlsxg0Yxp4ypNi2FEk+vzsry6hBf2tGR6KCmRVJRExAt8B/gA0AC8IiLrVXV7VLVbgDZVXSgia4F7gI+ISC2wFlgCzAZ+KyKnu20S2bwHuFdV14nI913b3xtGH0Hgb1X1NRGZBrwqIk/HjHtC0x0InVL7nv5Ta5+tHO3oZdPORjbubOSF+ma6+kP4fR4uml/KJy6uYtXiCuaWjL8TOw0jERctKOVrv9lJc2cfZVm+Py4VT2k5UK+qewFEZB2wBoh+uK8B/p97/SjwbXG+Wq4B1qlqH7BPROpde8SzKSI7gFXAR906D7p2v5duH6r6EnAYQFWPu7bnxIx7QtPVFzyl9j2uqI13RykUVt5sODYgRNsOdQAwuyiP686dw6rFFVy8oIwpflsfMiYmy6tLANiyr5U/OGtWhkczNKmI0hzgYNTrBuCCRHVUNSgi7UCpW/5yTNs57nU8m6XAMVUNxqk/nD4AEJEq4Fxgc7wPKCK3ArcCzJs3L16VcUnnKYpSt+sp6Ticv2vvCfC7XU1s2tnIs7uaaO3qxyNw/mnT+YfVi1i1uMKm5YxJw1lzipiS42Xz3pYJIUrx/mpjn1KJ6iQqj5f6eKj6w+nDaSQyFfgl8FlV7YhTF1W9D7gPYNmyZePvCZyAU/WUek9x+m8sUVV2He1k485GNtU18uq7bYTCyvT8HFYuqmDlonIuO72c4nzbP2RMPnK8Hs4/bTqb97VmeihJSUWUGoC5Ua8rgUMJ6jSIiA8oAlqTtI1X3gwUi4jP9Zai66fdh4jk4AjSz1T1Vyl81gnF8V5HlPy+4R1/0N3vtM9Wle7pD/HS3mZHiHY28d6xHgBqZxXy6csWcPniCs6ZW2zRcoaBc2LxN367i2Pd/Vn95SwVUXoFqHGj4t7DCSr4aEyd9cDNwEvA9cBGVVURWQ/8XES+gROEUANswfFuTrLpttnk2ljn2nxsOH24600PADtU9Rvp3piJQEdPAID8Ya6VtLvts2n2rqGte2Bt6MU9LfQFw+T7vaxYWMZtqxZy+aIKZhblZXqYhpF1LK8uQRVe2d/GB7L4NNqkouSu39wGbMAJ3/6hqm4TkbuAraq6Hufh/1M3yKAVR2Rw6z2CE1wQBD6jqiGAeDbdLm8H1onIV4HXXduk24eIXAJ8DHhbRN5wbfxfVX1ieLdq/BHxlPJz0hel7v4gvYHM71MKhsK8+m4bG+sa2bSzkV1HOwE4rTSfG5fPY9XiCi6YX2KbWA0jCWfPLcbv87B5b8v4FiUA90H+REzZl6Kue4EbErS9G7g7FZtu+V4GI/Siy9PqQ1WfJ/5606Sho9fxdPKG4Sm1dPYPXI+1o9TS2cezdU1srGvkd7uaON4bJMcrLK8u4Y+XzeXyxRXMLyuwIAXDSIO8HC/nzC1my/7sXleyjA4TmNYuR1iG40UMd8PtcIicPbTRnZZ7s+EYqlA+LZdrzpzJqsUVrFhYZrnlDOMUubC6hG9vqud4byBr/55MlCYwh9yF/0hmh3TY72Y2mFeSPyqLSp19QZ7f3cwmN1qu8XgfIrC0spjPXnE6qxZXsGR2IR4LUjCMEWN5dSnhjfW8+m4bKxdVZHo4cTFRmsAcOtYLQGAYR1DUN3Xi9QjzSvIHpgFPlYOt3fx2x1E27mzk5b0tBELKtFwfl55ezuWLnbDtbN9tbhjjmfNOK8bnETbvazVRMsaWcFjZcdjZlhUahqe0dX8rZ84pIsc7fE8lGArz+sFjPLOjkWd2HGV3oxOksKC8gE+uqObyRRUsq5pOjnd4IeuGYaRHvt/H0soiNu/N3jx4JkoTlG2HOjjeF8QjEEzzWPPGjl5eO3CMv7hsPtsPdaQ1e9cbCPFsXSMbth1lU10jx7oD+DzCBfNLWLt8HlcsrqCqrCDNT2MYxkixvLqU+3+/l57+UFam1jJRmqD86MV9+L0erqytYOv+trTa3ve7vYTCyvXnz+Urh7cnPU+pp98Rol+/fZiNOxvp7g8xPT+HVYsruGLxDN5/etm4OFzMMCYDF8wv4fvP7eG1A22sWFiW6eGchInSBOShLQf41Wvv8ReXLaCrL0gwjem7p7cf5Ycv7OPG5XOpLivA65GEntbbDe2se+UA6984xPG+IKUFfq47dw7XnjWLC6pL8Nm0nGFkHctOm45HYPPeFhMlY3TpC4b45yd28uMX93PZ6eV87gM1fO03O1M6rC8UVu7//V7+dUMdZ1UW84VrawHI9XnoD57Y/oX6Zr75zG627Gsl1+fh2qWzuP68Si6YX2opfQwjy5mWl8OS2UW8nKV58EyUJgjP1jXyT/+7nX3NXXxqRTV3XLMYv8+DbwhPB5w9Qs/tauKfn9hJ3dHjXHPmTO65filTc51fjVyfd+Cwv9aufr74P+/w67cPM6Mwly9+sJbrz6+kaIpNzRnGeOKC6hJ+8vK79AZC5A0j48toYqI0ztm6v5X/2FjPc7uamF9WwE8+tZxLTy8feN/n9RCMExLeGwjxv28e4kcv7Gf74Q7mleTznY+exx+cNfOETAl+n4e+YJgj7b189D9fpqGth7+76nT+9P3zs+6X2TCM1Lhgfin3P7+PNw8e44L5pZkezgmYKI1DAqEwz+xo5Ecv7GPzvlZKCvzcec1iPrmi+qSM4DkeIRBSVBURofF4Lz97+QA/2/wuzZ39nD5jKl/7w7P4w/Mq42YTz/V56OkP8umfvcrRjl5+/mcXsKyqZKw+qmEYo8DyqhJEYPO+VhMlY/jUN3byi1cP8stXG2ju7GdWUR5f/GAtNy6fS74//n9lJNigtauff3t6F7/YepBASFm1uIJPrahmxcLSIXPI5eZ46OoP8fqBY/z7R84xQTKMCUBRfg6LZkxj874WnIMVsgcTpSxGVXmroZ0N247w1Paj1Dc6WRZWLa7gxuVzubSmPGmEm8/d/PpH33uRg2093Lh8Lp9aUc388qkpjSGSN6+mYiprzpl9ah/IMIys4cL5pax75QD9wfCwz1wbDUyUsoxAKMzmva08tf0IT207ypGOXrwe4YLqEj524Wlcc+ZMKgpTPy8ox+P8su1v6T5pvSkVzptXDMBfX1FjWbkNYwJxQXUJP35xP2+/1875p03P9HAGMFHKArr7g/xuVxMbth3lmR1H6egNkpfj4bLTy/n72kVccUbFsE+KjOzYPruyKG1BAli5qIItX7iCiml2cJ5hTCSWVztT8Zv3tZgoGdDW1c9vdxxlw7aj/H53E33BMMX5OXygdiZXL5nB+2vKRyQFyIqFZVROn8KXP7xk2DZMkAxj4lE6NZeFFVPZvLeVv1yZ6dEMYqI0hjR29LJh2xGe3HaEl/e2Egors4vyuHH5PK5aMoPlVSOfBaG6rIDnb181ojYNw5gYXFBdwmNvHCIYCmdNBhYTpVGmNxBiw7Yj/GJrAy/saUYV5pcX8BeXzefqJTM5a06RrdUYhpERLphfys82H2D74Q6WVhZnejiAidKocay7nx+/uJ8fv7ifY90BKqdP4a9X1fDBpbOomTEt08MzDMPgQndd6aU9LSZKE5VQWPn55nf5lw11HO8NcuUZM/jkiiouml9qp6gahpFVVBTmUTurkKe3H+XPL1uQ6eEAkNIkooisFpE6EakXkTvivJ8rIg+7728Wkaqo9+50y+tE5OpkNkWk2rWx27XpH+k+Rotj3f184kdb+OJj2zi7spgnP/t+7r95GSsWlpkgGYaRlaw+cyavHmij8XhvpocCpCBKIuIFvgNcA9QCN4pIbUy1W4A2VV0I3Avc47atBdYCS4DVwHdFxJvE5j3AvapaA7S5tke6jxGnraufG77/Epv3tvLPf3gWP71lOYtnFo5Wd4ZhGCPC1UtmogqPv3k400MBUvOUlgP1qrpXVfuBdcCamDprgAfd60eBK8RZvV8DrFPVPlXdB9S79uLadNuscm3g2rxuJPtI7bakR38wzCd+/Arvtnbz40++jxuXz7PgBcMwxgWnz5jK8qoSvrVxN4eO9WR6OCmtKc0BDka9bgAuSFRHVYMi0g6UuuUvx7Sd417Hs1kKHFPVYJz6I9XHSYjIrcCt7stOEWkBmuPVTcaKu4fTKmspY5j3YQJi92IQuxcOE+4+zPnysJuWAaeNxBhSEaV4X/ljD+hJVCdReTwPbaj6I9nHyYWq9wH3RV6LyFZVXRav7mTC7sMgdi8GsXvhYPdhEPdeVI2ErVSm7xqAuVGvK4FDieqIiA8oAlqHaJuovBkodm3E9jVSfRiGYRhZSiqi9ApQ40bF+XGCCtbH1FkP3OxeXw9sVFV1y9e6kXPVODnStySy6bbZ5NrAtfnYSPaR2m0xDMMwMkHS6Tt3/eY2YAPgBX6oqttE5C5gq6quBx4Afioi9Tjey1q37TYReQTYDgSBz6hqCCCeTbfL24F1ka+arAAAA89JREFUIvJV4HXXNiPcRzLuS15lUmD3YRC7F4PYvXCw+zDIiN0LcZwNwzAMw8g82ZGBzzAMwzAwUTIMwzCyCBOlKMY6LVEmEJEfikijiLwTVVYiIk+7qZ2eFpHpbrmIyLfc+/GWiJwX1eZmt/5uEbk5Xl/ZjIjMFZFNIrJDRLaJyN+45ZPxXuSJyBYRedO9F//klo9Yyq/xhJsR5nURedx9PVnvw34ReVtE3hCRrW7Z6P99qKr9c9bVvMAeYD7gB94EajM9rlH4nJcC5wHvRJX9C3CHe30HcI97/QfAb3D2gl0IbHbLS4C97s/p7vX0TH+2NO/DLOA893oasAsnHdVkvBcCTHWvc4DN7md8BFjrln8f+LR7/ZfA993rtcDD7nWt+3eTC1S7f0/eTH++YdyPzwM/Bx53X0/W+7AfKIspG/W/D/OUBhmztESZRFV/hxO9GE10CqfY1E4/UYeXcfaQzQKuBp5W1VZVbQOexsk7OG5Q1cOq+pp7fRzYgZMJZDLeC1XVTvdljvtPGbmUX+MGEakErgXud1+PZOqzicCo/32YKA0SL53SnAR1JxozVPUwOA9roMItT3RPJtS9cqddzsXxECblvXCnrN4AGnEeHHtIMeUXEJ3ya7zfi38H/gEIu69TTn3GxLoP4HwxeUpEXhUnFRuMwd+Hnac0SCrplCYb6aZ2GneIyFTgl8BnVbVDEifSndD3Qp29feeISDHw38AZ8aq5PyfkvRCRDwKNqvqqiKyMFMepOqHvQxQrVPWQiFQAT4vIziHqjti9ME9pkMmcluio62rj/mx0yyd0CicRycERpJ+p6q/c4kl5LyKo6jHgWZx1gZFK+TVeWAF8WET240zfr8LxnCbbfQBAVQ+5PxtxvqgsZwz+PkyUBpnMaYmiUzjFpnb6uBtZcyHQ7rrsG4CrRGS6G31zlVs2bnDn/h8AdqjqN6Lemoz3otz1kBCRKcCVOGtsI5Xya1ygqneqaqU6iUXX4nyuP2GS3QcAESkQkWmRa5zf63cYi7+PTEd4ZNM/nAiSXTjz6V/I9HhG6TM+BBwGAjjfYm7BmQd/Btjt/ixx6wrOQYl7gLeBZVF2PoWzgFsPfDLTn2sY9+ESnGmEt4A33H9/MEnvxVKclF5vuQ+eL7nl83EepvXAL4BctzzPfV3vvj8/ytYX3HtUB1yT6c92CvdkJYPRd5PuPrif+U3337bI83As/j4szZBhGIaRNdj0nWEYhpE1mCgZhmEYWYOJkmEYhpE1mCgZhmEYWYOJkmEYhpE1mCgZhmEYWYOJkmEYhpE1/P9AdvL9XGINDwAAAABJRU5ErkJggg==\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, 1.5e-6)\n",
"# plt.yscale('log')\n",
"# plt.xlim(770, 785)\n",
"plt.savefig('test.png')\n",
"# print(jpsi_width)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"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": 14,
"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": 15,
"metadata": {},
"outputs": [],
"source": [
"total_f.update_integration_options(draws_per_dim=200000, mc_sampler=None)\n",
"# inte = total_f.integrate(limits = (x_min, x_max), norm_range=False)\n",
"# inte_fl = zfit.run(inte)\n",
"# print(inte_fl/4500)\n",
"# print(pdg[\"jpsi_BR\"]/pdg[\"NR_BR\"], inte_fl*pdg[\"psi2s_auc\"]/pdg[\"NR_auc\"])"
]
},
{
"cell_type": "code",
"execution_count": 16,
"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",
"## Mixture distribution for sampling"
]
},
{
"cell_type": "code",
"execution_count": 17,
"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": 18,
"metadata": {},
"outputs": [],
"source": [
"# total_f._sample_and_weights = UniformSampleAndWeights"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"# 0.00133/(0.00133+0.213+0.015)*(x_max-3750)/(x_max-x_min)"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [],
"source": [
"# zfit.settings.set_verbosity(10)"
]
},
{
"cell_type": "code",
"execution_count": 21,
"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": 22,
"metadata": {},
"outputs": [],
"source": [
"# with open(r\"data/zfit_toys/toy_0/0.pkl\", \"rb\") as input_file:\n",
"# sam = pkl.load(input_file)\n",
"# print(sam[:10])\n",
"\n",
"# with open(r\"data/zfit_toys/toy_0/1.pkl\", \"rb\") as input_file:\n",
"# sam2 = pkl.load(input_file)\n",
"# print(sam2[:10])\n",
"\n",
"# print(np.sum(sam-sam2))"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"# print(\"Time to generate full toy: {} s\".format(int(time.time()-start)))\n",
"\n",
"# total_samp = []\n",
"\n",
"# for call in range(calls):\n",
"# with open(r\"data/zfit_toys/toy_0/{}.pkl\".format(call), \"rb\") as input_file:\n",
"# sam = pkl.load(input_file)\n",
"# total_samp = np.append(total_samp, sam)\n",
"\n",
"# total_samp = total_samp.astype('float64')\n",
"\n",
"# data2 = zfit.data.Data.from_numpy(array=total_samp[:int(nevents)], obs=obs)\n",
"\n",
"# data3 = zfit.data.Data.from_numpy(array=total_samp, obs=obs)\n",
"\n",
"# print(total_samp[:nevents].shape)"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [],
"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": 25,
"metadata": {},
"outputs": [],
"source": [
"# sampler = total_f.create_sampler(n=nevents)\n",
"# nll = zfit.loss.UnbinnedNLL(model=total_f, data=sampler, fit_range = (x_min, x_max))\n",
"\n",
"# # for param in pdf.get_dependents():\n",
"# # param.set_value(initial_value)\n",
"\n",
"# sampler.resample(n=nevents)\n",
"\n",
"# # Randomise initial values\n",
"# # for param in pdf.get_dependents():\n",
"# # param.set_value(random value here)\n",
"\n",
"# # Minimise the NLL\n",
"# minimizer = zfit.minimize.MinuitMinimizer(verbosity = 10)\n",
"# minimum = minimizer.minimize(nll)"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [],
"source": [
"# jpsi_width"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [],
"source": [
"# plt.hist(sample, weights=1 / prob(sample))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Fitting"
]
},
{
"cell_type": "code",
"execution_count": 28,
"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": 29,
"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": 30,
"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": 31,
"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": 32,
"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": 33,
"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": 34,
"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": 35,
"metadata": {},
"outputs": [],
"source": [
"# 0.15**2*4.2/1000\n",
"# result.hesse()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Constraints"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [],
"source": [
"# 1. Constraint - Real part of sum of Psi contrib and D contribs\n",
"\n",
"sum_list = []\n",
"\n",
"sum_list.append(ztf.to_complex(jpsi_s) * tf.exp(tf.complex(ztf.constant(0.0), jpsi_p)) * ztf.to_complex(jpsi_w / (tf.pow(jpsi_m,3))))\n",
"sum_list.append(ztf.to_complex(psi2s_s) * tf.exp(tf.complex(ztf.constant(0.0), psi2s_p)) * ztf.to_complex(psi2s_w / (tf.pow(psi2s_m,3))))\n",
"sum_list.append(ztf.to_complex(p3770_s) * tf.exp(tf.complex(ztf.constant(0.0), p3770_p)) * ztf.to_complex(p3770_w / (tf.pow(p3770_m,3))))\n",
"sum_list.append(ztf.to_complex(p4040_s) * tf.exp(tf.complex(ztf.constant(0.0), p4040_p)) * ztf.to_complex(p4040_w / (tf.pow(p4040_m,3))))\n",
"sum_list.append(ztf.to_complex(p4160_s) * tf.exp(tf.complex(ztf.constant(0.0), p4160_p)) * ztf.to_complex(p4160_w / (tf.pow(p4160_m,3))))\n",
"sum_list.append(ztf.to_complex(p4415_s) * tf.exp(tf.complex(ztf.constant(0.0), p4415_p)) * ztf.to_complex(p4415_w / (tf.pow(p4415_m,3))))\n",
"sum_list.append(ztf.to_complex(DDstar_s) * tf.exp(tf.complex(ztf.constant(0.0), DDstar_p)) * (ztf.to_complex(1.0 / (10.0*tf.pow(Dstar_m,2)) + 1.0 / (10.0*tf.pow(D_m,2)))))\n",
"sum_list.append(ztf.to_complex(Dbar_s) * tf.exp(tf.complex(ztf.constant(0.0), Dbar_p)) * ztf.to_complex(1.0 / (6.0*tf.pow(Dbar_m,2))))\n",
"\n",
"sum_ru_1 = ztf.to_complex(ztf.constant(0.0))\n",
"\n",
"for part in sum_list:\n",
" sum_ru_1 += part\n",
"\n",
"sum_1 = tf.math.real(sum_ru_1)\n",
"constraint1 = zfit.constraint.GaussianConstraint(params = sum_1, mu = ztf.constant(1.7*10**-8), \n",
" sigma = ztf.constant(2.2*10**-8))\n",
"\n",
"# 2. Constraint - Abs. of sum of Psi contribs and D contribs\n",
"\n",
"sum_2 = tf.abs(sum_ru_1)\n",
"constraint2 = tf.cond(tf.greater_equal(sum_2, 5.0e-8), lambda: 100000., lambda: 0.)\n",
"\n",
"# 3. Constraint - Maximum eta of D contribs\n",
"\n",
"constraint3_0 = tf.cond(tf.greater_equal(tf.abs(Dbar_s), 0.2), lambda: 100000., lambda: 0.)\n",
"\n",
"constraint3_1 = tf.cond(tf.greater_equal(tf.abs(DDstar_s), 0.2), lambda: 100000., lambda: 0.)\n",
"\n",
"# 4. Constraint - Formfactor multivariant gaussian covariance fplus\n",
"\n",
"Cov_matrix = [[ztf.constant( 1.), ztf.constant( 0.45), ztf.constant( 0.19), ztf.constant(0.857), ztf.constant(0.598), ztf.constant(0.531), ztf.constant(0.752), ztf.constant(0.229), ztf.constant(0,117)],\n",
" [ztf.constant( 0.45), ztf.constant( 1.), ztf.constant(0.677), ztf.constant(0.708), ztf.constant(0.958), ztf.constant(0.927), ztf.constant(0.227), ztf.constant(0.443), ztf.constant(0.287)],\n",
" [ztf.constant( 0.19), ztf.constant(0.677), ztf.constant( 1.), ztf.constant(0.595), ztf.constant(0.770), ztf.constant(0.819),ztf.constant(-0.023), ztf.constant( 0.07), ztf.constant(0.196)],\n",
" [ztf.constant(0.857), ztf.constant(0.708), ztf.constant(0.595), ztf.constant( 1.), ztf.constant( 0.83), ztf.constant(0.766), ztf.constant(0.582), ztf.constant(0.237), ztf.constant(0.192)],\n",
" [ztf.constant(0.598), ztf.constant(0.958), ztf.constant(0.770), ztf.constant( 0.83), ztf.constant( 1.), ztf.constant(0.973), ztf.constant(0.324), ztf.constant(0.372), ztf.constant(0.272)],\n",
" [ztf.constant(0.531), ztf.constant(0.927), ztf.constant(0.819), ztf.constant(0.766), ztf.constant(0.973), ztf.constant( 1.), ztf.constant(0.268), ztf.constant(0.332), ztf.constant(0.269)],\n",
" [ztf.constant(0.752), ztf.constant(0.227),ztf.constant(-0.023), ztf.constant(0.582), ztf.constant(0.324), ztf.constant(0.268), ztf.constant( 1.), ztf.constant( 0.59), ztf.constant(0.515)],\n",
" [ztf.constant(0.229), ztf.constant(0.443), ztf.constant( 0.07), ztf.constant(0.237), ztf.constant(0.372), ztf.constant(0.332), ztf.constant( 0.59), ztf.constant( 1.), ztf.constant(0.897)],\n",
" [ztf.constant(0.117), ztf.constant(0.287), ztf.constant(0.196), ztf.constant(0.192), ztf.constant(0.272), ztf.constant(0.269), ztf.constant(0.515), ztf.constant(0.897), ztf.constant( 1.)]]\n",
"\n",
"def triGauss(val1,val2,val3,m = Cov_matrix):\n",
"\n",
" mean1 = ztf.constant(0.466)\n",
" mean2 = ztf.constant(-0.885)\n",
" mean3 = ztf.constant(-0.213)\n",
" sigma1 = ztf.constant(0.014)\n",
" sigma2 = ztf.constant(0.128)\n",
" sigma3 = ztf.constant(0.548)\n",
" x1 = (val1-mean1)/sigma1\n",
" x2 = (val2-mean2)/sigma2\n",
" x3 = (val3-mean3)/sigma3\n",
" rho12 = m[0][1]\n",
" rho13 = m[0][2]\n",
" rho23 = m[1][2]\n",
" w = x1*x1*(rho23*rho23-1) + x2*x2*(rho13*rho13-1)+x3*x3*(rho12*rho12-1)+2*(x1*x2*(rho12-rho13*rho23)+x1*x3*(rho13-rho12*rho23)+x2*x3*(rho23-rho12*rho13))\n",
" d = 2*(rho12*rho12+rho13*rho13+rho23*rho23-2*rho12*rho13*rho23-1)\n",
" \n",
" fcn = -w/d\n",
" chisq = -2*fcn\n",
" return chisq*10000.\n",
"\n",
"constraint4 = triGauss(bplus_0, bplus_1, bplus_2)\n",
"\n",
"# 5. Constraint - Abs. of sum of light contribs\n",
"\n",
"sum_list_5 = []\n",
"\n",
"sum_list_5.append(rho_s*rho_w/rho_m)\n",
"sum_list_5.append(omega_s*omega_w/omega_m)\n",
"sum_list_5.append(phi_s*phi_w/phi_m)\n",
"\n",
"\n",
"sum_ru_5 = ztf.constant(0.0)\n",
"\n",
"for part in sum_list_5:\n",
" sum_ru_5 += part\n",
"\n",
"constraint5 = tf.cond(tf.greater_equal(tf.abs(sum_ru_5), 0.02), lambda: 100000., lambda: 0.)\n",
"\n",
"# 6. Constraint on phases of Jpsi and Psi2s for cut out fit\n",
"\n",
"constraint6_0 = zfit.constraint.GaussianConstraint(params = jpsi_p, mu = ztf.constant(pdg[\"jpsi_phase_unc\"]),\n",
" sigma = ztf.constant(jpsi_phase))\n",
"\n",
"constraint6_1 = zfit.constraint.GaussianConstraint(params = psi2s_p, mu = ztf.constant(pdg[\"psi2s_phase_unc\"]), \n",
" sigma = ztf.constant(psi2s_phase))\n",
"\n",
"#List of all constraints\n",
"\n",
"constraints = [constraint1, constraint2, constraint3_0, constraint3_1, constraint4, constraint6_0, constraint6_1]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Analysis"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"WARNING:tensorflow:From C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\core\\sample.py:163: to_int32 (from tensorflow.python.ops.math_ops) is deprecated and will be removed in a future version.\n",
"Instructions for updating:\n",
"Use tf.cast instead.\n",
"Toy 0: Generating data...\n"
]
},
{
"ename": "InvalidArgumentError",
"evalue": "Incompatible shapes: [799999] vs. [800000]\n\t [[node create_sampler/while/truediv_1 (defined at C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\core\\sample.py:270) ]]\n\nCaused by op 'create_sampler/while/truediv_1', defined at:\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\runpy.py\", line 193, in _run_module_as_main\n \"__main__\", mod_spec)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\runpy.py\", line 85, in _run_code\n exec(code, run_globals)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\ipykernel_launcher.py\", line 16, in <module>\n app.launch_new_instance()\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\traitlets\\config\\application.py\", line 658, in launch_instance\n app.start()\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\ipykernel\\kernelapp.py\", line 505, in start\n self.io_loop.start()\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tornado\\platform\\asyncio.py\", line 148, in start\n self.asyncio_loop.run_forever()\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\asyncio\\base_events.py\", line 539, in run_forever\n self._run_once()\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\asyncio\\base_events.py\", line 1775, in _run_once\n handle._run()\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\asyncio\\events.py\", line 88, in _run\n self._context.run(self._callback, *self._args)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tornado\\ioloop.py\", line 690, in <lambda>\n lambda f: self._run_callback(functools.partial(callback, future))\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tornado\\ioloop.py\", line 743, in _run_callback\n ret = callback()\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tornado\\gen.py\", line 781, in inner\n self.run()\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tornado\\gen.py\", line 742, in run\n yielded = self.gen.send(value)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\ipykernel\\kernelbase.py\", line 378, in dispatch_queue\n yield self.process_one()\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tornado\\gen.py\", line 225, in wrapper\n runner = Runner(result, future, yielded)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tornado\\gen.py\", line 708, in __init__\n self.run()\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tornado\\gen.py\", line 742, in run\n yielded = self.gen.send(value)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\ipykernel\\kernelbase.py\", line 365, in process_one\n yield gen.maybe_future(dispatch(*args))\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tornado\\gen.py\", line 209, in wrapper\n yielded = next(result)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\ipykernel\\kernelbase.py\", line 272, in dispatch_shell\n yield gen.maybe_future(handler(stream, idents, msg))\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tornado\\gen.py\", line 209, in wrapper\n yielded = next(result)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\ipykernel\\kernelbase.py\", line 542, in execute_request\n user_expressions, allow_stdin,\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tornado\\gen.py\", line 209, in wrapper\n yielded = next(result)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\ipykernel\\ipkernel.py\", line 294, in do_execute\n res = shell.run_cell(code, store_history=store_history, silent=silent)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\ipykernel\\zmqshell.py\", line 536, in run_cell\n return super(ZMQInteractiveShell, self).run_cell(*args, **kwargs)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\IPython\\core\\interactiveshell.py\", line 2848, in run_cell\n raw_cell, store_history, silent, shell_futures)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\IPython\\core\\interactiveshell.py\", line 2874, in _run_cell\n return runner(coro)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\IPython\\core\\async_helpers.py\", line 67, in _pseudo_sync_runner\n coro.send(None)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\IPython\\core\\interactiveshell.py\", line 3049, in run_cell_async\n interactivity=interactivity, compiler=compiler, result=result)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\IPython\\core\\interactiveshell.py\", line 3214, in run_ast_nodes\n if (yield from self.run_code(code, result)):\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\IPython\\core\\interactiveshell.py\", line 3296, in run_code\n exec(code_obj, self.user_global_ns, self.user_ns)\n File \"<ipython-input-37-cb6da38b83d8>\", line 17, in <module>\n sampler = total_f.create_sampler(n=event_stack)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\core\\basemodel.py\", line 821, in create_sampler\n limits=limits, n=n, name=name)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\core\\basemodel.py\", line 841, in _create_sampler_tensor\n sample = self._single_hook_sample(n=n, limits=limits, name=name)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\core\\basemodel.py\", line 884, in _single_hook_sample\n return self._hook_sample(n=n, limits=limits, name=name)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\core\\basepdf.py\", line 491, in _hook_sample\n samples = super()._hook_sample(limits=limits, n=n, name=name)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\core\\basemodel.py\", line 887, in _hook_sample\n return self._norm_sample(n=n, limits=limits, name=name)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\core\\basemodel.py\", line 891, in _norm_sample\n return self._limits_sample(n=n, limits=limits, name=name)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\core\\basemodel.py\", line 895, in _limits_sample\n return self._call_sample(n=n, limits=limits, name=name)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\core\\basemodel.py\", line 905, in _call_sample\n return self._fallback_sample(n=n, limits=limits)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\core\\basemodel.py\", line 938, in _fallback_sample\n sample_and_weights_factory=self._sample_and_weights)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\core\\sample.py\", line 355, in accept_reject_sample\n back_prop=False)[1] # backprop not needed here\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\ops\\control_flow_ops.py\", line 3556, in while_loop\n return_same_structure)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\ops\\control_flow_ops.py\", line 3087, in BuildLoop\n pred, body, original_loop_vars, loop_vars, shape_invariants)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\ops\\control_flow_ops.py\", line 3022, in _BuildLoop\n body_result = body(*packed_vars_for_body)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\core\\sample.py\", line 270, in sample_body\n prob_weights_ratio = probabilities / weights_clipped\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\ops\\math_ops.py\", line 812, in binary_op_wrapper\n return func(x, y, name=name)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\ops\\math_ops.py\", line 920, in _truediv_python3\n return gen_math_ops.real_div(x, y, name=name)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\ops\\gen_math_ops.py\", line 6896, in real_div\n \"RealDiv\", x=x, y=y, name=name)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\framework\\op_def_library.py\", line 788, in _apply_op_helper\n op_def=op_def)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\util\\deprecation.py\", line 507, in new_func\n return func(*args, **kwargs)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\framework\\ops.py\", line 3300, in create_op\n op_def=op_def)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\framework\\ops.py\", line 1801, in __init__\n self._traceback = tf_stack.extract_stack()\n\nInvalidArgumentError (see above for traceback): Incompatible shapes: [799999] vs. [800000]\n\t [[node create_sampler/while/truediv_1 (defined at C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\core\\sample.py:270) ]]\n",
"output_type": "error",
"traceback": [
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[1;31mInvalidArgumentError\u001b[0m Traceback (most recent call last)",
"\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 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[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 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[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 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[1;31mInvalidArgumentError\u001b[0m: Incompatible shapes: [799999] vs. [800000]\n\t [[{{node create_sampler/while/truediv_1}}]]",
"\nDuring handling of the above exception, another exception occurred:\n",
"\u001b[1;31mInvalidArgumentError\u001b[0m Traceback (most recent call last)",
"\u001b[1;32m<ipython-input-37-cb6da38b83d8>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m 31\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mcall\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mcalls\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 32\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 33\u001b[1;33m \u001b[0msampler\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mresample\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mn\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mevent_stack\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 34\u001b[0m \u001b[0ms\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0msampler\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0munstack_x\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 35\u001b[0m \u001b[0msam\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[0ms\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\\core\\data.py\u001b[0m in \u001b[0;36mresample\u001b[1;34m(self, param_values, n)\u001b[0m\n\u001b[0;32m 640\u001b[0m \u001b[1;32mraise\u001b[0m \u001b[0mRuntimeError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"Cannot set a new `n` if not a Tensor-like object was given\"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 641\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mn_samples\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mload\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mvalue\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mn\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0msession\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msess\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 642\u001b[1;33m \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[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msample_holder\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0minitializer\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 643\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_initial_resampled\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;32mTrue\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 644\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 1346\u001b[0m \u001b[1;32mpass\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 1347\u001b[0m \u001b[0mmessage\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0merror_interpolation\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0minterpolate\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmessage\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_graph\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1348\u001b[1;33m \u001b[1;32mraise\u001b[0m \u001b[0mtype\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0me\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnode_def\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mop\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[0m\u001b[0;32m 1349\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 1350\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m_extend_graph\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\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;31mInvalidArgumentError\u001b[0m: Incompatible shapes: [799999] vs. [800000]\n\t [[node create_sampler/while/truediv_1 (defined at C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\core\\sample.py:270) ]]\n\nCaused by op 'create_sampler/while/truediv_1', defined at:\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\runpy.py\", line 193, in _run_module_as_main\n \"__main__\", mod_spec)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\runpy.py\", line 85, in _run_code\n exec(code, run_globals)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\ipykernel_launcher.py\", line 16, in <module>\n app.launch_new_instance()\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\traitlets\\config\\application.py\", line 658, in launch_instance\n app.start()\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\ipykernel\\kernelapp.py\", line 505, in start\n self.io_loop.start()\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tornado\\platform\\asyncio.py\", line 148, in start\n self.asyncio_loop.run_forever()\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\asyncio\\base_events.py\", line 539, in run_forever\n self._run_once()\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\asyncio\\base_events.py\", line 1775, in _run_once\n handle._run()\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\asyncio\\events.py\", line 88, in _run\n self._context.run(self._callback, *self._args)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tornado\\ioloop.py\", line 690, in <lambda>\n lambda f: self._run_callback(functools.partial(callback, future))\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tornado\\ioloop.py\", line 743, in _run_callback\n ret = callback()\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tornado\\gen.py\", line 781, in inner\n self.run()\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tornado\\gen.py\", line 742, in run\n yielded = self.gen.send(value)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\ipykernel\\kernelbase.py\", line 378, in dispatch_queue\n yield self.process_one()\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tornado\\gen.py\", line 225, in wrapper\n runner = Runner(result, future, yielded)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tornado\\gen.py\", line 708, in __init__\n self.run()\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tornado\\gen.py\", line 742, in run\n yielded = self.gen.send(value)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\ipykernel\\kernelbase.py\", line 365, in process_one\n yield gen.maybe_future(dispatch(*args))\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tornado\\gen.py\", line 209, in wrapper\n yielded = next(result)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\ipykernel\\kernelbase.py\", line 272, in dispatch_shell\n yield gen.maybe_future(handler(stream, idents, msg))\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tornado\\gen.py\", line 209, in wrapper\n yielded = next(result)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\ipykernel\\kernelbase.py\", line 542, in execute_request\n user_expressions, allow_stdin,\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tornado\\gen.py\", line 209, in wrapper\n yielded = next(result)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\ipykernel\\ipkernel.py\", line 294, in do_execute\n res = shell.run_cell(code, store_history=store_history, silent=silent)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\ipykernel\\zmqshell.py\", line 536, in run_cell\n return super(ZMQInteractiveShell, self).run_cell(*args, **kwargs)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\IPython\\core\\interactiveshell.py\", line 2848, in run_cell\n raw_cell, store_history, silent, shell_futures)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\IPython\\core\\interactiveshell.py\", line 2874, in _run_cell\n return runner(coro)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\IPython\\core\\async_helpers.py\", line 67, in _pseudo_sync_runner\n coro.send(None)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\IPython\\core\\interactiveshell.py\", line 3049, in run_cell_async\n interactivity=interactivity, compiler=compiler, result=result)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\IPython\\core\\interactiveshell.py\", line 3214, in run_ast_nodes\n if (yield from self.run_code(code, result)):\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\IPython\\core\\interactiveshell.py\", line 3296, in run_code\n exec(code_obj, self.user_global_ns, self.user_ns)\n File \"<ipython-input-37-cb6da38b83d8>\", line 17, in <module>\n sampler = total_f.create_sampler(n=event_stack)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\core\\basemodel.py\", line 821, in create_sampler\n limits=limits, n=n, name=name)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\core\\basemodel.py\", line 841, in _create_sampler_tensor\n sample = self._single_hook_sample(n=n, limits=limits, name=name)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\core\\basemodel.py\", line 884, in _single_hook_sample\n return self._hook_sample(n=n, limits=limits, name=name)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\core\\basepdf.py\", line 491, in _hook_sample\n samples = super()._hook_sample(limits=limits, n=n, name=name)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\core\\basemodel.py\", line 887, in _hook_sample\n return self._norm_sample(n=n, limits=limits, name=name)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\core\\basemodel.py\", line 891, in _norm_sample\n return self._limits_sample(n=n, limits=limits, name=name)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\core\\basemodel.py\", line 895, in _limits_sample\n return self._call_sample(n=n, limits=limits, name=name)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\core\\basemodel.py\", line 905, in _call_sample\n return self._fallback_sample(n=n, limits=limits)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\core\\basemodel.py\", line 938, in _fallback_sample\n sample_and_weights_factory=self._sample_and_weights)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\core\\sample.py\", line 355, in accept_reject_sample\n back_prop=False)[1] # backprop not needed here\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\ops\\control_flow_ops.py\", line 3556, in while_loop\n return_same_structure)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\ops\\control_flow_ops.py\", line 3087, in BuildLoop\n pred, body, original_loop_vars, loop_vars, shape_invariants)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\ops\\control_flow_ops.py\", line 3022, in _BuildLoop\n body_result = body(*packed_vars_for_body)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\core\\sample.py\", line 270, in sample_body\n prob_weights_ratio = probabilities / weights_clipped\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\ops\\math_ops.py\", line 812, in binary_op_wrapper\n return func(x, y, name=name)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\ops\\math_ops.py\", line 920, in _truediv_python3\n return gen_math_ops.real_div(x, y, name=name)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\ops\\gen_math_ops.py\", line 6896, in real_div\n \"RealDiv\", x=x, y=y, name=name)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\framework\\op_def_library.py\", line 788, in _apply_op_helper\n op_def=op_def)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\util\\deprecation.py\", line 507, in new_func\n return func(*args, **kwargs)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\framework\\ops.py\", line 3300, in create_op\n op_def=op_def)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\framework\\ops.py\", line 1801, in __init__\n self._traceback = tf_stack.extract_stack()\n\nInvalidArgumentError (see above for traceback): Incompatible shapes: [799999] vs. [800000]\n\t [[node create_sampler/while/truediv_1 (defined at C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\core\\sample.py:270) ]]\n"
]
}
],
"source": [
"# zfit.run.numeric_checks = False \n",
"\n",
"Ctt_list = []\n",
"Ctt_error_list = []\n",
"\n",
"nr_of_toys = 2\n",
"nevents = int(pdg[\"number_of_decays\"])\n",
"nevents = pdg[\"number_of_decays\"]\n",
"event_stack = 1000000\n",
"# zfit.settings.set_verbosity(10)\n",
"calls = int(nevents/event_stack + 1)\n",
"\n",
"total_samp = []\n",
"\n",
"start = time.time()\n",
"\n",
"sampler = total_f.create_sampler(n=event_stack)\n",
"\n",
"for toy in range(nr_of_toys):\n",
" \n",
" ### 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 = (x_min, x_max), constraints = constraints)\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",
" #plotting the result\n",
" \n",
" plotdirName = 'data/plots'.format(toy)\n",
" \n",
" if not os.path.exists(plotdirName):\n",
" os.mkdir(plotdirName)\n",
"# print(\"Directory \" , dirName , \" Created \")\n",
" calcs_test = zfit.run(probs)\n",
" res_y = zfit.run(jpsi_res(test_q))\n",
" plt.clf()\n",
" plt.plot(test_q, calcs_test, label = 'pdf')\n",
" plt.legend()\n",
" plt.ylim(0.0, 6e-6)\n",
" plt.savefig(plotdirName + '/toy_fit{}.png'.format(toy))\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": []
},
{
"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
}
