{
"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",
"epsilon = 0.3\n",
"\n",
"# # Full spectrum\n",
"\n",
"obs_toy = zfit.Space('q', limits = (x_min, x_max))\n",
"\n",
"# Jpsi and Psi2s cut out\n",
"\n",
"obs1 = zfit.Space('q', limits = (x_min + epsilon, jpsi_mass - 50. - epsilon))\n",
"obs2 = zfit.Space('q', limits = (jpsi_mass + 50. + epsilon, psi2s_mass - 50. - epsilon))\n",
"obs3 = zfit.Space('q', limits = (psi2s_mass + 50. + epsilon, x_max - epsilon))\n",
"\n",
"obs_fit = 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_toy, 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",
"total_f_fit = total_pdf(obs=obs_fit, 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",
"# 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+17YcXAAAgAElEQVR4nO29eXzc1XX3/z4zo92WZG3eZCPZFjYyGBMcGzBhsQmYbOZpoTF5SElCS9qGX54kXYAnTfI0CW1o0tCmIQuFJCQlGELS4hIHQ7CBsNkYDBgvsuVdeNEuWfss9/fH9440Gs8qS5rR6LxfL780c+d+z73ztTSfOeeee64YY1AURVGUdMCV6gkoiqIoShAVJUVRFCVtUFFSFEVR0gYVJUVRFCVtUFFSFEVR0gZPqieQbpSVlZmqqqpUT0NR0p59p06Tm+Vmbkn+sPaOXi9HW3uoqZhCbpY7RbNz2N/YRbbbxTml+fgChj0nOplVnEdpQXZK55WJvPHGG83GmPKztaOiFEZVVRXbt29P9TQUJe25+jvPc8HsIr5380XD2p9+9yR/8Z9v8NjnL2fxrKIUzc5hzb++yJySfP7jT5fR1j3ARd94lq9+tJZPr6xO6bwyERE5Mhp2NHynKMqI8AcMLjmz3W0bA4FxnlAUglN02Xn5A7o3M51RUVIUZUQEjBn8oA/FE/zwT4ON+caA2CkOimUazEuJjoqSoigjIhAwuORMURrySFLvKhkMYn0ltwTnlcoZKfHQNSVFUUZEwAx90IeSTh/+oZ6Sy34FH29Pyev10tDQQF9f37iOO1bk5uZSWVlJVlbWmNhXUVIUZUT4jRn8oA/FnUZrN4aQ8J2kZl4NDQ1MnTqVqqoqJIKITySMMbS0tNDQ0EB19dgki2j4TlGUEWFM5PBdWomSMYNCkKp59fX1UVpaOuEFCUBEKC0tHVOvT0VJUZQRETBEESXnZ9okOtjHIoJIahIdMkGQgoz1e1FRUhRlRERPCXc+VgLp4Ckx/EPULZIWHpwSHRUlRVFGRLSU8ODajS8NPvyNMYTO0OWStPDg0pHnn3+ej3zkIwD09/dzzTXXsHTpUh577LFxnYcmOiiKMiKip4Q7P9PBIwlNdABHMNPBg0t3duzYgdfr5a233hr3sRPylERkjYjUiUi9iNwV4fUcEXnMvr5VRKpCXrvbtteJyHXxbIpItbWx39rMjjWGiJSKyBYR6RKR70eZ/wYReTexW6IoSiIEzFDyQCjptEk1dE0JnLmlQ6r6eHP48GEWLVrErbfeypIlS7jxxhvp6enh6aefZtGiRVx++eX85je/AaCxsZFbbrmFt956i6VLl3LgwIFxnWtcT0lE3MD9wAeBBuB1EdlgjNkd0u02oM0Ys0BE1gH3Ah8XkVpgHbAYmAX8XkTOtddEs3kvcJ8xZr2I/Mja/mG0MYA+4CvA+fZf+Pz/COhK6q4oihIXvzFEWvMOVnRIi/AdZtiakitFiQ5B/uF/drH7eOeo2qydVcjXPro4br+6ujoeeughVq5cyWc+8xm++93v8uMf/5jNmzezYMECPv7xjwNQUVHBgw8+yHe+8x2eeuqpUZ1rIiTiKS0H6o0xB40xA8B6YG1Yn7XAw/bxE8BqcX4T1gLrjTH9xphDQL21F9GmvWaVtYG1eUOsMYwx3caYl3DEaRgiMgX4EvDNBN6noihJYIyJuHk2GNJLhzBZZE8p9fNKBXPmzGHlypUA3HLLLWzfvp3q6mpqamoQEW655ZYUz9AhkTWl2cCxkOcNwIpofYwxPhHpAEpt+2th1862jyPZLAXajTG+CP2jjdEcY+7fAP4F6In1BkXkduB2gLlz58bqqiiKxR9lTSm99imFZd+lONEhEY9mrAhP5e7o6EjLVPVEPKVIsw7/X43WZ7TaE53H0IRElgILjDH/Fa3PoBFjHjDGLDPGLCsvP+vjQBQl4zHG2H1KZ77mTquCrMNDjG6X4Penfl6p4OjRo7z66qsAPProo1xzzTUcOnRocM3o0UcfTeX0BklElBqAOSHPK4Hj0fqIiAcoAlpjXButvRkotjbCx4o2RjQuBS4WkcPAS8C5IvJ8zHeqKEpCBJ0gd4Q6Q2nlKREWvpPJmxJ+3nnn8fDDD7NkyRJaW1v54he/yAMPPMCHP/xhLr/8cs4555xUTxFILHz3OlAjItXAeziJC58I67MBuBV4FbgR2GyMMSKyAfiliHwXJ9GhBtiG83tyhk17zRZrY721+WSsMaJN2hjzQ5wECWym3lPGmKsSeL+KosQhKDjuSLXvUlRjLhKhBVnB2aeUDmtdqcDlcvGjH/1oWNuaNWvYu3fvGX2vuuoqrrrqqnGa2XDiipJdv7kD2AS4gZ8YY3aJyNeB7caYDcBDwC9EpB7He1lnr90lIo8DuwEf8DljjB8gkk075J3AehH5JrDD2ibaGNbWYaAQyBaRG4Brw7IDFUUZRYIZbLE8pbRICQ85ugJSv6akxCehzbPGmI3AxrC2r4Y87gNuinLtPcA9idi07QdxsvPC22ONURVn/oeJkC6uKMrI8MXylIIp4WmwdhPuKU3WMkNVVVW8++7E2KqpZYYURUma4Ad7rEP+0sNTihC+S8G8Yqw0TDjG+r2oKCmKkjRBUfLEOg49DTySM1LCU+Ap5ebm0tLSkhHCFDxPKTc3d8zG0Np3iqIkzVCiQ/TNs2lR0SFSQdZxLjNUWVlJQ0MDTU1N4zvwGBE8eXasUFFSFCVpEkp0SAdRImxNyTX+YcWsrKwxO6U1E9HwnaIoSRMz0UHSbPMsqQ3fKcmhoqQoStIE4iQ6iKTJmhLpkeigJI6KkqIoSRP0lDzuyLXT0sUjOaMga5rMS4mOipKiKEkTKyUc0ueEV6f2XcjRFZO4SvhEQUVJUZSkCYbAPBESHZz29CjnEz4Dt2j4Lt1RUVIUJWmC1RoiJTqA8+GfDinhmOHe3GQ+T2mioKKkKErSxEoJh/QpfBoIO7rCCSumbj5KfFSUFEVJmlgp4eCE79JiTYnhiQ4el+APOLtnO/u8vLQ/1hmhSipQUVIUJWkSSnRIA0/pjIKsLhkMPf7dr97hloe20nS6P0WzUyKhoqQoStLES3RIl9Rrw/Dsuyz30FrXW8faAejo9aZkbkpkVJQURUmaoLcRRZMcjyQdRMmEh+9c+Gzxu9wsZ/Kn+1SU0gkVJUVRkiZuSrg7XTwlhqmSxy14BwXVeaHfN84VWpWYqCgpipI0iSQ6pMMhf+Ep4VkuFz6b6BBsHVBRSitUlBRFSZpYte8AstwuvON9RkQEAmFHV2R5hsQyuNakopReqCgpipI0Q4f8Rf4ISRdRCi/I6nENzSt4FNRAGsxTGUJFSVGUpAmG76IlOnjc6ZLoMPzoiqyQNaVgu3pK6UVCoiQia0SkTkTqReSuCK/niMhj9vWtIlIV8trdtr1ORK6LZ1NEqq2N/dZmdqwxRKRURLaISJeIfD/ETr6I/FZE9orILhH5VvK3R1GUSMRLdMhypamn5A5ZUwp6SipKaUVcURIRN3A/cD1QC9wsIrVh3W4D2owxC4D7gHvttbXAOmAxsAb4gYi449i8F7jPGFMDtFnbUccA+oCvAH8TYfrfMcYsAi4CVorI9fHer6Io8Ymb6OBOj0SH8JTwLJfjKYVWD+9PA/FUhkjEU1oO1BtjDhpjBoD1wNqwPmuBh+3jJ4DV4vyPrwXWG2P6jTGHgHprL6JNe80qawNr84ZYYxhjuo0xL+GI0yDGmB5jzBb7eAB4Exi7g+UVZRIRL9HB43bhTYPwHTDMVfJYFfUHjGbfpSmJiNJs4FjI8wbbFrGPMcYHdAClMa6N1l4KtFsb4WNFGyMuIlIMfBR4Lsrrt4vIdhHZ3tTUlIhJRZnUxE10cMngJtVUYUxw7WiILCtKvoAZXA/r9/nHeWZKLBIRpUhfhcK/AkXrM1rtic7jDETEAzwKfM8YczBSH2PMA8aYZcaYZeXl5fFMKsqkx59IokOKw3fBerCusDJDAF5/YDDRwetLE49OARITpQZgTsjzSuB4tD5WBIqA1hjXRmtvBoqtjfCxoo0RjweA/caYf02gr6IoCeCPW9HBhTeQWk8pmIwxPCXceeLzm0FhHfCrp5ROJCJKrwM1NisuGydxYUNYnw3ArfbxjcBm4/jOG4B1NnOuGqgBtkWzaa/ZYm1gbT4ZZ4yoiMg3ccTrCwm8T0VREiReSnhWGlR0iBRiCa4pef2BwexAXVNKLzzxOhhjfCJyB7AJcAM/McbsEpGvA9uNMRuAh4BfiEg9jveyzl67S0QeB3YDPuBzxhg/QCSbdsg7gfVWUHZY20Qbw9o6DBQC2SJyA3At0Al8GdgLvGkzbb5vjHkw+dukKEoowUQHd4xEh9SvKTk/Q6c4GL4LmEFh9aZBlqAyRFxRAjDGbAQ2hrV9NeRxH3BTlGvvAe5JxKZtP4iTnRfeHmuMqihTj/wXoyjKWRG/ooMwkHJPaXg5IRiar88fGPSQtCBreqEVHRRFSZp4iQ5ZIZtUU0Wk4H6WJxi+Mxq+S1NUlBRFSZq4iQ4uV8rXlIIMC98FEx0CgcHwnda+Sy9UlBRFSZr4npKkvMzQ4JoSZ26e9fkNXushedVTSitUlBRFSZqg4GTFOOQv1QVZg2tKrtCU8JB9SkEPST2l9EJFSVGUpPH5DS4ZOr01HI/LhT9giLNrY0wJRMq+cw2tKQ2G79RTSitUlBRFSRpvIDAYCovEUOWE1IlSUBAjVXQY8AVCNs+qKKUTKkqKoiSNz28GkwYiMbh2k8IMvCFP6cw1pV7vUBUH9ZTSCxUlRVGSxueP7SkFy/mkh6c01Bb0lFSU0hcVJUVRksYbMIMf8JEYrMadwtBY0FNyRdg82zvgG2zT8F16oaKkKErS+PyBqHuUIDTLLZWJDmd6SjlZzpy7+4c8pVSnrivDUVFSYtLR4x1cEFaUID6/GRSeSGSFFD5NFUNVwofmmWMrOnT1h3hKGr5LK1SUlKj0DPi48OvP8I2ndqd6Kkqa4YTv4mffpXKvUqTzlHKz3MCQKOVmuVSU0gwVJSUqwRDHU++EH5+lTHac8F2M7DtXOqwpRQjfWU/pdJ8jSgXZHhWlNENFSVGUpPH6TdrvUxoshRTBU+q2nlJ+jlsTHdIMFSVFUZLGFwjEzL4b9JRSuE8p0nlKHpfgkqHwXUG2hwF/IKWVJ5ThqCgpipI0Pr+JGb4LHhGRytBYMHznDpmniJCb5R4UpbxsN8akdu1LGY6KkqIoSeONs3k2Jy1EyfnpCjsdN8fjoitkTQk0LTydUFFSFCVpvP7Y4bugKPWnQaJD+IntOZ4hTyk/21lj0mSH9EFFSVGUpPEFTMzNs9lBUfKmck3pzEQHcNLAg6I0JcfxlFSU0gcVJUVRksbrj71PKcfjeCD9Pn/UPmNN9PDd8DUlgH4VpbQhIVESkTUiUici9SJyV4TXc0TkMfv6VhGpCnntbtteJyLXxbMpItXWxn5rMzvWGCJSKiJbRKRLRL4fNq+LRWSnveZ7IuGOvKIoI8GXYPguHRIdwvMxQjfMFuTomlK6EVeURMQN3A9cD9QCN4tIbVi324A2Y8wC4D7gXnttLbAOWAysAX4gIu44Nu8F7jPG1ABt1nbUMYA+4CvA30SY/g+B24Ea+29NvPerKEp8fIHY+5SCNeZS6YEEs9HDv4sGvTgIWVNSUUobEvGUlgP1xpiDxpgBYD2wNqzPWuBh+/gJYLX1StYC640x/caYQ0C9tRfRpr1mlbWBtXlDrDGMMd3GmJdwxGkQEZkJFBpjXjVOcPnnIbYURTkLvP5AzPOUctypD4tF85SCggma6JCOJCJKs4FjIc8bbFvEPsYYH9ABlMa4Nlp7KdBubYSPFW2MWPNuiDNvAETkdhHZLiLbm5qaYphUFAXiF2Qd8pRSt6YUqfYdDPeUpuZmARq+SycSEaVIv3nhO82i9Rmt9kTnkciczmw05gFjzDJjzLLy8vIYJhVFAadSQ6zwXbY7jdaUwqYZTG4AKLSiNB4e3Yv7mvj0T7fRO5A6oZ4IJCJKDcCckOeVQHiFzsE+IuIBioDWGNdGa28Giq2N8LGijRFr3pVx5q0oygjwxjkO3eUSstySFuG78DWlqbmewceFeeOXEv7tTXVsqWti2+FYH1tKIqL0OlBjs+KycRIXNoT12QDcah/fCGy26zgbgHU2c64aJ9lgWzSb9pot1gbW5pNxxoiIMeYEcFpELrFrVX8aYktRlLMg3nHo4ITJUrlPKVpKeNA7gqHw3XiI0omOXgAONnWN+VgTGU+8DsYYn4jcAWwC3MBPjDG7ROTrwHZjzAbgIeAXIlKP472ss9fuEpHHgd2AD/icMcYPEMmmHfJOYL2IfBPYYW0TbQxr6zBQCGSLyA3AtcaY3cBfAj8D8oDf2X+Kopwl3kDsNSVw0sIH/KlcU4qc6BDqKeXZquFjXc3cGDN4XMbx9t4xHWuiE1eUAIwxG4GNYW1fDXncB9wU5dp7gHsSsWnbD+Jk54W3xxqjKkr7duD8SK8p8TExl+yUyYzPHyArRkUHcEQpLT2lvCFPKVh5YqzFs9frHwxlHm/vi9N7cqMVHZTo2D9qreqvhOLzBwiYoQ/0aGR7XGmypjS8vTDEUxqvckgdvd7Bxy3d/WM61kRHRUmJi9bBUEIJCk1OHFHK8bhTW2YowiF/MDx8l2/Dd73esZ1nZ69v8HFbtzdGT0VFSYmKOkhKJIJJAfE8pZyQcj6pIBi+c4ctKpUU5Aw+zs9xRKlnjNO0g55S5bQ8WnsGxnSsiY6KkhIVDdspkQiW5AndhBqJbHd6hO/CEx3mlxcAcE5pPtluFx6XDB6PPlYERam6rIC27gE96TYGCSU6KJMTTXRQIhFcf0nEU+pLaaJDtH1KWXxj7WKWzpmGiJCf7R43T6mqtIA/7G/mdL9vWGq6MoSKkqIoSRHMVEtkTSl0gX+8iVZmCOCTl1YNPs7P9tAzMLaeUmdQlMocL621a0BFKQoavlOiohEGJRL9ia4ppTwlPHL4Lpz8HDfd4+Yp5QPoulIMVJSUqAxpkqbfKUMkKkq5WW760vCQv3AKsj1jXo+uo9fL1BwPZVOcJIt2FaWoqCgpURlajFWXSRliIMGU8Lxsd0qLj0bbpxROfrZ7zBMdOvu8FOZlUVKQDUCrpoVHRUVJUZSkSHSfUn6Wm+7+dCgzFMdTyvGMeaJDZ6+XorwspllRautWTykaKkpKVIbWlDR8pwwx5CnFTgnPz/HQ6/UPbmIdbxIN3+Vlu+ke40SHDitKBdlustyia0oxUFFSFCUpEt08GzzVNVXrSokmOhSMQ5ixs9dHYZ4HEWFafrZ6SjFQUVIUJSmCpYOy4xxdERSlVIXwgp5S+D6lcPKzPXSNx5qSTQEvKcimVUUpKipKSlSGwnea6KAMMRi+y4onSs42yFQlO0Q7uiKcorwsuvp9+McwzHi6zzd4dtO0/GzaNHwXFRUlJSom5kn0ymQlWGYoUU+pxzu2Xkg0/FEKsoZTnJ+FMUMbXMdiHl39vsFTbksKsmnr0ey7aKgoKYqSFMENsTlZsRMd8lIcvguKUnhB1nCK8x0Ppn2MRKnLHu436CkVZOmaUgxUlJSoaEUHJRIJe0rBYyFSFL5LWJTybJr2GIXUOvscsQue4xQM36UqKzHdUVFSoqJ/Mkok+u3ZQ1lxjkMvyHE+hMe6rlw0/PZblSdBT6ljjEJqQVEKXVMKmKF2ZTgqSoqiJEW/P0COxxU3qy0YvhvrA/SikXj4zvGU2nvHyFOyB/yFrikBmoEXBRUlJSp65osSib4B/6DgxCLVKeE+f6LhO8eDGasTYYfCd8E1pbENF050EhIlEVkjInUiUi8id0V4PUdEHrOvbxWRqpDX7rbtdSJyXTybIlJtbey3NrPPYowvisguEXlXRB4Vkdzkbs/kZjD3TpPvlBB6BvyD60WxyM9KbfguuHk2nigV5mUhMnaJDqdtosPgPqX8YKkhDd9FIq4oiYgbuB+4HqgFbhaR2rButwFtxpgFwH3AvfbaWmAdsBhYA/xARNxxbN4L3GeMqQHarO2RjDEb+DywzBhzPuC2/ZQkUYdJCaXH6yc3AU9pMHyXokQHXyC4phT7Y87tEorzsmju6h+TeQRTzYPhu+AalpYaikwintJyoN4Yc9AYMwCsB9aG9VkLPGwfPwGsFifgvBZYb4zpN8YcAuqtvYg27TWrrA2szRtGOAY4hxjmiYgHyAeOJ/B+FYuKkRKJvgH/YGguFtkeF9ke15hXS4jG4D6lBD7lphfm0tjZNybzCHpKU3KGrylpWnhkEhGl2cCxkOcNti1iH2OMD+gASmNcG629FGi3NsLHSmoMY8x7wHeAo8AJoMMY80ykNygit4vIdhHZ3tTUFPVGTD4SK/2vTC56BvzkJRC+AydklaosM3+CnhLAjKJcTo6RKHX0einIduOxKfT52W6yPS71lKKQiChF+kgK/w4drc9otSc9hohMw/GiqoFZQIGI3BKhL8aYB4wxy4wxy8rLyyN1mdSox6SE0uv1k2dLCMWjMM8zmH023vgCiZUZAphRmMvJjrEJ37X1DAwmN4BTi69Ei7JGJRFRagDmhDyv5Mww2GAfGyorAlpjXButvRkotjbCx0p2jGuAQ8aYJmOMF/gNcFkC71exqBgpkegd8JMXp+5dkNR6SgHcLombug5O+K6lux+vf/SPb2/pHqA0RJTACeE1d6koRSKR36zXgRqbFZeNkyywIazPBuBW+/hGYLNx8ok3AOts5lw1UANsi2bTXrPF2sDafHKEYxwFLhGRfLv2tBrYk9htUUCz75TI9Hh9g8VW41GUlzVmNeXi4Q/Ez7wLMqMoF2Og8fToe0ut3f3DPCWAmUW5nOwYm3DhRCeuKNn1mzuATTgf6o8bY3aJyNdF5GO220NAqYjUA18C7rLX7gIeB3YDTwOfM8b4o9m0tu4EvmRtlVrbIxljK05CxJvATvteHxjBPVIUJYTegQC5ia4p5WXR2ZeqRIcA7gS/Uc0scnaLHG/vHfV5tHYNDCY3BBnLNayJTkJfd4wxG4GNYW1fDXncB9wU5dp7gHsSsWnbDzKUPRfaPpIxvgZ8LdI1Snw0fKdEonfAl1D2HTj13lLlKfkCJm6JoSDVZQUAHGru5v1VJaM2B2NMxPDdzKJcWrsH6PP6Exb4yYJWdFCiYrT6nRKGMcZJdEjKU/KmpDpIIGBwx6nPF2R2cR5ZbuFgU/eIx4pEz4Cffl+AkoKcYe0zivIAOKXe0hmoKClRUU9JCaffFyBgSKjMEDiJDl6/oc87+gkE8fAFTMLhO4/bxdySfA41dyU9zv5Tp1n69Wf4p41nLlm32GSG0ilnekoAJ3Rd6QxUlBRFSZg+W1w1cU/JWSFIRQaeP2ASTnQAmFc+hUPNyXtKv3jtCJ19Pn784kFawqpCHO9w1qhmWc8oyAwrSprscCYqSkpU1FNSwumxJYOS8ZRg7E51jYU/iTUlgAUVjij1+5Iri7TreCdT7VlJT+86Oey1YOLErOLhZTdnFNrEio7RT6yY6KgoKVEJrilpRrgSJFgyKFgyJx7BOm+pOP7bHzC4khClC2YX4fUb6k6eTvgaYwz1jV189MJZVJXm8/S70URpuKdUkOOhpCCbY609CY81WVBRUuKiDpMS5PTggXWJiVKpXeAPD2uNB8lk34EjSgBvN3QkfE1z1wAdvV4WlE/husUzePVACx0hXuF77b2UFGRHzLCrLisYUbgw01FRUqKi4TslnOCeo+ApqvEoswv8zSkoqeM3ya0pVU7Lo6Qgm50N7QlfU9/oJEYsqJjCtYtn4AsYnq9rHPb6/PKCiNdWlaooRUJFSYmLhu+UIF2DZwMl5ikFN402j0GlhHj4/cmJkoiwdE4x2w61JnxNfdOQKF00p5iKqTlssutKxjihwHOnT4147bzyAk519qfsvKl0RUVJiYs6TEqQwWMYEhQlj9vFtPwsWrpTE75zJ3JuRQiXLyjjcEtPwms9Bxq7KMh2M7MoF5dL+GDtdJ6va6LP66ehrZfOPh+LZkQWpapSx4M63KzrSqGoKClR0fCdEs7QmlJi4TuAsik5NJ9OQfguEEhqTQnginPLAPjD/uaE+tc3djG/Yspg0ddrF8+gZ8DPH/Y388oBx8aKeaURrw1WkTjQlPzeqExGRUmJimbfKeF09ftwCRQkmBIOzsbRsTrVNRZevyHbk9xH3PzyKcwqyuW5PacS6l/f2MWC8imDzy+dV8r0why+99x+Htl6lNnFedRUTIl47fyKAjwuYc+JzqTmmOmoKCmKkjCn+3xMyfEkdBxEkLIpObSkINFhwBcgK8EyQ0FEhA8vmckL+5rinnd0us/Lyc4+5oeITrbHxZc/XMvO9zp4p6GDv7hqftR7leNxUzN9KruOqyiFklhgWJmUaPhOCaezz5tU6A6C4bvx95QG/AEKs5ObK8ANF83mP/5wiN/uPMEtl5wTtV9o5l0oH7twFmVTsjnd5+Pa2ukxxzp/ViGb9zZijElK6DMZ9ZSUqKgmKeF09fkS3qMUZGZRLqf7fYPrUeOF1x8gO0lPCaB2ZiGLZkzlP187ErOQ7P5TjihFyq67bH4Z1y2eEVdoFs8qpKV7gFOdsUW7u9/Hd5+pmxRliVSUFEVJGMdTSk6UZk9zqhm8NwZnFcXCCd8l/xEnIvzZB+ax9+TpmAkP+06dJsfjFHIdKRdUFgPw1rG2mP1+/WYD39tczzd/u3vEY00UVJSUqKTiuAElvWnr9jItPzt+xxBm2xI7Da3jK0pe/8hECZwQ3PTCHP7tuf1R/w72NXYxv3xKUnuhwrlgdhF5WW5eOxh7b9SOo86G3nffS7zaxERFRUmJikqSEk5rz8AZxzDEI1We0kiy74Jke1x88ZpzeeNIG0+9c+KM1wMBw86GdhbPKjyrOWZ7XCyrmsZrB1ti9jtoKz8cbe0ZrNSeqagoKYqSEMYY2roHkvaUygpyyPa4xj98dxaeEsBNy+ZQO7OQbzy1m9awTLwDTV209Xh5f/XZn1J7ybxS9p48fcYYQYwxHGrqoigvi4CBwy2ZXZpIRUmJikbvlFA6+3z4AmawdFCiuFzC7OI83msb/zWlkcwu7xIAACAASURBVCQ6BHG7hG/ftIT2Hi9/+6u38YecLvt8XRMAl1RH3hibDJfNd2y8uK8p4uut3QN09vn4oM3kO9CooqRMWlSVlCGC+3aSFSWAOSX54/4N/2zWlIIsnlXElz98Hs/tbeTv//td/AGDP2D41RvHuGB2EXNLR57kEOTCyuE188IJFm29emEFoJ4SACKyRkTqRKReRO6K8HqOiDxmX98qIlUhr91t2+tE5Lp4NkWk2trYb21mn8UYxSLyhIjsFZE9InJpcrdHUZQgrT2OKE0bgSjVVEzhQFMXgcD4fdHx+gMjXlMK5dbLqvirq+bz6LajfPzHr/L59TvYd6qLP79i3ijMksGaeS/sa4q4XnS4xamNVzurkOmFORxsmuSiJCJu4H7geqAWuFlEasO63Qa0GWMWAPcB99pra4F1wGJgDfADEXHHsXkvcJ8xpgZos7aTHsNe82/A08aYRcCFwJ5Eb4yi4TtlOK1d1lNKck0JHFHq8wbGbV3JGIPXb87aUwryt9ct5Ns3LqGhrZdN757ks1fM46NLZo6KbYDrbM28YFgwlCMt3bhtCLSqtEA9JWA5UG+MOWiMGQDWA2vD+qwFHraPnwBWi7NrbC2w3hjTb4w5BNRbexFt2mtWWRtYmzeMZAwRKQSuAB4CMMYMGGMSPyhF0eCdMozWswjf1Ux3qh7sb0z8VNezYcAfABgVTwmcvUs3LZvDq3evYu831nD3h84b1QoMl80vpWJqDr/afuyM14609DCrOJdsj4vqsgIOZ/gZTIn8j80GQu9Ug22L2McY4wM6gNIY10ZrLwXarY3wsZIdYx7QBPxURHaIyIMiEvG0LRG5XUS2i8j2pqbIi42KMtlpPO1UEyifmpP0tQvKnaoHwSoIY82AzxGlZGvfxUNE8IyS9xWKx+3ixosr2VLXeEbVhiMt3ZxT4nx0VZcV0NI9MOx020wjkbsb6X81/Et0tD6j1T6SMTzA+4AfGmMuArqBM9bDAIwxDxhjlhljlpWXl0fqMinR8J0SyomOPqblZ0U82jseRflZzCzKHbfio31eR5TyRjDXVPHx988hYODRbUcH24wxHGzu5hybUFFVFjyDKXO9pUREqQGYE/K8EjgerY+IeIAioDXGtdHam4FiayN8rJGM0WCM2Wrbn8ARKSVBgjvZtU6kAnCqs48ZRXkjvn7pnGJ2xCmnM1r0DjgJA3nZE6fm9DmlBVxz3nQefvUwXf1OsKihrZfTfT4WzyoChs5gyuR1pURE6XWgxmbFZeMkFWwI67MBuNU+vhHYbJxPtA3AOps5Vw3UANui2bTXbLE2sDafHMkYxpiTwDERWWivWQ1kfuGoMUA9JgUcT2lGYfKhuyAXzS3mWGsvTeNQMbzXZrFNJE8J4I5VC2jv8fLzVw8D8E6DU1YoWDlibkk+IkNp4plIXFGy6zd3AJtwstceN8bsEpGvi8jHbLeHgFIRqQe+hA2TGWN2AY/jiMHTwOeMMf5oNq2tO4EvWVul1nbSY9hr/j/gERF5B1gK/GOyN2gyo1qkhHK2ntL75k4DYMfRsfeWegYcTyM/icMI04Glc4pZvaiC72+u51hrDy/sa2RqrodaK0q5WW5mFeVltCgl5NsaYzYCG8PavhryuA+4Kcq19wD3JGLTth/Eyc4Lbx/JGG8ByyJdo8Qn6CFp+E7p9/lp7hpgRmHuiG2cP7uIHI+LVw+2cO3iGaM4uzMZ9JQmmCgB/MPaxVx334t88qGtnOjo46MXzhqW2p7pGXha0UGJi4bvlBPtTkbYzOKRi1JulptL55eyZW/jaE0rKoNrShMsfAdQOS2fH39yGZ19PqYX5vKFa2qGvV5dVsCh5u6MreI/cVYBlXHHaABPsRyyC+vBhfaRcvXCCr5Wt4tDzd1nbSsWQU9pooXvglxeU8Ybf38NwBn7oarKCujs89HaPUDplJGv8aUr6ikp0dHwnWI5YsNFVaVnJySrFjn12zbuPPM4iNGkx3pKI0lfTxdEJOIG3eoyJz08UzPwVJQURYnL4ZYeCrLdlCV5llI4c0ryWVFdwq+2HxvT8FMwfDdRPaVYVJc51TEONfekeCZjg4qSEhUN3ilBDjV3U1VWMCqldW5aNofDLT28Gudgu7Oh01Y8mJqbNWZjpIrKaXlkuWXcSjaNNypKSlQydB1VGQEHmrpGbQ3oI0tmUjYlh/u31I+KvUh09HrJy3KPWu27dCLL7eK8mYW8cywzj0bPvP8xRVFGlY5eLw1tvZw38+yO/g6Sm+Xms1fM4+X6Fl49MDbeUkevl+L8zPOSgiypLOLd9zrG9SiQ8UJFSYmKZt8pAHtOOPXqglUFRoNbLjmHyml5fPm/d0Y8Q+hsae/1UpSXyaJUzOl+HwczcL+SipISFQ3fKQC7bRHV2lEUpbxsN//4vy7gYFM333hq9Kt/dfR6KcxgUbqwshiAt49l3mk8KkqKosTk3fc6KJ+aQ8XUkW+cjcQV55bz2Svn8cjWozzw4oFRtd3eM0BxBovSgoopFGS7x63A7Xiim2eVqKijpABsPdTKsnOmjYntv7tuEcdae/jHjXs53efjC9eci9t19hl+pzr7uWRe6SjMMD1xu4RlVSW8drA11VMZddRTUqKSqWVMlMQ51trDe+29rKguGRP7bpfwvXUXcdPFlfz75nr+9CdbOdZ6dvtvegZ8dPR6mVE0up5dunHJvFLqG7sGD1/MFFSUFEWJytZDzjfxFWPodXjcLv75xiV8648u4M0j7az+7gt8e9PewePXkyV4cuvMDBelS+c7/ydbM8xbUlFSoqJ+kvLcnlOUT81h4fSpYzqOiLBu+Vye++srWbN4BvdvOcBl33qOrz75LgeakjtC/UiL42nNLs4fi6mmDefPKmRKjmdMNyGnAhUlJTqqSpOaPq+fF/Y1cW3tdFyjsM6TCLOK8/jezRfx7Bev4KNLZvHotqOs/pcX+MR/vMbGnSfw+gNxbdSdciodjLWQphqP28Xy6hJeqW9O9VRGFRUlJSrBfUqCVmSdjLywr4meAT/XjfHZR5GomT6Vb990IS/ftYq/ufZcjrT08FePvMnKb23mu8/Ucby9N+q1777XwcyiXIoyePNskCtqyjjc0pNR5yupKClRCdgvpeP0JVlJMx57/RgVU3MG1y5SQcXUXO5YVcOLf3c1D926jMWzCvn3LfVcfu9m/vzn23lhX9OwqgY+f4CX65u5bH5ZyuY8nly10Km6/nzd2J9RNV5oSrgSlYDNvhuNIpzKxOK99l6er2vkc1cvGHbqaapwu4TV501n9XnTOdbawy+3HeXx14/x7O5TzC3J5xMr5vIny+awaddJ2nq8XH/++Ht3qaCqrIB5ZQVsqWviUyurUz2dUUFFSYlKBpbVUhLkxy8cwGWTD9KNOSX53LlmEV+4poZNu07xn68d4Vu/28s/P72XgIHlVSWD5zZNBq5eVMEvXjtC74B/Qh7/Ho6KkhKV4D4lV+q/KCvjyHvtvazfdoybls1hdnFeqqcTlRyPm49dOIuPXTiL/adO85sd7zElx8Otl1WNW2JGOnD1wgoeeukQrx5sZtWi6amezlmT0MeNiKwRkToRqReRuyK8niMij9nXt4pIVchrd9v2OhG5Lp5NEam2NvZbm9kjHcO+5haRHSLyVOK3RYEhT0kTHSYX3/if3YjAHasWpHoqCVMzfSp3rlnE565ewJScyfVd+/3V08jPdrN5b2asK8UVJRFxA/cD1wO1wM0iUhvW7TagzRizALgPuNdeWwusAxYDa4AfWJGIZfNe4D5jTA3QZm0nPUbI3P4PsCex26GEMrSmlOKJKOPGxp0neHrXSb5wzblp7SUpQ+R43KxcUMaWvU0ZUYUlEU9pOVBvjDlojBkA1gNrw/qsBR62j58AVouzOr4WWG+M6TfGHALqrb2INu01q6wNrM0bRjgGIlIJfBh4MLHboYQSFCWXqtKk4EBTF3c+8Q5LKov4sw9kxqL5ZGHVogrea+9l36nkNhqnI4mI0mzgWMjzBtsWsY8xxgd0AKUxro3WXgq0WxvhYyU7BsC/An8HxNxxJyK3i8h2Edne1NQUq+ukwgyG75RMp+l0P3/+8HayPC5+eMvFaZFxpyTO1TY1PBNCeIn85kX6TAr3EaP1Ga32pMcQkY8AjcaYNyK8PryzMQ8YY5YZY5aVl5fH6z5pGDzkT1Upo2nu6ud/P/gaJzr6+PEnL9aw3QRkRlEui2cVsnnvqVRP5axJRJQagDkhzyuB49H6iIgHKAJaY1wbrb0ZKLY2wsdKdoyVwMdE5DBOeHCViPxnAu9XsQQ3z6omZS57T3ay9vsvc7S1h4c+tYz3V41NNXBl7Fm9qII3jrTR3jOyQrbpQiKi9DpQY7PisnGSCjaE9dkA3Gof3whsNs6K2wZgnc2cqwZqgG3RbNprtlgbWJtPjmQMY8zdxphKY0yVtb/ZGHNLgvdFQTfPZjLGGJ54o4E//sEr+AIBHv/spZOmCkKmcvWiCgLGKQ81kYmbO2mM8YnIHcAmwA38xBizS0S+Dmw3xmwAHgJ+ISL1ON7LOnvtLhF5HNgN+IDPGWP8AJFs2iHvBNaLyDeBHdY2IxlDOTuCa0qTaMvHpKCxs4+//+93eWb3KZZXl/C9dRdl/NlDk4ELK4spLchm895G1i4NX/afOCSU0G+M2QhsDGv7asjjPuCmKNfeA9yTiE3bfhCbPRfWnvQYIa8/Dzwf7XUlMoOekgbwMoI+r5+HXjrED7bU4/Ubvvyh87jt8upJtdE0k3G5hKsWVvD7Pafw+QN4JmiyyuTaZaYkxeDmWf3MmtB4/QH+a8d7fO+5/TS09XJt7XT+74fOo6qsINVTU0aZ1edV8Os3G9hxrH3Crg+qKClR0TWliU2/z88TbzTww+cP0NDWy/mzC/nnP17CZQt07ShTubymDI9LeG5Po4qSknmYwfCdMpFoOt3Po9uO8sjWI5zq7GfpnGK+sfZ8rlpYrl8wMpzC3CyWV5ewZW8jd12/KNXTGREqSkpUNHw3sXj7WDsPv3KYp945wYA/wAdqyvj2jRfygZoyFaNJxKpFFXzzt3toaOuhctrEOxJeRUmJitHad2nP6T4vv33nBOtfP8Zbx9opyHZz8/I5/OllVcwvn5Lq6Skp4GorSlv2NvLJS6tSPZ2kUVFSoqJVwtMTYwzbDrXy+PYGNu48Qa/Xz4KKKfy/j9byxxdXMjU3848BV6Izr6yAqtJ8nlNRUjKNoYKsKZ6IAsDJjj5+/WYDv9p+jMMtPUzJ8XDDRbP5k2WVLJ1TrCE6BXASk65eVMEjW4/SM+AjP3tifcxPrNkq48pgFXz9sEsZp/u8PLPrFE++fZyX9jcRMLCiuoTPr65hzfkzJtwHjjI+rF40nZ++fJhX6lu4pnZiHfynv9FKVAKafZcS+n1+Xqhr4sm3j/P73afo9wWonJbHX121gBsvrtT9RUpclleXUJDtZnNdo4qSkjn47KKSR+N3Y04gYNh6qJUNb7/Hxp0n6ej1UlKQzcffP4e1S2fzvrkanlMSJ9vj4gM15WzZ24gxZkL97qgoKVHx+R1RcqsojQnGGHYd72TD28fZ8NZxTnb2kZ/t5rrFM/jY0llcvqBMzzVSRsyqRRU8veske06cpnZWYaqnkzAqSkpUfPbsCj15dvQICtHGnSfYuPMEh1t68LiEK88t5+4PLeKDtdN1nUgZFa5c6JwN91J9k4qSkhkEw3fmjDMdlWQICtFvd57gd1aI3C7h0nml/PkV87j+/JmUFGSneppKhjG9MJcFFVN4qb6F26+Yn+rpJIyKkhIVn9/xlIxqUtKECtHGnSc4YoXosvmlfPbK+VxbO53SKTmpnqaS4Vy+oIz1rx+l3+cnx+NO9XQSQkVJiYrXr2qUDMYY3n1vSIiOtg4J0V9eOZ9rF89Qj0gZV1YuKONnrxxmx9F2LplXmurpJISKkhIV/2D4TomGP2DYcbSNTbtOsmnXKY62OmtEly0o43NXz+fa2hlMUyFSUsSKeSW4XcLL9c0qSsrEJ5jooAyn3+fnlfoWntl9kmd3n6K5a4Bst4tL55dyx9UL+GDtdBUiJS0ozM3iwsoiXqpv5q+vXZjq6SSEipISlcHwnbpKnO7zsqWuiWd2neT5uia6+n1MyfFw1cJyrls8g6sWlmvNOSUtWbmgjPu31NPZ56VwAvyOqihlODuOtlE5LZ/yqckvqvsnefZd4+k+fr+7kU27TvLKgWa8fkPZlGw+euFMrl08g8vml06YxWNl8rJyQRn/vrme1w60cO3iGameTlxUlDKc//WDV6iYmsO2L1+T9LVe/+QL3x1q7ubZ3c760JtH2zAG5pbk86nLqrhu8QwumjtNNxMrE4qL5haTl+Xm5frmzBElEVkD/BvgBh40xnwr7PUc4OfAxUAL8HFjzGH72t3AbYAf+LwxZlMsmyJSDawHSoA3gU8aYwaSHUNE5tj+M4AA8IAx5t+SvUETmYD1dBpP94/o+n5f5qeE9/v8bDvUypa9TWypa+RQczcAi2cV8sVrzuXaxdNZOH3qhCrToiih5HjcLK8u4eUDLameSkLEFSURcQP3Ax8EGoDXRWSDMWZ3SLfbgDZjzAIRWQfcC3xcRGqBdcBiYBbwexE5114Tzea9wH3GmPUi8iNr+4cjGMMH/LUx5k0RmQq8ISLPhs07o+nx+s/q+t6Bs7s+XTnV2ceWvY1s3tvIy/XNdA/4yfa4uHReKZ+6rIpViyqYUzLxTuxUlGhcOr+Ub/1uL81d/ZSl+f64RDyl5UC9MeYggIisB9YCoR/ua4H/Zx8/AXxfnK+Wa4H1xph+4JCI1Ft7RLIpInuAVcAnbJ+Hrd0fJjuGMeZV4ASAMea0tT07bN4ZTXe/76yu77WiNtEdJX/A8HZD+6AQ7TreCcCsolxuuGg2qxZVcNn8MvKydX1IyUyWV5cAsO1QKx+6YGaKZxObRERpNnAs5HkDsCJaH2OMT0Q6gFLb/lrYtbPt40g2S4F2Y4wvQv+RjAGAiFQBFwFbI71BEbkduB1g7ty5kbpMSLrOUpR6rKdkJmD8rqPXy4v7mtiyt5Hn9zXR2j2AS+Dic6bxd2sWsmpRhYbllEnDBbOLyMtys/VgS0aIUqS/2vBPqWh9orVHKn0cq/9IxnAuEpkC/Br4gjGmM0JfjDEPAA8ALFu2bOJ9AkfhbD2lvrMM/40nxhj2nepi895GttQ18saRNvwBw7T8LK5aWMFVC8u58txyivN1/5Ay+chyu7j4nGlsPdSa6qnEJRFRagDmhDyvBI5H6dMgIh6gCGiNc22k9magWEQ81lsK7Z/0GCKShSNIjxhjfpPAe80oTvc5opTtGdnxBz0DzvXpqtK9A35ePdjsCNHeJt5r7wWgdmYhf3nlfK5eVMHSOcWaLacoOCcWf/f3+2jvGUjrL2eJiNLrQI3NinsPJ6ngE2F9NgC3Aq8CNwKbjTFGRDYAvxSR7+IkIdQA23C8mzNs2mu2WBvrrc0nRzKGXW96CNhjjPlusjcmE+js9QKQP8K1kg57fTpF7xraegbXhl450EK/L0B+tpuVC8q4Y9UCrl5YwYyi3FRPU1HSjuXVJRgDrx9u44NpfBptXFGy6zd3AJtw0rd/YozZJSJfB7YbYzbgfPj/wiYZtOKIDLbf4zjJBT7gc8YYP0Akm3bIO4H1IvJNYIe1TbJjiMjlwCeBnSLylrXxf40xG0d2qyYeQU8pPyt5UeoZ8NHnTf0+JZ8/wBtH2thc18iWvY3sO9UFwDml+dy8fC6rFlWwYl6JbmJVlDhcOKeYbI+LrQdbJrYoAdgP8o1hbV8NedwH3BTl2nuAexKxadsPMpShF9qe1BjGmJeIvN40aejsczyd3BF4Si1dA4OPx9tRaunq5/m6JjbXNfLiviZO9/nIcgvLq0v4k2VzuHpRBfPKCjRJQVGSIDfLzdI5xWw7nN7rSlrRIYNp7XaEZSRexEg33I6E4NlDm21Y7u2GdoyB8qk5XH/+DFYtqmDlgjKtLacoZ8kl1SV8f0s9p/u8afv3pKKUwRy3C//Byg7JcNhWNphbkj8mi0pd/T5e2t/MFpst13i6HxFYUlnMF1afy6pFFSyeVYhLkxQUZdRYXl1KYHM9bxxp46qFFameTkRUlDKY4+19AHhHcARFfVMXbpcwtyR/MAx4thxr7eH3e06xeW8jrx1swes3TM3xcMW55Vy9yEnbTvfd5ooykXnfOcV4XMLWQ60qSsr4EggY9pxwtmX5R+ApbT/cyvmzi8hyj9xT8fkD7DjWznN7Gnluzyn2NzpJCvPLC/j0ymquXljBsqppZLlHlrKuKEpy5Gd7WFJZxNaD6VsHT0UpQ9l1vJPT/T5cAr4kjzVv7OzjzaPt/MWV89h9vDOp6F2f18/zdY1s2nWKLXWNtPd48biEFfNKWLd8LqsXVVBVVpDku1EUZbRYXl3Kg384SO+APy1La6koZSg/feUQ2W4X19RWsP1wW1LXPvDiQfwBw40Xz+EbJ3bHPU+pd8ARot/uPMHmvY30DPiZlp/FqkUVrF40nQ+cWzYhDhdTlMnAinkl/OiFA7x5tI2VC8pSPZ0zUFHKQB7ddpTfvPkef3HlfLr7ffiSCN89u/sUP3n5EDcvn0N1WQFul0T1tHY2dLD+9aNseOs4p/t9lBZkc8NFs/nwBTNZUV2CR8NyipJ2LDtnGi6BrQdbVJSUsaXf5+efNu7lZ68c5spzy/niB2v41u/2JnRYnz9gePAPB/n2pjouqCzmyx+uBSDH42LAN/z6l+ub+bfn9rPtUCs5HhcfXjKTG99XyYp5pVrSR1HSnKm5WSyeVcRraVoHT0UpQ3i+rpF/+J/dHGru5jMrq7nr+kVke1x4Yng64OwRemFfE/+0cS91p05z/fkzuPfGJUzJcX41cjzuwcP+WrsH+Mp/v8tvd55gemEOX/lILTdeXElRnobmFGUisaK6hJ+/doQ+r5/cEVR8GUtUlCY42w+38u+b63lhXxPzygr4+WeWc8W55YOve9wufBFSwvu8fv7n7eP89OXD7D7RydySfO7/xPv40AUzhlVKyPa46PcFONnRxyf+4zUa2nr5m2vP5c8+MC/tfpkVRUmMFfNKefClQ7x9rJ0V80pTPZ1hqChNQLz+AM/taeSnLx9i66FWSgqyufv6RXx6ZfUZFcGzXILXbzDGICI0nu7jkdeO8sjWIzR3DXDu9Cl8648u4I/eVxmxmniOx0XvgI+/fOQNTnX28cs/X8GyqpLxequKoowBy6tKEIGth1pVlJSRU9/Yxa/eOMav32iguWuAmUW5fOUjtdy8fA752ZH/K4PJBq3dA/zLs/v41fZjeP2GVYsq+MzKalYuKI1ZQy4ny0X3gJ8dR9v5148vVUFSlAygKD+LhdOnsvVQC87BCumDilIaY4zhnYYONu06yTO7T1Hf6FRZWLWogpuXz+GKmvK4GW4eu/n1j3/4Csfaerl5+Rw+s7KaeeVTEppDsG5eTcUU1i6ddXZvSFGUtOGSeaWsf/0oA77AiM9cGwtUlNIMrz/A1oOtPLP7JM/sOsXJzj7cLmFFdQmfvOQcrj9/BhWFiZ8XlOVyftkOt/Scsd6UCO+bWwzA51fXaFVuRckgVlSX8LNXDrPzvQ4uPmdaqqcziIpSGtAz4OPFfU1s2nWK5/acorPPR26WiyvPLedvaxey+ryKEZ8UGdyxfWFlUdKCBHDVwgq2fXk1FVP14DxFySSWVzuh+K2HWlSUFGjrHuD3e06xadcp/rC/iX5fgOL8LD5YO4PrFk/nAzXlo1ICZOWCMiqn5fG1jy0esQ0VJEXJPEqn5LCgYgpbD7byV1elejZDqCiNI42dfWzadZKnd53ktYOt+AOGWUW53Lx8Ltcuns7yqtGvglBdVsBLd64aVZuKomQGK6pLePKt4/j8gbSpwKKiNMb0ef1s2nWSX21v4OUDzRgD88oL+Isr53Hd4hlcMLtI12oURUkJK+aV8sjWo+w+0cmSyuJUTwdQURoz2nsG+Nkrh/nZK4dp7/FSOS2Pz6+q4SNLZlIzfWqqp6coisIldl3p1QMtKkqZij9g+OXWI/zzpjpO9/m45rzpfHplFZfOK9VTVBVFSSsqCnOpnVnIs7tP8dkr56d6OgAkFEQUkTUiUici9SJyV4TXc0TkMfv6VhGpCnntbtteJyLXxbMpItXWxn5rM3u0xxgr2nsG+NRPt/GVJ3dxYWUxT3/hAzx46zJWLihTQVIUJS1Zc/4M3jjaRuPpvlRPBUhAlETEDdwPXA/UAjeLSG1Yt9uANmPMAuA+4F57bS2wDlgMrAF+ICLuODbvBe4zxtQAbdb2aI8x6rR1D3DTj15l68FW/umPLuAXty1n0YzCsRpOURRlVLhu8QyMgafePpHqqQCJeUrLgXpjzEFjzACwHlgb1mct8LB9/ASwWpzV+7XAemNMvzHmEFBv7UW0aa9ZZW1gbd4wmmMkdluSY8AX4FM/e50jrT387NPv5+blczV5QVGUCcG506ewvKqE723ez/H23lRPJ6E1pdnAsZDnDcCKaH2MMT4R6QBKbftrYdfOto8j2SwF2o0xvgj9R2uMMxCR24Hb7dMuEWkBmiP1jcfKe0ZyVdpSxgjvQwai92IIvRcOGXcfZn9txJeWAeeMxhwSEaVIX/nDD+iJ1idaeyQPLVb/0RzjzEZjHgAeCD4Xke3GmGWR+k4m9D4MofdiCL0XDnofhrD3omo0bCUSvmsA5oQ8rwSOR+sjIh6gCGiNcW209mag2NoIH2u0xlAURVHSlERE6XWgxmbFZeMkFWwI67MBuNU+vhHYbIwxtn2dzZyrxqmRvi2aTXvNFmsDa/PJ0RwjsduiKIqipIK44Tu7fnMHsAlwAz8xxuwSka8D240xG4CHgF+ISD2O97LOXrtLRB4HdgM+4HPGGD9AJJt2yDuB9SLyTWCHtc0ojxGPTCcqJAAAA8RJREFUB+J3mRTofRhC78UQei8c9D4MMWr3QhxnQ1EURVFST3pU4FMURVEUVJQURVGUNEJFKYTxLkuUCkTkJyLSKCLvhrSViMiztrTTsyIyzbaLiHzP3o93ROR9IdfcavvvF5FbI42VzojIHBHZIiJ7RGSXiPwf2z4Z70WuiGwTkbftvfgH2z5qJb8mErYizA4Reco+n6z34bCI7BSRt0Rku20b+78PY4z+c9bV3MABYB6QDbwN1KZ6XmPwPq8A3ge8G9L2z8Bd9vFdwL328YeA3+HsBbsE2GrbS4CD9uc0+3haqt9bkvdhJvA++3gqsA+nHNVkvBcCTLGPs4Ct9j0+Dqyz7T8C/tI+/ivgR/bxOuAx+7jW/t3kANX278md6vc3gvvxJeCXwFP2+WS9D4eBsrC2Mf/7UE9piHErS5RKjDEv4mQvhhJawim8tNPPjcNrOHvIZgLXAc8aY1qNMW3Aszh1BycMxpgTxpg37ePTwB6cSiCT8V4YY0yXfZpl/xlGr+TXhEFEKoEPAw/a56NZ+iwTGPO/DxWlISKVU5odpW+mMd0YcwKcD2ugwrZHuycZda9s2OUiHA9hUt4LG7J6C2jE+eA4QIIlv4DQkl8T/V78K/B3QMA+T7j0GZl1H8D5YvKMiLwhTik2GIe/Dz1PaYhEyilNNpIt7TThEJEpwK+BLxhjOiV6Id2MvhfG2du3VESKgf8CzovUzf7MyHshIh8BGo0xb4jIVcHmCF0z+j6EsNIYc1xEKoBnRWRvjL6jdi/UUxpiMpclOmVdbezPRtue0SWcRCQLR5AeMcb8xjZPynsRxBjTDjyPsy4wWiW/JgorgY+JyGGc8P0qHM9pst0HAIwxx+3PRpwvKssZh78PFaUhJnNZotASTuGlnf7UZtZcAnRYl30TcK2ITLPZN9fatgmDjf0/BOwxxnw35KXJeC/KrYeEiOQB1+CssY1Wya8JgTHmbmNMpXEKi67DeV//m0l2HwBEpEBEpgYf4/xev8t4/H2kOsMjnf7hZJDsw4mnfznV8xmj9/gocALw4nyLuQ0nDv4csN/+LLF9BeegxAPATmBZiJ3P4Czg1gOfTvX7GsF9uBwnjPAO8Jb996FJei+W4JT0esd+8HzVts/D+TCtB34F5Nj2XPu83r4+L8TWl+09qgOuT/V7O4t7chVD2XeT7j7Y9/y2/bcr+Hk4Hn8fWmZIURRFSRs0fKcoiqKkDSpKiqIoStqgoqQoiqKkDSpKiqIoStqgoqQoiqKkDSpKiqIoStqgoqQoiqKkDf8/Z+T2JUHicG0AAAAASUVORK5CYII=\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",
"constraint1 = tf.pow((sum_1-ztf.constant(1.7*10**-8))/ztf.constant(2.2*10**-8),2)/ztf.constant(2.)\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\n",
"\n",
"# constraint4 = triGauss(bplus_0, bplus_1, bplus_2)\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",
"constraint4_0 = tf.pow((bplus_0-mean1)/sigma1,2)/ztf.constant(2.)\n",
"constraint4_1 = tf.pow((bplus_1-mean2)/sigma2,2)/ztf.constant(2.)\n",
"constraint4_2 = tf.pow((bplus_2-mean3)/sigma3,2)/ztf.constant(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",
"\n",
"# constraint6_0 = zfit.constraint.GaussianConstraint(params = jpsi_p, mu = ztf.constant(pdg[\"jpsi_phase_unc\"]),\n",
"# sigma = ztf.constant(jpsi_phase))\n",
"# constraint6_1 = zfit.constraint.GaussianConstraint(params = psi2s_p, mu = ztf.constant(pdg[\"psi2s_phase_unc\"]),\n",
"# sigma = ztf.constant(psi2s_phase))\n",
"\n",
"constraint6_0 = tf.pow((jpsi_p-ztf.constant(jpsi_phase))/ztf.constant(pdg[\"jpsi_phase_unc\"]),2)/ztf.constant(2.)\n",
"constraint6_1 = tf.pow((psi2s_p-ztf.constant(psi2s_phase))/ztf.constant(pdg[\"psi2s_phase_unc\"]),2)/ztf.constant(2.)\n",
"\n",
"# zfit.run(constraint6_0)\n",
"\n",
"# ztf.convert_to_tensor(constraint6_0)\n",
"\n",
"#List of all constraints\n",
"\n",
"constraints = [constraint1, constraint2, constraint3_0, constraint3_1, constraint4_0, constraint4_1, constraint4_2,\n",
" constraint6_0, constraint6_1]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Reset params"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [],
"source": [
"def reset_param_values(): \n",
" jpsi_m.set_value(jpsi_mass)\n",
" jpsi_s.set_value(jpsi_scale)\n",
" jpsi_p.set_value(jpsi_phase)\n",
" jpsi_w.set_value(jpsi_width)\n",
" psi2s_m.set_value(psi2s_mass)\n",
" psi2s_s.set_value(psi2s_scale)\n",
" psi2s_p.set_value(psi2s_phase)\n",
" psi2s_w.set_value(psi2s_width)\n",
" p3770_m.set_value(p3770_mass)\n",
" p3770_s.set_value(p3770_scale)\n",
" p3770_p.set_value(p3770_phase)\n",
" p3770_w.set_value(p3770_width)\n",
" p4040_m.set_value(p4040_mass)\n",
" p4040_s.set_value(p4040_scale)\n",
" p4040_p.set_value(p4040_phase)\n",
" p4040_w.set_value(p4040_width)\n",
" p4160_m.set_value(p4160_mass)\n",
" p4160_s.set_value(p4160_scale)\n",
" p4160_p.set_value(p4160_phase)\n",
" p4160_w.set_value(p4160_width)\n",
" p4415_m.set_value(p4415_mass)\n",
" p4415_s.set_value(p4415_scale)\n",
" p4415_p.set_value(p4415_phase)\n",
" p4415_w.set_value(p4415_width)\n",
" rho_m.set_value(rho_mass)\n",
" rho_s.set_value(rho_scale)\n",
" rho_p.set_value(rho_phase)\n",
" rho_w.set_value(rho_width)\n",
" omega_m.set_value(omega_mass)\n",
" omega_s.set_value(omega_scale)\n",
" omega_p.set_value(omega_phase)\n",
" omega_w.set_value(omega_width)\n",
" phi_m.set_value(phi_mass)\n",
" phi_s.set_value(phi_scale)\n",
" phi_p.set_value(phi_phase)\n",
" phi_w.set_value(phi_width)\n",
" Dstar_m.set_value(Dstar_mass)\n",
" DDstar_s.set_value(0.0)\n",
" DDstar_p.set_value(0.0)\n",
" D_m.set_value(D_mass)\n",
" Dbar_m.set_value(Dbar_mass)\n",
" Dbar_s.set_value(0.0)\n",
" Dbar_p.set_value(0.0)\n",
" tau_m.set_value(pdg['tau_M'])\n",
" Ctt.set_value(0.0)\n",
" b0_0.set_value(0.292)\n",
" b0_1.set_value(0.281)\n",
" b0_2.set_value(0.150)\n",
" bplus_0.set_value(0.466)\n",
" bplus_1.set_value(-0.885)\n",
" bplus_2.set_value(-0.213)\n",
" bT_0.set_value(0.460)\n",
" bT_1.set_value(-1.089)\n",
" bT_2.set_value(-1.114)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Analysis"
]
},
{
"cell_type": "code",
"execution_count": 38,
"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",
"WARNING:tensorflow:From C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow_probability\\python\\distributions\\categorical.py:263: multinomial (from tensorflow.python.ops.random_ops) is deprecated and will be removed in a future version.\n",
"Instructions for updating:\n",
"Use tf.random.categorical instead.\n",
"Toy 0: Generating data...\n",
"Toy 0: Data generation finished\n",
"Toy 0: Loading data...\n",
"Toy 0: Loading data finished\n",
"Toy 0: Fitting pdf...\n",
"------------------------------------------------------------------\n",
"| FCN = 291.3 | Ncalls=753 (753 total) |\n",
"| EDM = 4.03E-05 (Goal: 5E-06) | up = 0.5 |\n",
"------------------------------------------------------------------\n",
"| Valid Min. | Valid Param. | Above EDM | Reached call limit |\n",
"------------------------------------------------------------------\n",
"| True | True | False | False |\n",
"------------------------------------------------------------------\n",
"| Hesse failed | Has cov. | Accurate | Pos. def. | Forced |\n",
"------------------------------------------------------------------\n",
"| False | True | True | True | False |\n",
"------------------------------------------------------------------\n",
"Function minimum: 291.34660987562074\n",
"----------------------------------------------------------------------------------------------\n",
"| | Name | Value | Hesse Err | Minos Err- | Minos Err+ | Limit- | Limit+ | Fixed |\n",
"----------------------------------------------------------------------------------------------\n",
"| 0 | bplus_2 | 2.0 | 3.3 | | | -2 | 2 | |\n",
"| 1 | p4415_p | -4.1 | 1.8 | | |-6.28319 | 6.28319 | |\n",
"| 2 | p4040_s | 2.0 | 2.0 | | |0.00501244| 2.01499 | |\n",
"| 3 | Ctt | -0.5 | 0.5 | | | -0.5 | 0.5 | |\n",
"| 4 | bplus_0 | 0.27 | 0.15 | | | -2 | 2 | |\n",
"| 5 | p4160_s | 3.3 | 2.5 | | | 0.71676 | 3.68324 | |\n",
"| 6 | DDstar_s | 0.30 | 0.55 | | | -0.3 | 0.3 | |\n",
"| 7 | p4415_s | 2.4 | 1.4 | | |0.126447 | 2.35355 | |\n",
"| 8 | Dbar_p | 6 | 10 | | |-6.28319 | 6.28319 | |\n",
"| 9 | p4160_p | 2.4 | 1.9 | | |-6.28319 | 6.28319 | |\n",
"| 10| p3770_s | 4.1 | 3.0 | | |0.918861 | 4.08114 | |\n",
"| 11| p3770_p | -3.8 | 1.7 | | |-6.28319 | 6.28319 | |\n",
"| 12| jpsi_p | 5.0 | 0.8 | | |-6.28319 | 6.28319 | |\n",
"| 13| psi2s_p | 1.80 | 0.03 | | |-6.28319 | 6.28319 | |\n",
"| 14| Dbar_s | 0.30 | 0.59 | | | -0.3 | 0.3 | |\n",
"| 15| bplus_1 | -0.48 | 0.28 | | | -2 | 2 | |\n",
"| 16| p4040_p | 5.4 | 1.9 | | |-6.28319 | 6.28319 | |\n",
"| 17| DDstar_p | -6 | 9 | | |-6.28319 | 6.28319 | |\n",
"----------------------------------------------------------------------------------------------\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"| | bplus_2 p4415_p p4040_s Ctt bplus_0 p4160_s DDstar_s p4415_s Dbar_p p4160_p p3770_s p3770_p jpsi_p psi2s_p Dbar_s bplus_1 p4040_p DDstar_p |\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"| bplus_2 | 1.000 0.002 0.000 -0.000 -0.024 -0.001 -0.000 -0.000 0.000 0.001 0.000 -0.001 0.002 0.000 -0.000 0.024 -0.001 -0.000 |\n",
"| p4415_p | 0.002 1.000 -0.002 -0.008 0.157 -0.212 0.001 0.005 0.092 0.237 -0.000 0.020 -0.063 0.001 0.003 0.164 -0.108 0.008 |\n",
"| p4040_s | 0.000 -0.002 1.000 0.000 0.002 -0.000 -0.000 0.000 0.001 -0.004 -0.000 0.001 0.001 0.000 -0.000 -0.003 0.002 -0.000 |\n",
"| Ctt | -0.000 -0.008 0.000 1.000 -0.060 0.010 0.000 -0.000 -0.010 -0.010 -0.000 -0.015 0.012 0.000 0.000 0.030 -0.021 0.000 |\n",
"| bplus_0 | -0.024 0.157 0.002 -0.060 1.000 0.134 0.001 -0.003 -0.129 0.005 0.001 -0.107 0.068 0.003 0.003 -0.536 -0.098 0.000 |\n",
"| p4160_s | -0.001 -0.212 -0.000 0.010 0.134 1.000 -0.005 -0.001 0.142 -0.561 -0.001 0.164 0.027 0.005 -0.002 -0.033 0.010 -0.001 |\n",
"| DDstar_s | -0.000 0.001 -0.000 0.000 0.001 -0.005 1.000 0.000 0.001 0.005 0.000 -0.002 0.001 0.000 -0.000 -0.015 -0.001 0.001 |\n",
"| p4415_s | -0.000 0.005 0.000 -0.000 -0.003 -0.001 0.000 1.000 -0.001 -0.005 -0.000 -0.001 -0.000 -0.000 0.000 0.003 0.001 0.000 |\n",
"| Dbar_p | 0.000 0.092 0.001 -0.010 -0.129 0.142 0.001 -0.001 1.000 0.016 -0.001 0.274 0.209 0.000 0.003 0.371 0.160 0.003 |\n",
"| p4160_p | 0.001 0.237 -0.004 -0.010 0.005 -0.561 0.005 -0.005 0.016 1.000 0.000 -0.059 -0.079 -0.002 0.004 0.216 -0.119 0.007 |\n",
"| p3770_s | 0.000 -0.000 -0.000 -0.000 0.001 -0.001 0.000 -0.000 -0.001 0.000 1.000 0.003 0.000 -0.000 0.000 -0.002 0.001 0.000 |\n",
"| p3770_p | -0.001 0.020 0.001 -0.015 -0.107 0.164 -0.002 -0.001 0.274 -0.059 0.003 1.000 -0.027 -0.003 0.005 0.243 0.114 0.007 |\n",
"| jpsi_p | 0.002 -0.063 0.001 0.012 0.068 0.027 0.001 -0.000 0.209 -0.079 0.000 -0.027 1.000 -0.003 0.003 -0.264 0.078 0.007 |\n",
"| psi2s_p | 0.000 0.001 0.000 0.000 0.003 0.005 0.000 -0.000 0.000 -0.002 -0.000 -0.003 -0.003 1.000 0.000 0.005 0.001 0.000 |\n",
"| Dbar_s | -0.000 0.003 -0.000 0.000 0.003 -0.002 -0.000 0.000 0.003 0.004 0.000 0.005 0.003 0.000 1.000 -0.006 0.003 0.000 |\n",
"| bplus_1 | 0.024 0.164 -0.003 0.030 -0.536 -0.033 -0.015 0.003 0.371 0.216 -0.002 0.243 -0.264 0.005 -0.006 1.000 -0.094 0.001 |\n",
"| p4040_p | -0.001 -0.108 0.002 -0.021 -0.098 0.010 -0.001 0.001 0.160 -0.119 0.001 0.114 0.078 0.001 0.003 -0.094 1.000 0.006 |\n",
"| DDstar_p | -0.000 0.008 -0.000 0.000 0.000 -0.001 0.001 0.000 0.003 0.007 0.000 0.007 0.007 0.000 0.000 0.001 0.006 1.000 |\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"Hesse errors: OrderedDict([(<zfit.Parameter 'bplus_2' floating=True>, {'error': 3.325119189670308}), (<zfit.Parameter 'p4415_p' floating=True>, {'error': 1.7815174671180474}), (<zfit.Parameter 'p4040_s' floating=True>, {'error': 1.9926966547547778}), (<zfit.Parameter 'Ctt' floating=True>, {'error': 0.5114351017394924}), (<zfit.Parameter 'bplus_0' floating=True>, {'error': 0.14819864700462704}), (<zfit.Parameter 'p4160_s' floating=True>, {'error': 2.4730090801780316}), (<zfit.Parameter 'DDstar_s' floating=True>, {'error': 0.5456802813189129}), (<zfit.Parameter 'p4415_s' floating=True>, {'error': 1.4247502221049095}), (<zfit.Parameter 'Dbar_p' floating=True>, {'error': 9.742999517381843}), (<zfit.Parameter 'p4160_p' floating=True>, {'error': 1.9168225450258651}), (<zfit.Parameter 'p3770_s' floating=True>, {'error': 2.9995956487093722}), (<zfit.Parameter 'p3770_p' floating=True>, {'error': 1.7105494289800758}), (<zfit.Parameter 'jpsi_p' floating=True>, {'error': 0.7843308347897375}), (<zfit.Parameter 'psi2s_p' floating=True>, {'error': 0.03332981923505862}), (<zfit.Parameter 'Dbar_s' floating=True>, {'error': 0.5881037384446894}), (<zfit.Parameter 'bplus_1' floating=True>, {'error': 0.27508756809324075}), (<zfit.Parameter 'p4040_p' floating=True>, {'error': 1.8523591244884532}), (<zfit.Parameter 'DDstar_p' floating=True>, {'error': 9.488064745244461})])\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\ipykernel_launcher.py:160: UserWarning: Creating legend with loc=\"best\" can be slow with large amounts of data.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Toy 1/10\n",
"Time taken: 4 min, 1 \n",
"Projected time left: 36 min, 10 s\n",
"Toy 1: Generating data...\n"
]
},
{
"ename": "KeyboardInterrupt",
"evalue": "",
"output_type": "error",
"traceback": [
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[1;31mKeyboardInterrupt\u001b[0m Traceback (most recent call last)",
"\u001b[1;32m<ipython-input-38-a8ff80ea8d66>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m 37\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 38\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 39\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 40\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 41\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 1332\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m_do_call\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfn\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 1333\u001b[0m \u001b[1;32mtry\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1334\u001b[1;33m \u001b[1;32mreturn\u001b[0m \u001b[0mfn\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 1335\u001b[0m \u001b[1;32mexcept\u001b[0m \u001b[0merrors\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mOpError\u001b[0m \u001b[1;32mas\u001b[0m \u001b[0me\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 1336\u001b[0m \u001b[0mmessage\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mcompat\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mas_text\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0me\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmessage\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32m~\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\client\\session.py\u001b[0m in \u001b[0;36m_run_fn\u001b[1;34m(feed_dict, fetch_list, target_list, options, run_metadata)\u001b[0m\n\u001b[0;32m 1317\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_extend_graph\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 1318\u001b[0m return self._call_tf_sessionrun(\n\u001b[1;32m-> 1319\u001b[1;33m options, feed_dict, fetch_list, target_list, run_metadata)\n\u001b[0m\u001b[0;32m 1320\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 1321\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m_prun_fn\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mhandle\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfeed_dict\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfetch_list\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32m~\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\client\\session.py\u001b[0m in \u001b[0;36m_call_tf_sessionrun\u001b[1;34m(self, options, feed_dict, fetch_list, target_list, run_metadata)\u001b[0m\n\u001b[0;32m 1405\u001b[0m return tf_session.TF_SessionRun_wrapper(\n\u001b[0;32m 1406\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_session\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0moptions\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfeed_dict\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfetch_list\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtarget_list\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1407\u001b[1;33m run_metadata)\n\u001b[0m\u001b[0;32m 1408\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 1409\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m_call_tf_sessionprun\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mhandle\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfeed_dict\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfetch_list\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;31mKeyboardInterrupt\u001b[0m: "
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ8AAAD8CAYAAACo9anUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3de5hc1Xnn++9bVX3R/dJIMkgCCUtYtBzARmAnEBuDbWSbWCSBIGawmWPOYUjsZGLPeIDJhIwZmBN8PGYmx1dseIbBGQsMONY42PgCxMEXYcmIiwSyWhKgRoAuLbWkbvWlqt75Y6+qripV1S7dunfRv8/z9FO7Vu291qot9X57rb3W2ubuiIiIjKbUWFdARETGHwUfEREZdQo+IiIy6hR8RERk1Cn4iIjIqFPwERGRUddQ8DGz5Wa2ycy6zOymKp+3mdn94fM1Zrag5LObQ/omM7s0Lk8zWxjy2BzybG2gjLPM7JdmtsHMnjOz9qM5GSIiMjpig4+ZpYEvAx8COoGrzayzYrfrgL3uvgi4E7gjHNsJrASWAsuBr5hZOibPO4A73X0xsDfkXa+MDPAt4AZ3XwpcBAwf4XkQEZFR1EjL53ygy923uvsQsApYUbHPCuDesP0gcImZWUhf5e6D7r4N6Ar5Vc0zHHNxyIOQ5+UxZXwQeNbdnwFw9z3unmv8FIiIyGjLNLDPXGB7yftu4F219nH3rJn1Ah0h/VcVx84N29Xy7AD2uXu2yv61yjgDcDN7FJhFFOw+X/klzOx64HqASZMmnbtkyZIGvrqIAAznnBdf38+syW3sOjjIvBkTmDGxdayrVXRwMMu23X0AvG3OFFozup19Iqxbt263u886Hnk1EnysSlrlmjy19qmVXu1/Rr3965WRAS4EzgP6gZ+a2Tp3/2nZju53AXcBLFu2zNeuXVslOxGp5vXeAd79//6UP73orXz1iS3cceXZXHHuvLGuVtEvtuzmX3xjDQDf/+z7OLVj4hjX6M3JzF4+Xnk18udBNzC/5P08YEetfcI9mGlAT51ja6XvBqaHPCrLqlfGP7n7bnfvBx4B3tnA9xKRBnn4GzAV/gRM8pqQftjfxpJEjQSfXwOLwyi0VqIBBKsr9lkNXBu2rwAe8+h/52pgZRipthBYDDxVK89wzOMhD0Ke34sp41HgLDObGILSe4GNjZ8CEYlTiDUpi6JP4i7viauQxIntdgv3Vz5FdJFPA/e4+wYzuxVY6+6rgbuB+8ysi6g1sjIcu8HMHiAKBlngk4XBANXyDEXeCKwys9uAp0Pe1Cljr5l9kSigOfCIu//jMZ0VESlT7Ps2K09IoAQ3yqREI/d8cPdHiLqzStNuKdkeAK6sceztwO2N5BnStxKNhqtMr1fGt4iGW4vICVTtxqtEhoeH6e7uZmBgYKyrcsza29uZN28eLS0tJ6yMhoKPiIxvhXs8Iw2f5DYvxqpm3d3dTJkyhQULFoy0EJuQu7Nnzx66u7tZuHDhCStH4xFFJNZh93wSFnuSUJ2BgQE6OjqaOvBA1LXa0dFxwltwCj4i0rBU8m/5jOlIvGYPPAWj8T0UfEQkVuF6bglt+UjzUfARkYa9Sf6wH9eeeOIJLrvsMgAGBwd5//vfzznnnMP9998/qvXQgAMRiVUYYGBY2fukKG2JJatmyfb0008zPDzM+vXrR71stXxEJNbIgIPy95IsL730EkuWLOHaa6/lrLPO4oorrqC/v58f/vCHLFmyhAsvvJCHH34YgJ07d3LNNdewfv16zjnnHLZs2TKqdVXLR0RiFWJNYlc4KJGEwPi5/72BjTv2H9c8O0+Zyt/8wdLY/TZt2sTdd9/NBRdcwCc+8Qm++MUv8vWvf53HHnuMRYsWcdVVVwEwe/ZsvvnNb/KFL3yB73//+8e1ro1Qy0dEYlXO80nEFb5E0roBx9L8+fO54IILALjmmmtYu3YtCxcuZPHixZgZ11xzzRjXMKKWj4jEOmx5nUQb+0DUSAvlRKn8N+rt7U3kv5taPiLSsGLDZ0xrcbiENcTG1CuvvMIvf/lLAL797W/z/ve/n23bthXv6Xz7298ey+oVKfiISKxmGnCQ5LqNhjPPPJN7772Xs846i56eHj796U9z11138ZGPfIQLL7yQ0047bayrCKjbTUQaEp7nkypMMk3WFT5ZtRlbqVSKr33ta2Vpy5cv58UXXzxs34suuoiLLrpolGpWTi0fEYl12AoHY1iXOEmum4xQ8BGRWPWeZS/JsWDBAp5//vmxrkZDFHxEpGGW0Hs+SekGTEo9jtVofA8FHxGJlfjHaJcYq+t/e3s7e/bsafoAVHieT3t7+wktRwMORCRWYRLnyGi3ZF1gk1CbefPm0d3dza5du8a6Kses8CTTE0nBR0RiVQ44SLKxWu2gpaXlhD75881G3W4iEqsYfMa2GrUloekjR0TBR0RijXS7Jf9hckmum4xQ8BGRhhVHuyWsqZG0+kg8BR8RiXXYaLcEX+uTXDcZoeAjIg0bafkkiwJO81HwEZFYGu0mx5uCj4jEOnyezxhWpoqk1UfiKfiISMOMwgoHutrLsWko+JjZcjPbZGZdZnZTlc/bzOz+8PkaM1tQ8tnNIX2TmV0al6eZLQx5bA55ttYrw8wWmNkhM1sffsrXEheRY5b05/mUVidpdZPqYoOPmaWBLwMfAjqBq82ss2K364C97r4IuBO4IxzbCawElgLLga+YWTomzzuAO919MbA35F2zjGCLu58Tfm44ojMgIrGa6zHa0gwaafmcD3S5+1Z3HwJWASsq9lkB3Bu2HwQuseh/6QpglbsPuvs2oCvkVzXPcMzFIQ9CnpfHlCEiJ1hhLbek/sYlba05iddI8JkLbC953x3Squ7j7lmgF+ioc2yt9A5gX8ijsqxaZQAsNLOnzeyfzOz3G/hOInIECpf2kXk+utjLsWlkYdFqf+tU/s+rtU+t9GpBr97+9cp4DTjV3feY2bnAP5jZUnffX1ZBs+uB6wFOPfXUKlmJSC265yPHWyMtn25gfsn7ecCOWvuYWQaYBvTUObZW+m5gesijsqyqZYQuvT0A7r4O2AKcUfkl3P0ud1/m7stmzZrVwNcWkUpJnWQqzaeR4PNrYHEYhdZKNIBgdcU+q4Frw/YVwGMetctXAyvDSLWFwGLgqVp5hmMeD3kQ8vxevTLMbFYYwICZnR7K2Nr4KRCReIV7Pk2wvI5CY1OI7XZz96yZfQp4FEgD97j7BjO7FVjr7quBu4H7zKyLqMWzMhy7wcweADYCWeCT7p4DqJZnKPJGYJWZ3QY8HfKmVhnAe4BbzSwL5IAb3L3n6E+JiFRK+iMVkhwMpbqGHibn7o8Aj1Sk3VKyPQBcWePY24HbG8kzpG8lGg1XmV61DHd/CHgo9kuIyFE7bMBBglsXCkTNQSsciEis5K9qnbgKSQwFHxFpmAYcyPGi4CMisQ6b15Owpk9pdZJVM6lFwUdEYpVOtkvqKgfSXBR8RCSWV0z1TlrronySadJqJ9Uo+IhIrMLoNgsPVdD1XY6Vgo+IxCs+yTSaaJq0oda659N8FHxEpGEWftTykWOl4CMisUpjTRIHHJS2xBQYm4OCj4jEKi6vU1zhQOTYKPiISKzigAOLBh0krXVRXp+EVU6qUvARkVhlC4tastd2k+ag4CMiDYtaPiSucZGw6kgDFHxEJFbSBxyUSlqXoFSn4CMisbxiiYOkXd+1qkHzUfARkVjF0FMccJDci31yayalFHxEJF7JgAMzdW3JsVPwEZFYI0Otw9puY1uduhQYm4OCj4g0LGr5JHzEgTQFBR8RiZXwZ8mVLyyatMpJVQo+IhJrZHmdsLBoojvepBko+IhIrJGB1lH0SVrjomxh0TGshzROwUdEYhW6soorHIgcIwUfETkiSRxwkLSWmMRT8BGRWJXX9iTf1E9w1aSEgo+IxCobcGDJu6+igNN8FHxEpAHhng+W+MdoayRec1DwEZFY5S0fS9wFPlm1kUY0FHzMbLmZbTKzLjO7qcrnbWZ2f/h8jZktKPns5pC+ycwujcvTzBaGPDaHPFvjygifn2pmB83s3x3pSRCR+soXFk04RaKmEBt8zCwNfBn4ENAJXG1mnRW7XQfsdfdFwJ3AHeHYTmAlsBRYDnzFzNIxed4B3Onui4G9Ie+aZZS4E/hBo19cRBo38iRTK3ufFEkeACHVNdLyOR/ocvet7j4ErAJWVOyzArg3bD8IXGLReMwVwCp3H3T3bUBXyK9qnuGYi0MehDwvjykDM7sc2ApsaPyri0ijCt1sqYQOOCiV5LrJiEaCz1xge8n77pBWdR93zwK9QEedY2uldwD7Qh6VZVUtw8wmATcCn6v3JczsejNba2Zrd+3aFfOVRaRUvuSeD1jyWj5jXQE5Yo0En2pdvJX/1rX2OV7p9cr4HFE33cEqn4/s6H6Xuy9z92WzZs2qt6uIVBhZ4cBCANLlXo5NpoF9uoH5Je/nATtq7NNtZhlgGtATc2y19N3AdDPLhNZN6f61yngXcIWZfR6YDuTNbMDdv9TAdxORBhRaOilL6FDrslWtx64a0rhGWj6/BhaHUWitRAMIVlfssxq4NmxfATzm0Z9Kq4GVYaTaQmAx8FStPMMxj4c8CHl+r14Z7v777r7A3RcA/w34Lwo8IsdXvtDyodD1JnJsYls+7p41s08BjwJp4B5332BmtwJr3X01cDdwn5l1EbVGVoZjN5jZA8BGIAt80t1zANXyDEXeCKwys9uAp0Pe1CpDRE680pZP6fskStocJKmukW433P0R4JGKtFtKtgeAK2sceztweyN5hvStRKPhKtNrllGyz3+q97mIHJ182arWSZxkmqz6SDytcCAiscommSbweT6lklw3GaHgIyKxCqPdigMOxrY6h1HAaT4KPiISq3SeTxKf51NKcag5KPiISKykDzhIWHWkAQo+IhKrdKg1JPsGv9Z5aw4KPiISa2TAQVjhIGHXd8Wb5qPgIyKxvHSodcIXFpXmoOAjIrHKl9dJ3oCD0m5ABcbmoOAjIrHyxaHW0XvdV5FjpeAjIrGKQ62xRHa7lcXCpFVOqlLwEZFYxXs+KZK5qrU0HQUfEYk18hjt6L5P0mJPecMnabWTahR8RCTWyGO0o263vJo+cowUfEQkVr5ktFvKLHkDDkrqk7SqSXUKPiISK18xzyefH+MKSdNT8BGRWMV7Pla455Pc5oVaPs1BwUdEYhVHu2GYWbEbLikSVh1pgIKPiMQaWeEg+kncPR9pOgo+IhJr5Hk+0YCDxLV8SuqTsKpJDQo+IhJrZKh19KOh1nKsFHxEJFZpy4dEtny86rYkl4KPiMRydwpPz9Y9HzkeFHxEJJb7yCO0o0mmY1yhClpXtPko+IhIrLx78Sk+uucjx4OCj4jEckZaPtE8n2QFn7LRbsmqmtSg4CMisfIV93ySNuBAmo+Cj4jEcqcYfIzkLSzqdd5JMjUUfMxsuZltMrMuM7upyudtZnZ/+HyNmS0o+ezmkL7JzC6Ny9PMFoY8Noc8W+uVYWbnm9n68POMmf3h0Z4MEanO3UcGHKTU8pFjFxt8zCwNfBn4ENAJXG1mnRW7XQfsdfdFwJ3AHeHYTmAlsBRYDnzFzNIxed4B3Onui4G9Ie+aZQDPA8vc/ZxQxtfNLHNkp0FE6sk7JQMOEtjySVh9JF4jLZ/zgS533+ruQ8AqYEXFPiuAe8P2g8AlZmYhfZW7D7r7NqAr5Fc1z3DMxSEPQp6X1yvD3fvdPRvS21GbW+S4Kx1qncSFRUspDjWHRoLPXGB7yfvukFZ1nxAIeoGOOsfWSu8A9pUEk9KyapWBmb3LzDYAzwE3lBxfZGbXm9laM1u7a9euBr62iBTkfaTpk8RJpgmrjjSgkeBjVdIq/6lr7XO80uvWw93XuPtS4DzgZjNrP2xH97vcfZm7L5s1a1aVrESklrJ7Pgls+ZQ+XyhhVZMaGgk+3cD8kvfzgB219gn3W6YBPXWOrZW+G5hecs+mtKxaZRS5+wtAH/D2Br6XiDQomucTbSdxkmnCqiMNaCT4/BpYHEahtRINIFhdsc9q4NqwfQXwmEft8tXAyjBSbSGwGHiqVp7hmMdDHoQ8v1evjJBHBsDMTgPeBrzU8BkQkVjRPJ9C50PyWj6l9VEgag6xo8LcPWtmnwIeBdLAPe6+wcxuBda6+2rgbuA+M+siao2sDMduMLMHgI1AFviku+cAquUZirwRWGVmtwFPh7ypVQZwIXCTmQ0DeeDP3H330Z8SEakUDTiIthN5z0edbU2noSHJ7v4I8EhF2i0l2wPAlTWOvR24vZE8Q/pWotFwlelVy3D3+4D7Yr+EiBy1qGWR4IVFyx4ml7DKSVVa4UBEGuAjLZ9U8u75SPNR8BGRWPl80hcWLX2Y3BhWRBqm4CMiscoXFk1et1vSBkBIPAUfEYmVL13hgOR1uyWsOtIABR8RiZV3JxWuFilL3kROTTJtPgo+IhIrl3fSZSscJOsSn7DqSAMUfEQkVtTyKRlwkB/jClUojT1Jm4Mk1Sn4iEisvJe2fJJ3gU9afSSego+IxMrlE76waEl9ktYlKNUp+IhIrFyekm635F3gSwccJK1LUKpT8BGRWHl30uFqYWaJG1FWvryONAMFHxGJVT7aLXn3WEprk7RWmVSn4CMisUpHuyX9no+aPs1BwUdEYlW2fJLWuihtiSWtblKdgo+IxDp8nk+yLvBq+DQfBR8RiZXPU2z5mCVvRQG1fJqPgo+IxMq5ky6555O0y3vZaLekVU6qUvARkVi5fOmAg+S1LrS8TvNR8BGRWNHyOtF2EhcWLa1PsmomtSj4iEis0uV1LOFDrZM2GEKqU/ARkViV3W5J7tpKbs2klIKPiMQqX9U6iS2fkcd8J61uUp2Cj4jEyjvF0W7JXFh0ZCh4kltlMkLBR0Ri5fPlk0yTdn0vnQSbtLpJdQo+IhIrVzLardDCyCWof8s9uhcF5Y9XkORS8BGRWKUDDgqPVkhU8GEkKCaoWlKHgo+IxMqXLCyaTkWXjSTd94laPup2ayYNBR8zW25mm8ysy8xuqvJ5m5ndHz5fY2YLSj67OaRvMrNL4/I0s4Uhj80hz9Z6ZZjZB8xsnZk9F14vPtqTISLVlS6vkwmv2UQ1MUZaZkkKilJbbPAxszTwZeBDQCdwtZl1Vux2HbDX3RcBdwJ3hGM7gZXAUmA58BUzS8fkeQdwp7svBvaGvGuWAewG/sDdfwe4FrjvyE6BiMQpfYx2IQjlcsm5yHvJaDxpDo20fM4Hutx9q7sPAauAFRX7rADuDdsPApeYmYX0Ve4+6O7bgK6QX9U8wzEXhzwIeV5erwx3f9rdd4T0DUC7mbU1egJEJF7pPJ9MutDyyY9llcrkfWQFBq1w0BwaCT5zge0l77tDWtV93D0L9AIddY6tld4B7At5VJZVq4xSfww87e6DlV/CzK43s7VmtnbXrl0xX1lESkUX92g7lfjRbtIMGgk+1dqylf++tfY5Xumx9TCzpURdcf+6yn64+13uvszdl82aNavaLiJSQ+lot8I9n1yC7q04I91uuufTHBoJPt3A/JL384AdtfYxswwwDeipc2yt9N3A9JBHZVm1ysDM5gHfBT7u7lsa+E4icgRyeS8GncJFPpuwez4a7dZcGgk+vwYWh1ForUQDCFZX7LOa6GY/wBXAYx6tcbEaWBlGqi0EFgNP1cozHPN4yIOQ5/fqlWFm04F/BG52958fyZcXkcZkc04mTPAp3PNJVLcb0dpulvBFT2VEbPAJ91c+BTwKvAA84O4bzOxWM/to2O1uoMPMuoDPADeFYzcADwAbgR8Cn3T3XK08Q143Ap8JeXWEvGuWEfJZBPy1ma0PP7OP8nyISBXD+TwtJU8yhWQNtXaPAo+hSabNIhO/C7j7I8AjFWm3lGwPAFfWOPZ24PZG8gzpW4lGw1WmVy3D3W8Dbov9EiJyVHJ5x52Rlk+YZJqolo87hoVHfCenXlKbVjgQkbqGc9GQ6kJ3W3GeT5KCDxS73RJULalDwUdE6ioEn5ZUoeWTwODjUZdbElfcluoUfESkrsKotmLLJ4GTTHPhkQqGBhw0CwUfEalrOF/odosuF0l8pEJh4dPono80AwUfEamr0PJpSfDCovmw8KmZltdpFgo+IlLXSLdbaPkUVhJI0EU+l4/u96jl0zwUfESkrkK3W0u6cmHR5Fzmo5ZPYZ5PcuoltSn4iEhdxZZPqtDySd48n8Kq1tEKB2NdG2mEgo+I1HXYPJ8ErnCQyxeCj2m0W5NQ8BGRugpBpiXJk0zDw+RSpkcqNAsFH2FgOMdz3b1jXQ1JqGyh5ZNK7sKiUcsnGnSgez7NQcFH+KvvPs8ffOlJXu8dGOuqSAINV04yTSV0kqmFlo9iT1NQ8BHWb98LwMHB4TGuiSRRtjjaLbmTTD3M8wHT2m5NQsFHROoaGe1W2fJJzlW+MOAgqlpy6iW1KfiISF3FhUUT/DC5vBOt7WaQoN5AqUPBR/R3otQ1kI2u5u0t0eWiEIQKQSkJonk+6Hk+TUTBR0TqGhzOAdCWSQPQmokuG0PZ5ASfXFhYVE8ybR4KPoKNdQUk0QZDkGkLLZ/W0PIZTFDwKXS7pVKWqDXnpDYFHxGpa6Cy5ZNOXssnH+b5ZFJGTmOtm4KCj6iHXOoqtnxCd1sqZbSkLVEtn1wYap1OWaJG4UltCj4iUtfIPZ+Ry0VbJp2slk+YZJpJpcjlFHyagYKP6J6P1DWYzdOWSWE28j+lNZNiKJcbw1qVy4d5Pmr5NA8FH1G3m9Q1mM3T3pIuS2tNpxLW8okmv2bSRk4TfZqCgo+I1DUwnCvrcoPQ8klQ8MnlHTPU8mkiCj5SpN9ZqWYwmy8Osy6Iut2SE3zyHs3zyaQsUSsvSG0KPlKkpeilmoHhHO2Zw7vdBoeTE3yGc04mnVLLp4ko+EiRusqlmoHh3GEtn7aWZLV8svk8Lekw2k3Bpyk0FHzMbLmZbTKzLjO7qcrnbWZ2f/h8jZktKPns5pC+ycwujcvTzBaGPDaHPFvrlWFmHWb2uJkdNLMvHe2JELQmllTVN5RjUmumLK01nUrUPJ/hbJ4WtXyaSmzwMbM08GXgQ0AncLWZdVbsdh2w190XAXcCd4RjO4GVwFJgOfAVM0vH5HkHcKe7Lwb2hrxrlgEMAH8N/Lsj/O5SQb1uUs3BgSyT2yqCT8IGHAznnUy6cM8nOfWS2hpp+ZwPdLn7VncfAlYBKyr2WQHcG7YfBC6xaFLACmCVuw+6+zagK+RXNc9wzMUhD0Kel9crw9373P1JoiAkx0D3fKSavqEsk9vLg09bwoJPNpentdDy0STTptBI8JkLbC953x3Squ7j7lmgF+ioc2yt9A5gX8ijsqxaZchxotgj1fQNZplU0fJpa0kX13wba7m8k3fIpFJhno/+IzeDRoJPtQnwlf+6tfY5XumN1qMmM7vezNaa2dpdu3Y1eti4opaPVHNgIMuUiuAzqTVN31C2xhGjq/BcoUzaSGvAQdNoJPh0A/NL3s8DdtTax8wywDSgp86xtdJ3A9NDHpVl1SqjIe5+l7svc/dls2bNavSw8SH8rup3VioN5/IMZvOHtXwmtWXoH0xGy6cwwKA1nSKTMoZ1z6cpNBJ8fg0sDqPQWokGEKyu2Gc1cG3YvgJ4zN09pK8MI9UWAouBp2rlGY55PORByPN7MWXIcaPTKeX6BqPWTeWAg8ltGfqGsiThV3A4W9ryMS0s2iQycTu4e9bMPgU8CqSBe9x9g5ndCqx199XA3cB9ZtZF1BpZGY7dYGYPABuBLPBJd88BVMszFHkjsMrMbgOeDnlTq4yQ10vAVKDVzC4HPujuG4/2pIxXavlIpYM1gs/E1gx5h0PDOSa2xl5GTqhCSycTWj4aat0cGvpf4+6PAI9UpN1Ssj0AXFnj2NuB2xvJM6RvJRoNV5ler4wFdb+A1FX4VdUTIKXSvv5hAKZOaClLn9wWrXjQN5iA4JMrdLuFlo/+HzcFrXAgxYEG+pWVSj19QwB0TG4tSy/cAyp0y42lbGHAQUotn2ai4CPFIdYa7SaVCsFn5qTqwedgAoJP6Wi3TDpVfC/JpuAjxWV1FHuk0p5Cy6cy+LQmp+UzEBY4ndCSpr0lxcBwLhEDIaQ+BR8pLiiq31eptLdviHTKmNpefs9nUuGeTwLm+hwKk10ntKaZ0JIm7yP3gSS5FHykSN1uUmlP3yAzJraQSpXP8Z4xMWoJ7e0bHotqlTk0FAWfia3p4hNXB7LJmIMktSn4SLGLQsFHKu3YN8BbprUfll4YgFC4JzSW+kPwaW9J01YIPglZ+kdqU/CR4vwexR6p9Oq+Q8ydPuGw9MltGVozKXb3DY5BrcoVAs2EljTt4XHfSXrQnVSn4CMjAw402FpKuDuv7j3E3OkTD/vMzOiY1Mqeg2Pf8im959Oulk/TUPCRkaHW+mNRSuztH+bQcI65Mw5v+UDU9ZaEbrfCPZ8JLdGAAxgZASfJpeAjIyscqN9NSmx+4wAAp8+aVPXzjklt7D449t1uheHeEzTgoKko+EhxwIFCj5TaFILPmW+ZWvXzU6a38+reQ6NZpap6Dw0zsTVNWyaa5wMjrSFJLgUfKXa7aWKelHrhtf1Mn9jCnKltVT8/deYk9vQNjfkqB3v7h5ke1p4rrDPXn4D5R1Kfgo8UlyPRklhSau1Lezlr3nSip9sf7rSOaCDCy3v6RrNah+k9NMS0MO9o+sQoCO3tH/v5R1Kfgo8wFIKPGj5SsPPAAJt3HuSCt9Z+Uv2pMwvBp3+0qlXVvpKWT2Hy6z4Fn8RT8JHiUiQacCAFP9m4E4ALF59Uc59FsyeTSRkbdvSOVrWq2nVwsDjpdUJrmrZMin39Yz8KT+pT8JHi808UfKTgod90s2j2ZDpPrj7YAKIVBZacPIVnto9d8Mnnndf2DZQNB58+sYW9Cj6Jp+AjImXWbN3Dupf3svK8+TXv9xScM38667fvYyg7NvNqdh0cZCiXZ17JKgwzJrbSk4A156Q+BZ83ieMxo1stHxkYzvE3qzcwZ2ob17z7tNj9LzpjNgcHs/xiy+5RqN3htu6KBjuc2jEyF+mU6RPo3ju296EknoLPmzRM9WgAABFYSURBVMBjL77Bkr/+IS+8tv+Ijy19dLZWOBjfhrJ5/nLVel58/QB/+8dnFSds1nPh4pOY3Jbh4d+8Ogo1PNzzr0ZdfktPGekePHXmRF7p6S9OHfjRhte5+AtPjPm9KSk3tg9fl+PiO2u7Adi88yBn1umjr2b/wEj3hNo949e6l3v463/YwMbX9nPLZZ28722zGzquvSXN1efP5+4nt/EXlyxm0ezJJ7im5Z56qYe50ydw0uSRuUindUykfyjHroODzJ7SznfWdbN1dx//9oFn+N9/fiEt6eT/ze3urN++j2e7ezk0nKMlneKkya3MmdrOvBkTmDt9QmyXaNIp+LwJFJY4aUkd+X/G0uVR1O02fvQNZnnx9f38amsPP9rwOs9093LS5Fbu+ti5fHDpW44or+vf81YeWNvNp+9fz6rr3118xPaJ1ntomJ/9dhcrz5tfll74A+zZ7b28v7OdzW8coDWd4sXXD/CFRzdx04eWNHzh7u0f5ltrXuZHG9/grLnT+A8fPpMJrfEtwmMxmM3x5//raX608Y2a+0xpy7Dk5Cmct2Am7z1jFueeNoNMEwTVUgo+bwKv9Q4AHNVM8y27RiYIjtcVDtydvEM2nyeXd7J5Jx9eC+9zOSfnTt49zIcaefS4F/MpfyR5ceUIvGwOVWW6hzp48bPS9PLPnOiD0vdleXj0fJtDwzkODWXpH8rRP5Rj/8AwOw8MsnP/AK/vH6B776FiPc6eN41bLutk5fnziysEHIlZU9r4wpVnc8O31rHyrl/xhSvP5m1vmXLE+Ryprz6xhcFsniuXlQefc+ZPpzWdYs22PVyw6CRe7unnLy5ezO6Dg3z9Z1vp2nmQT1y4kHNPm1G1a3E4l+fZ7n3847Ovc/+vX6FvKMfb507lW2teZsuug9zzr86L7ZJ0d/5h/as88tzr9PQN0ZZJcdLkNk6e1s7J09pZcvJUOk+ZetgTYgG++c/b+NHGN/jspW/jymXzmNyWYSibZ9eBQd7YP8jLPX28+NoBNuzo5es/28pXntjC7CltXHXefD72u6cxe8rhz19KIgWfJjecy9Md1tc6muCz9qWe4vZYrXDg7gzl8hwaytE3VH7RjNKyxe3+oRwDwzmGcnkGh/MM5XIMZfMMZvMVr4enD2bz5PL58qASft7sJrWmmTO1ndlT23jH/Blc8c75LD1lKmfNm8bsqcd+sfpA5xy+ds25/PsHn2H5f/8Z7z1jFh/onMP5C2ay4KRJx7WraziX554nt/H1n23hynPn8fa508o+b29J867TZ/LohjdYcc5c3OGMOVP4N5csZkHHJP7up5v56Ys7SaeMeTMmMH1iK22ZFIPDOXoPDbN97yFyeSedMi4762RueO9bOfPkqTz8m24+88Az/Pm3n+ar//KddVsat/3jC9z95DZOnTmReTMmMJjNs377Pn64YaBsZOCy02Zw5bJ5/PE755FJp3it9xBfeqyLS5fO4ZPvW1Tcb2IrTJ/YyuI5U7iQkblX+weGeXLzbr6zdjtferyLu5/cxr9+z1u54aLTacuc2BbasVLwaWLuzjf+eWvxfd8RBp/tPf08sLabs+ZN49nu3qNe4WA4l2dv/xA9fUP09g9zYCDLgcFh9h/KcmBgmP0D4fVQlv0D0ef7B4Y5OJCNAspw7ogDQGs6RWsmRVum/DXaTtOaTjFpUqZkvzStmRQtaSOdMjIpIxVe06lUeB35LF3yWTpF8TUVumvMjELHjRkU3kXbI+lUTa84tnB82KfweSG9eGzFe7PS7WiHlEUTLSe2ZKLX1ugxA5WPwT4RPtA5h5985r3c+8uXeWhdN09s2gVAS9qYP3Mis6e0MXtKOydNbmNKe4ZJbWkmtkavk1oztGRSJec+Vfw3yLnTN5hlX/8wa7bt4acv7OS13gEuXTqH/3z526vW5aNnn8JnH3yWe3/xEgBnzJlMKmX8P+85nY/97mn802938Wz3Pl7pOUTvoWEGh3NMn9jKvJkTueysU1hy8hR+f9Espk0caZn80TvncWAgy9+s3sCf/v1vuP3yt1cN3Gtf6uHuJ7dxzbtP5daPvr3s3Ls7uw4OsnHHftZv38cjz73GjQ89xzf+eRt/9eEz+c667eTd+Y8f6WzonE9tb+HDv3MyH/6dk9m2u4/P//BF7vzJb3l0w+v83dXvGPV7cEfCxmNXy7Jly3zt2rVjXY2jMjCc4zev7OXnXbv5ycadbHrjAB/onMPjL+7k//7907npQ0vqHn9gYJjnXu3lJxt38sDa7aRTxjc+vow/+fov+dxHl3Lt7y0o7pvN5Xmlp59X9x3itX0D7Og9xOu9A+zpiwJNT98Qew4Osn+gftBry6SYOqGFKe0ZprS3MLU9w9T2Fia3RRfIwkVoQkt0sZzYlmFi6Xa4gE5sjfZry6RG5WIqR8/d2bLrIM9297J550Fe3tPHzv2D7Do4yO4Dg/Qd5arTE1rSXLCog6vPP5WLl8yuee9mYDjH+77wBK/1DtCSNjZ8bjmtmePT+rr7yW3c8YMXAbjs7JO5atl8zlswk1TKGM7l+YP//0n2HxrmJ//2vbHdmO7Ojza+wd/+4EW27Y66wD976dvKWj1H6scb3+DGh54lm8vzjY8v412n114i6UiZ2Tp3X3Y88lLLJ6HyeeeNAwNs293HC68dYOOO/WzY0UvXzoNkQ5fAO+ZP5/NXnMUfvWMuy27/SVnL58DAMJt3HmTzGwfY/MZBfhu2C/eHWtMp3t85m89euqT4mOT9h4Z5tnsf//D0DtZs28PmnQcPmzx40uQ2TprcysxJrSw9ZSodk1qZOamNmZNamDmpjWkTWpg6IQoyUbDJJL75L8efmbFo9hQWza5+7yefd/qHc/QPZukbytE3mGU4ly/rDo1e8xjGpLYMk9syvHX2pIb+P7W3pPnPK97On/79Oi4/Z+5xCzwA1124kEuWzOaen2/joXXdPPybV5k7fQIrzjmF3kPDvPj6Ab52zbkN3T8zMy5d+hbe97bZPLium5TBn1TcwzpSH+icw5knX8DH73mKj93zFP/jX53H7y2qvUzSWFHLZwwV7te8vKePV3r6eWl3P6/09PHynn5e6elnsOTCP3tKG52nTKXz5KksWzCD8xd2MLlkVNEFf/sYMya18LY5U3nqpT1s7xl5zkpbJsWi2ZM5Y84UFs+ZzNvmTOFdp5cff8Z//AFD2Txm0JJOcd6CGSw9ZRqLZ09m/syJzJ0+gTlT24/rL7HIidZ7aJip7ZkTNiy5fyjLjze+wXeffpV/3rybXN75o3fO5b9eefaYD4Xe2zfEVXf9klf3HmLV9b/L78ybFn9QjOPZ8lHwOYGyuTyv9UYji7r39rM9vHbvPUR3Tz+v7x8ou8nf3pJiQcckTp05kdM6JnJqxyQWdExkyVumMmtK9WeqFLzn84/zSk8/U9oyvOeMWXSeMjUKNiF4pGO6qVZ8+ec8s30fn7hgIZ/+wGKmVBmFIyK17T44yPaefs6eNz0x3cKv9w7wx1/9BYPZPN/9s99jfliJ/GiNevAxs+XAfwfSwDfd/W8rPm8D/idwLrAHuMrdXwqf3QxcB+SAv3D3R+vlaWYLgVXATOA3wMfcfehoyqjleAafXN7p3tvPll0H2bqrjy27DvLS7n627+3ntd6BshvpZvCWqe3MnxGNgJk3YwLzZ05kwUmTOG3mRGZNaTvqv5b+x8+38WTXHv7LH1a/CRpn14FBDg5mWXhS9Ucmi0hz6tp5gD/6yi+YPbWdh274vbJBFEdqVIOPmaWB3wIfALqBXwNXu/vGkn3+DDjL3W8ws5XAH7r7VWbWCXwbOB84BfgJcEY4rGqeZvYA8LC7rzKzrwHPuPtXj7QMd695R/NYgs9gNscvtuzh55t380z3Pp57tZeB4ZHusRkTW1hw0qTiEMso0ETbp0yfoG4rERl1v9q6h4/dvYZzT5vBvZ84/6jvw472gIPzgS533xoKXwWsADaW7LMC+E9h+0HgSxb9Cb8CWOXug8A2M+sK+VEtTzN7AbgY+Bdhn3tDvl89ijJ+2eA5aNi6l3v4+N1P0TeUozWT4u2nTGXleafSefJUTp81idNnTWbmpNbjXayIyDF59+kd/H9XnM1f3r+e//Dw8/zXPzl7rKvUUPCZC2wved8NvKvWPu6eNbNeoCOk/6ri2Llhu1qeHcA+d89W2f9oyigys+uB68Pbg2a2BzimpXg3A989lgyS4ySO8Vy8iehcRHQeRrypzsUXgS9edVSHngTEL3XeoEaCT7WbEJV9dbX2qZVere+p3v5HU0Z5gvtdwF2F92a29ng1H5udzsUInYuIzsMInYtIOA8Ljld+jdyA6AZKB57PA3bU2sfMMsA0oKfOsbXSdwPTQx6VZR1pGSIiklCNBJ9fA4vNbKGZtQIrgdUV+6wGrg3bVwCPeTSSYTWw0szawii2xcBTtfIMxzwe8iDk+b2jLENERBIqttst3F/5FPAo0bDoe9x9g5ndCqx199XA3cB94WZ/D1EwIez3ANHghCzwycIotGp5hiJvBFaZ2W3A0yFvjqaMGHfF7zJu6FyM0LmI6DyM0LmIHNfzMC4nmYqIyNjSpBMRERl1Cj4iIjLqxmXwMbPlZrbJzLrM7Kaxrs+JYGb3mNlOM3u+JG2mmf3YzDaH1xkh3czs78L5eNbM3llyzLVh/81mdm21spLMzOab2eNm9oKZbTCzfxPSx9W5MLN2M3vKzJ4J5+FzIX2hma0J3+n+MACIMIDn/nAe1pjZgpK8bg7pm8zs0rH5RsfOzNJm9rSZfT+8H5fnwsxeMrPnzGy9ma0NaSf+9yN69O74+SEa4LAFOB1oBZ4BOse6Xifge74HeCfwfEna54GbwvZNwB1h+8PAD4jmTL0bWBPSZwJbw+uMsD1jrL/bEZ6Hk4F3hu0pRMs6dY63cxG+z+Sw3QKsCd/vAWBlSP8a8Kdh+8+Ar4XtlcD9Ybsz/M60AQvD71J6rL/fUZ6TzwD/C/h+eD8uzwXwEnBSRdoJ//0Yjy2f4nJB7j5EtIjpijGu03Hn7j8jGhVYagXRkkWE18tL0v+nR35FNNfqZOBS4Mfu3uPue4EfA8tPfO2PH3d/zd1/E7YPAC8QrYAxrs5F+D4Hw9uW8ONEy1k9GNIrz0Ph/DwIXGJWvpyVu28DSpfMahpmNg/4CPDN8N4Yp+eihhP++zEeg0+15YIOW47nTWqOu78G0UUZmB3Sa52TN9W5Ct0l7yD6q3/cnYvQzbQe2El0cdhCg8tZAaXLWTX1eQj+G/DvgcKqwA0v7cWb71w48CMzW2fRMmQwCr8f4/FJpg0txzPOHNPSRc3AzCYDDwF/6e77rfajK96058Kj+W/nmNl0omUJz6y2W3h9054HM7sM2Onu68zsokJylV3f9OciuMDdd5jZbODHZvZinX2P27kYjy2f8bwczxuhiUx43RnSj3QZpKZiZi1Egefv3f3hkDwuzwWAu+8DniDqsx+Py1ldAHzUzF4i6na/mKglNB7PBe6+I7zuJPqj5HxG4fdjPAafRpYLerMqXaKocumij4eRLO8GekNT+1Hgg2Y2I4x2+WBIaxqhb/5u4AV3/2LJR+PqXJjZrNDiwcwmAO8nuv817pazcveb3X2eR4tkriT6bv+ScXguzGySmU0pbBP9v36e0fj9GOuRFmPxQzRi47dEfd5/Ndb1OUHf8dvAa8Aw0V8l1xH1U/+U6GkQPwVmhn0N+HI4H88By0ry+QTRjdQu4P8a6+91FOfhQqLm/7PA+vDz4fF2LoCziJarejZcXG4J6acTXTC7gO8AbSG9PbzvCp+fXpLXX4Xzswn40Fh/t2M8LxcxMtpt3J2L8J2fCT8bCtfD0fj90PI6IiIy6sZjt5uIiIwxBR8RERl1Cj4iIjLqFHxERGTUKfiIiMioU/AREZFRp+AjIiKj7v8AY0s9ytkWjLwAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# zfit.run.numeric_checks = False \n",
"\n",
"fitting_range = 'cut'\n",
"\n",
"Ctt_list = []\n",
"Ctt_error_list = []\n",
"\n",
"nr_of_toys = 10\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",
"# clear_output(wait=True)\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",
" reset_param_values()\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",
"\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",
" if fitting_range == 'full':\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",
" \n",
" probs = total_f.pdf(test_q, norm_range=False)\n",
" calcs_test = zfit.run(probs)\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_full_range{}.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",
" \n",
" if fitting_range == 'cut':\n",
" \n",
" tot_sam_1 = np.where((total_samp >= x_min) & (total_samp <= (jpsi_mass - 50.)))\n",
" \n",
" tot_sam_2 = np.where((total_samp >= (jpsi_mass + 50.)) & (total_samp <= (psi2s_mass - 50.)))\n",
"\n",
" tot_sam_3 = np.where((total_samp >= (psi2s_mass + 50.)) & (total_samp <= x_max))\n",
"\n",
" tot_sam = np.append(tot_sam_1, tot_sam_2)\n",
" tot_sam = np.append(tot_sam, tot_sam_3)\n",
" \n",
" data = zfit.data.Data.from_numpy(array=tot_sam[:int(nevents)], obs=obs_fit)\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_fit.get_dependents():\n",
" param.randomize()\n",
"\n",
" nll = zfit.loss.UnbinnedNLL(model=total_f_fit, data=data, constraints = constraints)\n",
"\n",
" minimizer = zfit.minimize.MinuitMinimizer(verbosity = 5)\n",
" # minimizer._use_tfgrad = False\n",
" result = minimizer.minimize(nll)\n",
"\n",
" print(\"Function minimum:\", result.fmin)\n",
" print(\"Hesse errors:\", result.hesse())\n",
"\n",
" params = result.params\n",
" \n",
" if result.converged:\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",
" \n",
" probs = total_f_fit.pdf(test_q, norm_range=False)\n",
" calcs_test = zfit.run(probs)\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_cut_region{}.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('{0}/{1} fits converged'.format(len(Ctt_list), nr_of_toys))\n",
"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": [
"# sample from original values"
]
}
],
"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
}
