{
"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=-1.0, upper_limit=1.0)"
]
},
{
"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": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<tf.Tensor 'normalization/hook_integrate/hook_numeric_integrate/mul_1:0' shape=() dtype=float64>"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"total_f_fit.normalization(obs_toy)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Test if graphs actually work and compute values"
]
},
{
"cell_type": "code",
"execution_count": 12,
"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": 13,
"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/z4zo9FqSdbmTcaSbWFjgzHg2BBDsM1mEhqTBoLpA6EJfWjS0DTL04R0SZ+m0CekKaRtSIgLNGQphpLkh0PYYwMBjBdWr7JleZM3rZZkSSPNcn9/fL8zGo9nlSXNaHTer5de/s793nvuna+l+cy599xzxRiDoiiKomQCjnQPQFEURVGCqCgpiqIoGYOKkqIoipIxqCgpiqIoGYOKkqIoipIxuNI9gEyjoqLC1NTUpHsYipLx7DnRTV6Ok3PKCk4r7+zzcqi9l7qqIvJynGkancXe5lO4nQ5mlBfgCxh2Hetiamk+5YXutI4rG3nnnXdajTGVZ2tHRSmCmpoatm7dmu5hKErGs+xfNrCgupR/v/Wi08pf2H6cL/ziHZ788uXMn1qSptFZrPzB60wvK+A/P7uIjp4BLvqnl/n2H83jc0tr0zqubEREDg6HHZ2+UxRlSAQMOB1yRnmwLBAY7RFFJzhChz0uf0D3ZmYyKkqKogwJf8AgZ2oSruCHfwZszDeG0BhDYpkB41Jio6KkKMqQMMbgjKJKgx5J+l0lg0FsXyk4Vn/6h6XEQdeUFEUZEn5jcEQRpUz68A/3lBz2V/DR9pS8Xi9NTU14PJ5R7XekyMvLo7q6mpycnBGxr6KkKMqQCJhBrygcZwat3RjCpu8kPeNqampiwoQJ1NTUINHmO8cQxhja2tpoamqitnZkgkV0+k5RlCERCBiiaFJmiZIxISFI17g8Hg/l5eVjXpAARITy8vIR9fpUlBRFGRJ+Y2JE3w3eTzfGDEbfiQgi6Ql0yAZBCjLS70VFSVGUIWF5StFEyRG6n26s6bvBMTpFMsKDU2KjoqQoypAIGOIGOvgy4MPfGEP4CJ0OyQgPLhN59dVXueGGGwDo7+/n6quvZuHChTz55JOjOg4NdFAUZUgEjAlN1YUTjHLLBI8kPNABLFHKBA8u03nvvffwer28//77o953Up6SiKwUkXoRaRCRe6LczxWRJ+37m0SkJuzet+zyehG5LpFNEam1bey1bbrj9SEi5SKyQUROicgPY4x/nYhsT+6RKIqSDP6Y03eZs0k1fE0JgtN3aRtO2jhw4ABz587ljjvuYMGCBdx000309vbywgsvMHfuXC6//HJ+/etfA9Dc3Mxtt93G+++/z8KFC9m3b9+ojjWhpyQiTuAh4BqgCdgiIuuMMTvDqt0JdBhjZovIauB+4BYRmQesBuYDU4FXRORcu00sm/cDDxpj1orIw7btH8fqA/AAfw+cb/9Ejv+PgVMpPRVFURJiYoSEBzM6ZMT0Hea0NSWHQ9Iqlv/42x3sPNo1rDbnTS3mH/5ofsJ69fX1PProoyxdupTPf/7zPPDAA/zkJz9h/fr1zJ49m1tuuQWAqqoqHnnkEb7//e/z7LPPDutYkyEZT2kx0GCMaTTGDABrgVURdVYBj9vXTwNXifWbsApYa4zpN8bsBxpse1Ft2m1W2Dawbd4Yrw9jTI8x5g0scToNESkCvgbcm8T7VBQlBazNs2eWB72nTJgmC988C/aaUgaMKx1Mnz6dpUuXAnDbbbexdetWamtrqaurQ0S47bbb0jxCi2TWlKYBh8NeNwFLYtUxxvhEpBMot8vfjmg7zb6OZrMcOGmM8UWpH6uP1jhj/yfgX4HeeG9QRO4C7gI455xz4lVVFMXGH4ieZiiz9ikRSjMElmCmM9AhGY9mpIgM5e7s7MzIUPVkPKVoo478X41VZ7jKkx3H4IBEFgKzjTG/iVUnZMSYNcaYRcaYRZWVZ30ciKJkPcb+YI+b0SEj1pRMhKcEfn/6x5UODh06xMaNGwF44oknuPrqq9m/f39ozeiJJ55I5/BCJCNKTcD0sNfVwNFYdUTEBZQA7XHaxipvBUptG5F9xeojFpcBl4jIAeAN4FwReTXuO1UUJSmCXlC8QIeM8JSIEuiQAWKZDs477zwef/xxFixYQHt7O1/96ldZs2YNn/jEJ7j88suZMWNGuocIJDd9twWoE5Fa4AhW4MKfRNRZB9wBbARuAtYbY4yIrAP+W0QewAp0qAM2Y/2enGHTbrPBtrHWtvlMvD5iDdoY82OsAAnsSL1njTHLkni/iqIkIKg3UTM6pCnHXDQi15Qc4zgk3OFw8PDDD59WtnLlSnbv3n1G3WXLlrFs2bJRGtnpJBQle/3mbuBFwAk8ZozZISLfAbYaY9YBjwI/F5EGLO9ltd12h4g8BewEfMCXjDF+gGg27S6/CawVkXuB92zbxOrDtnUAKAbcInIjcG1EdKCiKMNIMIIt2pJERoWEhx1dAbp5diyQ1OZZY8xzwHMRZd8Ou/YAN8doex9wXzI27fJGrOi8yPJ4fdQkGP8BooSLK4oyNIKCEy/QwZcBazdnRN+N0zRDNTU1bN8+NrZqapohRVFSJt6akiOjPKUoGR3SMK44Kw1jjpF+LypKiqKkTPBQ2XibZzPBI7E8pYjpu1EeV15eHm1tbVkhTMHzlPLy8kasD819pyhKygTXZVxRRMmRwQlZHWlIM1RdXU1TUxMtLS2j2/EIETx5dqRQUVIUJWV8tqsU/TylDMroQPqn73JyckbslNZsRKfvFEVJmeAUWNyQ8AyYrrI8pdNz32XCtKISGxUlRVFSJp4oORzWCa+Z8OF/hqeUppNnleRRUVIUJWVCohQjd1qmhF6fcXSFekoZj4qSoigpE/xgdzmji5IjQzapWrnvIhKyqihlNCpKiqKkTLzpO7Ci8jTQQRkKKkqKoqSMP05Gh2B5ZoSEc2aaoQwYlxIbFSVFUVImmEIolqeUKYlPI4+usM5TSt94lMSoKCmKkjKh3Hdxpu8yYk2J0wMdXA7Bb++x6uzz8tqe7NjQmk2oKCmKkjK+BGtKmbIfKNpx6EEv76//5wPueGwzLd39aRqdEg0VJUVRUiYUfeeI/hGSMSHhnB59l+N0hAT1w6ZOwPKYlMxBRUlRlJQJZQmP8QnidGRSoMMgLqfgs5Pf5eZYg+/2qChlEipKiqKkTCJPyeXMFE+J01TJ5XDg9Z8eOdjvG+UMrUpcVJQURUmZwX1K0e+7wtZu0okx5rQzn3KcEkomGxSrARWljEJFSVGUlBkUpegfITlOB97RPiMiCtGn704/oFBFKbNQUVIUJWUS5b7LGFHi9Og7a/rOGleweCADxqkMoqKkKErKJAoJdzkzJdDh9KMr3K7BNSXR6buMJClREpGVIlIvIg0ick+U+7ki8qR9f5OI1ITd+5ZdXi8i1yWyKSK1to29tk13vD5EpFxENojIKRH5YZidAhH5nYjsFpEdIvLd1B+PoijRSJSQNceRqZ7S4JqSTt9lJglFSUScwEPA9cA84FYRmRdR7U6gwxgzG3gQuN9uOw9YDcwHVgI/EhFnApv3Aw8aY+qADtt2zD4AD/D3wP+JMvzvG2PmAhcBS0Xk+kTvV1GUxASzNThiTN+Fr92kkzPXlCxPyYRlm+jPAPFUBknGU1oMNBhjGo0xA8BaYFVEnVXA4/b108BVYu1YWwWsNcb0G2P2Aw22vag27TYrbBvYNm+M14cxpscY8waWOIUwxvQaYzbY1wPAu8DIHSyvKOOIYKoeV8zpOwfeDJi+A05zlXLs8foDg5tq+73+tAxLiU4yojQNOBz2uskui1rHGOMDOoHyOG1jlZcDJ20bkX3F6iMhIlIK/BHw+xj37xKRrSKytaVFc2EpSiKCzkWsNaUcx+Am1XRhQt7cYJnLjmH3BUyoXAMdMotkRCnab13kV6BYdYarPNlxnIGIuIAngH83xjRGq2OMWWOMWWSMWVRZWZnIpKKMe4KeUtxAhzRP3wUdtfBAhxx7DczrD4QcKK8vQzw6BUhOlJqA6WGvq4GjserYIlACtMdpG6u8FSi1bUT2FauPRKwB9hpjfpBEXUVRkiDoXMSfvssMTyky0AGsozeC72HAr9N3mUQyorQFqLOj4txYgQvrIuqsA+6wr28C1hvrN2IdsNqOnKsF6oDNsWzabTbYNrBtPpOgj5iIyL1Y4vWVJN6noihJEvSUHHGn79LrgUSbYslxWR95Xn8gNL2o0XeZhStRBWOMT0TuBl4EnMBjxpgdIvIdYKsxZh3wKPBzEWnA8l5W2213iMhTwE7AB3zJGOMHiGbT7vKbwFpbUN6zbROrD9vWAaAYcIvIjcC1QBfwt8Bu4F17UfOHxphHUn9MiqKE4wvlvovtKaV/Tcn6N9xTyrEzUHgDJhSyrqKUWSQUJQBjzHPAcxFl3w679gA3x2h7H3BfMjbt8kas6LzI8nh91MQYevS/GEVRzorBLOGxMjoIA2n3lILTd4NjDO6r8vkDoU206R6ncjqa0UFRlJTxJ/CUrHOLMsNTCicYfef1q6eUqagoKYqSMv6Ex6E70r6mFMQRZZ+SLxAYFCUNCc8oVJQURUkZf8SZRJHkOCXtaYYC0aLvgvuU/CYkml71lDIKFSVFUVImoaeUAQlZQ4EOYWWusH1KA+opZSQqSoqipIzfzoggsXLfORz4A6fnmBttQiHhYUN065pSxqOipChKyvgCJuZR6BCeOSGNohQlaWwwMGPAFwhlfFBRyixUlBRFSZlAwBBHk8JyzKXvAz+UZui0kHBrXH1hSVh1+i6zUFFSFCVlvH4T2ogajaBHkhme0mBZ0IM7TZTUU8ooVJQURUkZrz8QStkTjZxQlFv6PaXTp+9sT2nAFypTTymzUFFSFCVlfIFAzI2zEB7llj5PKRDFU8rNsT7yevoHPaV0h64rp6OipChKynj9JuQNRSPHOZj4NF0M7lMaVKVc27s71R/mKen0XUahoqQoSsp4/YHQ+kw0gvfSuVfJRJm+y8txAoOilJfjUFHKMFSUlJjsPdHNd367M617TZTMxOc3oUi2aATXbtK5phRKGhs+fWd7St0eS5QK3S4VpQxDRUmJyR2PbeaxN/dztNOT7qEoGYbXH39NKRP2KYXWlBxneko9tqdUkOvUQIcMQ0VJiUm8c+iV8Y3XH8AdJ/ou5CmlcZ9StOk7l0NwyOD0XaHbxYA/oLMBGYSKkpIQ/XNVIrEyOsTxlGzBSufUWLToOxEhL8cZEqV8txNj0rv2pZyOipKiKCkz4AvEXVPKzQhRsv51ROTny3U5OBW2pgQaFp5JqCgpCdHpOyUSX8CEkptGIyhK/RkREn56ea7LSY+9ebbAba0xabBD5qCipCREJzaUSHz+QGiDbDSC60393nSuKZ2ZkBWsMPCQp5RreUoqSpmDipKiKCkz4I+fJTzXZXkg/T5/zDojTezpOyfdYWtKAP0qShlDUqIkIitFpF5EGkTknij3c0XkSfv+JhGpCbv3Lbu8XkSuS2RTRGptG3ttm+54fYhIuYhsEJFTIvLDiHFdIiLb7Db/LrEOf1Hiog9NicTnD+B2xf7NCE3fZVigA5y+YbbQFiVdU8ocEoqSiDiBh4DrgXnArSIyL6LanUCHMWY28CBwv912HrAamA+sBH4kIs4ENu8HHjTG1AEdtu2YfQAe4O+B/xNl+D8G7gLq7J+Vid6vciY6fadEkug8pWCOubQGOthdR34XDXpxAAV2oIPuVcockvGUFgMNxphGY8wAsBZYFVFnFfC4ff00cJXtlawC1hpj+o0x+4EG215Um3abFbYNbJs3xuvDGNNjjHkDS5xCiMgUoNgYs9FYk8s/C7OlKMpZYEXfxfGUnOmfFovlKQUFE6AwVwMdMo1kRGkacDjsdZNdFrWOMcYHdALlcdrGKi8HTto2IvuK1Ue8cTclGDcAInKXiGwVka0tLS1xTCqKAtam2LjRdznB6bv0rSlF2zwLp3tKE/JyABWlTCIZUYr2dShyRidWneEqT3YcyYzpzEJj1hhjFhljFlVWVsYxOb7QtSQlFlbuuzjRd84MmL4LpRk6vTwY3ABQHBSlUZi+e7W+mc8+tpm+gfQJ9VggGVFqAqaHva4GjsaqIyIuoARoj9M2VnkrUGrbiOwrVh/xxl2dYNyKogyBAX8g7pqSwyHkOCUjpu8i15Qm5LlC18X5oxcS/v2X6nl9TwubD8T72FKSEaUtQJ0dFefGClxYF1FnHXCHfX0TsN5ex1kHrLYj52qxgg02x7Jpt9lg28C2+UyCPqJijDkGdIvIpfZa1WfDbCmKchb4/Cbu0RVgTZOlc59SrJDwoHcEozt9d9xObNzYcmrE+xrLuBJVMMb4RORu4EXACTxmjNkhIt8Bthpj1gGPAj8XkQYs72W13XaHiDwF7AR8wJeMMX6AaDbtLr8JrBWRe4H3bNvE6sO2dQAoBtwiciNwrTFmJ/BF4KdAPvC8/aMoylniCwTiHvIHVlh4eteUogc6hHtK+TnBkPCRjTE1xoSOyzjS0TeifY11EooSgDHmOeC5iLJvh117gJtjtL0PuC8Zm3Z5I1Z0XmR5vD5qYpRvBc6Pdk9RlKFhjMGb4DwlsEQpE3PfFecPekrBzBMD/pEVzz6vPzSVebRTRSkemtFBUZSUCGbUzomTJRysD/zMWFM6vbw4zFMarXRInX3e0HV7z8CI9jXWUVFSFCUlgtkPcuKcpwT2mlJa0wxFz30XvqZUYE/f9XlHdpxdfb7QdUePN05NRUVJSYgegKaEE1x/iXeeElh7ldI5fRdrn1JZoTt0XWBvnu0d4TDtoKdUPTGf9l71lOKhoqQoSkoEvZ/cHGfcem5nZkzfRWrnrKoiAGrKC3A7HbgcEjoefaQIilJtRSEdPQP6RS8OSQU6KOMbzWOrhBP0fnITTd/lOPBkQEh45O9vUa6Le288n4XTSxERCtzOUfOUaisK+cPeVro8PkrCAi6UQVSUlITotzolnP5kRcnlPG2Bf7SJ5SkB3HbpjNB1Ya6L3oGR9ZS67OdQU14IQEfPgIpSDHT6TlGUlAhGqiUWJUdGHvIXSb7bSc8oeUozygsAaNMIvJioKCmKkhKhNSVX/DWlvBwnnnRG39l6mEiUCt0uekdhTWlCrouKolwATmqwQ0xUlJSE6OydEk6y03f5bie9/ekPCU+0JDoaa0pdfV6K83NCkX+6Vyk2KkqKoqREKNAhJ/7HR0HOyH/YxyNWRodIrDWlERYlj5eS/Bwm2qLUoZ5STFSUFEVJiUFPKf70XUGuiz6vn0AgPa62iXF0RSTWmtLIT9+V5OdQ6Hbidjro6NUNtLFQUVIUJSWCa0ruBNN3Bfa5RelaV0raUxqFacauPh/F+S5EhNKCHDp0+i4mKkpKQnRNSQkn2ei7oCj1pGldKV5IeDgFbteIe0pdHm8ovVFZoVvXlOKgoqQkxMQ94FcZbwRPaU04fee2tkGm66TVWIf8RVKSn0O3x4d/BKcZuz2+0NlNEwvcuqYUBxUlRVFSot8bDAlPzlPq9Y6sFxKLWAlZIyktsMSia4Q2+vr8AU71+0Kn3KqnFB8VJSUhOn2nhBMMdEi0ppSf5uk726FLmDg2KEonR0iUTtl7oILTdxMLczTQIQ4qSoqipESy+5RCx0KkafrOb++edSYUpZEN0w6eOBs88baswM3J3oERnS4cy6goKQnRPx0lnH6fH6dDEp48W5hrfQiPdF65WAQPI0woSnYOus4R8l6CKYaCJ96WFrgJmJGbLhzrqCgpCdGErEo4A75AQi8JBqfv0rWBNpCsKNme0sm+UfKUglkdNNghKipKiqKkRH+SolSQZlEKekoJ15RsD2akToTt8tieUmhNyRZBFaWoJCVKIrJSROpFpEFE7olyP1dEnrTvbxKRmrB737LL60XkukQ2RaTWtrHXtuk+iz6+KiI7RGS7iDwhInmpPR4FdPpOOZ1+byBhkAMMhoSna/ouuGbjSCBKxfk5iIxcoENwmi60T6kgmP9Op++ikfA3S0ScwEPA9cA84FYRmRdR7U6gwxgzG3gQuN9uOw9YDcwHVgI/EhFnApv3Aw8aY+qADtv2UPqYBnwZWGSMOR9w2vWUFNHZOyWcXq+f/ASnzsKgp5S+QIfkPCWnQ5hY4Kb1VP+IjCM4fRcMCZ9YGPTM1FOKRjKe0mKgwRjTaIwZANYCqyLqrAIet6+fBq4Sa8faKmCtMabfGLMfaLDtRbVpt1lh28C2eeMQ+wDrEMN8EXEBBcDRJN6voihx6Bvwke9OfD5ojtOB2+UIhUSPNsHpu0T7lACqJuRyotMzIuMITt8V2YEfEwt0TSkeyYjSNOBw2OsmuyxqHWOMD+gEyuO0jVVeDpy0bUT2lVIfxpgjwPeBQ8AxoNMY81K0Nygid4nIVhHZ2tLSEvNBjF/UVVIG6R3wh7ygRBTn5YQ+lEebZD0lgMkleRzvGiFR6vNR6HaGohUL3E7cLod6SjFIRpSi/Y9GfkrFqjNc5Sn3ISITsbyoWmAqUCgit0WpizFmjTFmkTFmUWVlZbQqiqLYpCRK+S66+tK7ppQo+g5gcnEeJ0ZIlDp6B0LBDWClPSor0KwOsUhGlJqA6WGvqzlzGixUx54qKwHa47SNVd4KlNo2IvtKtY+rgf3GmBZjjBf4NfDRJN6vEoGuKSnh9A0kt6YEVl65dHpKDkmc+w5gUnEeracG8PqH//j2tp4BysNECayw8JFawxrrJCNKW4A6OyrOjRUssC6izjrgDvv6JmC9sTa3rANW25FztUAdsDmWTbvNBtsGts1nhtjHIeBSESmw156uAnYl91iUcFSTlHB6vb7QxthEFOflpG2TqC9gcCU6TMlmcokVmNvcPfxC0d7TH9qbFGRKSR7Hu1SUopHwf8xev7kbeBHrQ/0pY8wOEfmOiHzSrvYoUC4iDcDXgHvstjuAp4CdwAvAl4wx/lg2bVvfBL5m2yq3bQ+lj01YARHvAtvs97pmCM9IUZQw+gb8oY2xiSjOz6HLk76ErMlM3YElEgBHOvqGfRztpwYoK8w9rWxySR7HO4e/r2wgqa87xpjngOciyr4ddu0Bbo7R9j7gvmRs2uWNDEbPhZcPpY9/AP4hWhsleXT6Tgmnd8AfymuXiOI8V/o8JX/yolRbUQjAgdYeFteWDdsYjDHW9F3RmZ5SR68Xj9dPXpLPcrygGR0URUkaYwx93lQCHaw1pXSkqvIHAkmLUvXEAnKcQmNrz5D6inXke5/XT78vcMb03eSSfACOj1AY+lhGRUlJiB7ypwTxeAMYQ1L7lMBaU/L6DR7v8AcQJMJvTFLh4GBF6M0oL2R/66mU+9lzopsLv/MS//zcmUvWrd1WhF20NSWAYypKZ6CipCREp++UIMGUQamEhANpicDzB0zCFEPh1FYUsn8IntLPNx6k2+NjzeuNtEVE1B21142mleafVh4MrDjepetKkagoKYqSNMHkqkkHOuSN7Kmu8fAHkveUAGZXFbG/tYd+X2ppkXYe66LETur6/Pbjp907etISnamRolSsnlIsVJSUhKinpATps49CT9ZTCp7qmo6TVn2B5AMdAC6YVoLXb9h9rDvpNsYYGppPccOCKdRWFPLijuiiFJyuC1KY66Ks0M3h9t6k+xovqCgpipI0PXYeu8Ik15TK7VDoyGmt0cCfoigtqC4B4MMjnUm3aT01QGefl9lVRVw7fxIb97WFDvUDOHKyj/JCd9QIu6FOF2Y7KkpKQjTQQQkSeWBdIirsUOh0ZC9IVZSmleZTVujmw8Mnk27T0GwFRsyuKuK6+ZPxBQwbdjefdn9mZWHUtipK0VFRUhKi03dKkEFRykmqflmhGxHLoxhtUl1TEhEuml7KlgPtSbdpaBkUpYXVpVRNyA1N4RljqD/ezbmTJkRtW1tRyImu/pD3qVioKCmKkjShU1Tzk/OUXE7HiJ5VFA9fwCR1bEU4l9dVcKCtN+m1nn3Npyh0O5lcnIfDIVwzbxKv1rfg8fpp6uijy+Nj7pTiqG1DG3bb1FsKR0VJUZSk6bZFKVlPCawpvLSIkj+Ay5maKF1RZ50S8Ie9rUnVb2g+xayqolDS1+vmT6bP6+f1PS28tc+ycWmMDBFBUWpsUVEKR0VJiYnO2imRdPX5cAgUJhl9B1awQ1sapu+8foPbmdpH3KzKQqaV5vPKrhNJ1W9oPsXsyqLQ60tnljO5OI9/+/1efvH2Iaon5jO7qihq25mVhbgcwq5jXSmNMdtRUVISomtKSpBuj5cJeTlJHQcRpGJCblo8pQF/ALcrtY84EeGGBVN4fU9LwvOOuj1ejnd5mBUmOm6Xg7+74Tx2HO1i25FOvnDlrJjPKtfl5NxJE9h+VEUpnOQmhpVxSfBPSaPvlCBdHl/S60lBKorctIzAkRCJGPAFko4SDGfVwmn85PVGfrftGLdfOiNmvWDkXV2EJ3TDgqmUF+Zyqt/H1edVxe1r/tRi1u9uxhiTktBnM+opKTFRKVIi6fZ4mZCb/HoSWBtHewb8o55qyOsPpDx9B3DelAnMnTyBX2w8GDeR7N4TlihFi667bFY518yblFBozp9WQlvPACcSnK3U0+/j+y/Wc2wcHHehoqQkRKfvlCBdHl/K3se00gJgZM4qiod3CNN3YE3h/dkVM6k/0c3rcQIe9pzoJtflYHpZwZDHeP40a8Pue4c64tb71btN/HBDA/f9LvvPKVVRUmKiYqRE0tXnTSnyDmDaRCvv22iL0oAvQM4QPCWAT1441QpYeGVPTG9pT/MpZlcVpbRBN5ILppWQn+Nk0/74e6PeO2Rt6N2WQraJsYqKkpIQ1SYlSEfvAGWFKYqSnYz0yMnR9pTMkEXJ7XLw1WvqePfQSX774bEz7gcChg+bTjIvxh6kVPpZVDORtxvb4tYLnvN0qL0Xjze1hLFjDRUlJSHpOKBNyTyMMbT3nHm0dyIqitzkuhyjLkpDib4L56ZLpjN/ajH/9OzOM3L37Ws5xcleLx8ZhlNqL51Zzu7j3THzAxpj2N9yitKCHIwh61MTqSgpMdGoOyWcLo8Pr99QHnFgXSJEhGml+WmZvnOnuHk2HKdD+JebLqSz18tfP/0h/rDTZV+tbwHgspnlZx8AmiQAACAASURBVD3Oj86ybLy+tyXq/faeAbo8Pq4+bxJgCWI2o6KkJESlSQFC+3bKi1ITJYDqsoJRT6fj9Q99TSnIvKnF/P0N57F+dzN/+5tt+AMGf8Dw1NbDLKguOasghyAX2jnzXtoRfcNu0DNaPscKLz+gnhKIyEoRqReRBhG5J8r9XBF50r6/SURqwu59yy6vF5HrEtkUkVrbxl7bpvss+igVkadFZLeI7BKRy1J7POMbnbVTwmnvsaaXIo/2Toa6qiL2tZwiEBi9X6qhRt9FcvtlNdy9fDZrtxzm5off4stPvMfe5lP87ytmDsMoOSNnXiRBUZo3tZhJxbnsb83uM5gS/o+JiBN4CLgemAfcKiLzIqrdCXQYY2YDDwL3223nAauB+cBK4Eci4kxg837gQWNMHdBh2065D7vNvwEvGGPmAhcC2R9POQKoOClAKFVQeYprSmCJkscbGLV1pUDAnFWgQyRfv/Zc/vXmCznW6eGlncf5wpWzuGHBlGGxDbDyfCtnXnBaMJxD7b04HdYUaE15YdYncE3mf2wx0GCMaTTGDABrgVURdVYBj9vXTwNXibVrbBWw1hjTb4zZDzTY9qLatNussG1g27xxKH2ISDHwMeBRAGPMgDEm+YNSlDBUlZTB6buyIUzf1U2ysh7sbU7+VNezwRsIAAyLpwTWutinL6nmrXtWsPufruee6+cOawaGy2aWUzUhl//ZeviMewfbeplamofb5WBmZaFO3wHTgPAn1WSXRa1jjPEBnUB5nLaxysuBk7aNyL5S7WMm0AL8l4i8JyKPiEjU07ZE5C4R2SoiW1taoi82jkdUipRw2oJrSkOYvptdZWU92HNidBbpB3y2KA2TpxRERM5qX1IsXE4HNy+qZkN98xlZGw629TCjzProqikvpK1n4LTTbbONZP7Hov0PRH5exaozXOVD6cMFXAz82BhzEdADnLEeBmCMWWOMWWSMWVRZWRmtyrhGp+8UgOOdHkryc6Ie7Z2IkvwcppTksWOUko/22WszeSlkM083n1k0nYCBJzYPfsc2xtDY2kNNhRVQURM8gymLvaVkRKkJmB72uho4GquOiLiAEqA9TttY5a1AqW0jsq+h9NFkjNlklz+NJVJKkqgYKeEc6/QwpSRvyO0vOqeU9w/HT6czXHgGLE8pfwgCmi5mlBdy9XmTePytA6Fzqw6399Ht8TFvipWOaOY4OBgwGVHaAtTZUXFurKCCdRF11gF32Nc3AeuNteNyHbDajpyrBeqAzbFs2m022DawbT4zlD6MMceBwyIyx25zFbAziferRKDapAAc6+w7O1GaPpHD7X2jkjG812utABSMIU8J4C9XzKazz8vPNh4E4MMj1jL4/KlW5ojpZQWIZPfBgAlFyV6/uRt4ESt67SljzA4R+Y6IfNKu9ihQLiINwNewp8mMMTuAp7DE4AXgS8YYfyybtq1vAl+zbZXbtlPuw27zl8AvReRDYCHwz6k+IEU9JsXieKeHKXbKoKFw8YxSAN5NkHx0OOgbsD4CxpKnBHDh9FKuPq+KH65v4HB7L6/VtzAhz8U8W5TycpxMLcnPak8pqXS/xpjngOciyr4ddu0Bbo7R9j7gvmRs2uWNWNF5keVD6eN9YFG0NkoyqBopFh6vn7aeAaaehac0f2oJeTkONu5r47r5k4dxdGcSEqUx5ikB/OOq87n2gde4/dFNHO308MkLp54W2p7tEXia0UFJiOa+U453egCYXDJ0Tykvx8lHZ1Xwan3zcA0rJsFAh7HmKYGVwHbNZxfR7fExuTiPr15z7mn3a8oLaWztydq/Sz15VomDFdiYnb/6SiocareyCFRPHLooASyfU8n63c00tpxiZmVR4gZDpNf2lMbamlKQpbMr2Pp3VwOcsR+qpqKQbo+P9p4ByotS38ic6ainpMRB5UixCKa6mVkZdatf0iyfa+Vv+12U4yCGk1BI+Bj0lIKISNQNutkegaeipCQkS2cJlBTY39pDUa6LyrP8Zl49sYBLZ5bx9LtNIzr9NJbXlBIR3KuUrTnwVJSUmKgYKUH2tZxiZmXhsKTWufmS6Rxs6+WtffEPtjsbuuyMB8UpnpI7FqiemI/b6WDvidFJ2TTaqCgpCdFzlZT9rT3UVpzd1F2QTyyYQkVRLg9taBgWe9Ho7POSn+Mcttx3mUSO08HcKRP4sCk7j0bPvv8xZfhRTRrXdHm8NHX0UVc1PIEJeTlOvnDlTN7a18Zb+1qHxWYknX1eSguyz0sKsqC6hG1HOkf1KJDRQkVJiUn2/borQ2Gnna/u/Gklw2bzfy2ZwfSyfP7uN9ujniF0tpzs81KSn82iVMqpfh+Nrdl3Cq2KkpIQFafxzfYj1jTR/KnDJ0r5bif//KkLaGzt4R9/u3PYgx46+7wUZ7EoXVhtZcf44HD2TeGpKCkxydbNeUpqbD/SyeTiPConDO+emCvqKvnislk8sfkQP3m9cVhtd/Z6Kc1iUZpdVUSh2zkqKZtGG908qyREtWl8s/VgBxdOHz4vKZy/vnYOh9p7+e7zu+n2ePnaNXOG5byi410elswsG4YRZiZOh7CopoxN+9vTPZRhRz0lJSEafTd+OdzeS1NHH5fNLB8R+w6H8G+3LOSWRdN5aMM+bn90E4fbz27/Te+Aj84+L5PPIk/fWOCyWeU0NJ+iuduT7qEMKypKSkxUipSN9l6ij86uGLE+XE4H3/30BXzv0wt4//BJrnrgNb77/G7aTg3tiItgnr6zOWZjLHCp/UXh7cbs8pZUlJSE6PTd+OW1vS1UFOUOWzh4LESEz3xkOuu/vowbLpjCw6/t46PfXc/f/X/baGhOLcLsYJvlaU0rLRiJoWYM508tpijXxduNI7cJOR3ompISExWj8U3fgJ8Nu5v51EXThiWTQzJMLsnjgVsW8hfLZ/Gfr+/nqS1N/OLtQ1w6s4zbLp3BtfMmJ9wQu/u4lelgzqQJozHktOFyOlhcW8ZbDSOz1ytdqKekJES1aXzy2p5megf8fPyCKaPe9+yqCdx/0wLevGcFf33dHA6393H3f7/HR7+7nu+/WM+Rk30x224/2snUkjxKsnjzbJCP1VVwoK03q85XUlFSEqKh4eOTp7Y2UVGUy5La9EWxVU7I5UvLZ/P6N5bz2J8u4sLqEh56tYEr7l/Pnz2+hQ27m/GHZTXw+QO82dDKpbNGJjAj01gxdxLAqJxRNVro9J0SExWj8cvh9l421Dfzl8tn43Km/7ur0yGsmDuJFXMn0dTRyxObD/HklsO8squZ6on5/MmSc/jMoum8tOMEJ3u9fCIN3l06OKe8gJmVhayvb+FPl9amezjDgoqSkhCVpvHHI39oxCHCrUvOSfdQzqB6YgF/fd1c/uqqc3lxx3F+uekg33uhnn95sR5jYHFNGcvnVKV7mKPG8jlV/Pztg/QO+Chwj/2P9LH/DpQRQ8VofHK4vZf/3nyIWz4ynSlncfz5SON2OfijC6fyRxdOpaG5m1+/e4QJeTncftkMHMOwAXessHxOFY++sZ+N+9q46rxJ6R7OWZOUXy4iK0WkXkQaROSeKPdzReRJ+/4mEakJu/ctu7xeRK5LZFNEam0be22b7qH2Yd9zish7IvJs8o9FOQ1Vp3GDMYZ7f7cThwhfXlGX7uEkzeyqCXxj5Vy+uGwWRbnj67v2R2onUuB2siFL1pUSipKIOIGHgOuBecCtIjIvotqdQIcxZjbwIHC/3XYesBqYD6wEfmSLRDyb9wMPGmPqgA7bdsp9hI3tr4BdyT0OJRqa0WH88OyHx3hxxwm+es25WZ8RIVvIdTm5fHYFG3a3ZMU6cDKe0mKgwRjTaIwZANYCqyLqrAIet6+fBq4Sa2PDKmCtMabfGLMfaLDtRbVpt1lh28C2eeMQ+0BEqoFPAI8k9ziU0xj7v99KCjQ0d/M3v97GhdUl/Nnl2bFoPl5YMbeKIyf72HNi7B9lkYwoTQMOh71ussui1jHG+IBOoDxO21jl5cBJ20ZkX6n2AfAD4BtAIN4bFJG7RGSriGxtaWmJV3VckgVfvpQEtHT3879/9g5ul4Mf3XZJRkTcKcmzfK4V2PH73SfSPJKzJ5nfvGgrhpEfU7HqDFd5yn2IyA1AszHmnSj3T69szBpjzCJjzKLKyspE1ccNqkXjg+ZuD7f+59sc7/Twk9svYVpp5gY3KNGZVJzH+dOK2bB77K8rJSNKTcD0sNfVwNFYdUTEBZQA7XHaxipvBUptG5F9pdrHUuCTInIAa3pwhYj8Ion3q9gE56fVU8pedhzt5FMPvcWRjj7+63MfYVFN9h73kO2smFPFOwc76OgZSPdQzopkRGkLUGdHxbmxggrWRdRZB9xhX98ErDfWJ9o6YLUdOVcL1AGbY9m022ywbWDbfGYofRhjvmWMqTbG1Nj21xtjbkvyuShAcKO8alL2YYzhf7Ye5qYfb8QfMKy969JQ1mllbLJ8bhUBA6/vHdtLEAljJ40xPhG5G3gRcAKPGWN2iMh3gK3GmHXAo8DPRaQBy3tZbbfdISJPATsBH/AlY4wfIJpNu8tvAmtF5F7gPds2Q+lDOTsC6iJlJcc7PfzNb7axfnczi2vL+OGfXETVBI20G+tcWF1KeaGb3+9qZtXCyGX/sUNSAf3GmOeA5yLKvh127QFujtH2PuC+ZGza5Y3Y0XMR5Sn3EXb/VeDVWPeV6AQ1KRvCTBUr6/cjf2jkx6/tI2AM375hHn/60ZpxtdE0m3E4hGVzqnhl1wl8/sCYDVYZX7vMlJRQTyk78PoD/Oa9Izz48h6OdXq4bv4k/ubj5zGjvDDdQ1OGmavOq+JX7zbx7qGTLE5jIt2zQUVJiUlQlFSaxiYer5+nth7mJ681cuRkHwuqS/jBLQtZomtHWcvldRW4HBKamh2LqCgpMQkFOqgqjSlauvtZu/kQj288SOupfhbNmMi9nzqfZedWjtphfUp6KM7LYXFtGRt2N3PP9XPTPZwhoaKkKFmAMYZ3D53kZxsP8Ny2Y3j9ho+dW8lfLJvFktoyFaNxxIq5Vdz7u100dfRSPXHsHQmvoqQkgbpKmcqpfh/PfnCUn799kB1Hu5iQ6+K2S2dw+6UzmFlZlO7hKWkgKEobdjdz+2U16R5OyqgoKcoYIxAwbNrfzv+8c5jntx2nz+tnzqQJ3Pep87lx4TQKx1mWbOV0ZlYWUVNewO9VlJRsRdeUMoPD7b386t0mfvVuE4fb+5iQ6+LGi6Zy0yXTuficUp2iU0Isn1vFLzcdGpMH/42t0SppQTUpfZzq9/Hi9uM8/U4TGxvbEIGlsyr4+jVzuG7+ZPLdzsRGlHHHVXMn8V9vHuCthjaunje2Dv5TUVKUDMPj9fNqfQu//eAor+w6Qb8vwDllBXztmnP544unjcnFa2V0WVxbRqHbyfr6ZhUlJfvQ6buRxx8wbNzXxroPjvD89uN0e3yUF7pZ/ZHpfHLhVC4+Z6JOzylJ43Y5uKKukg27mzHGjKnfHRUlJSF68uzIYIzh/cMneeb9o/xu2zFauvspynVx3fzJrFo4lY/OKh+zqWKU9LNibhUv7DjOrmPdzJtanO7hJI2KkqKMIsYYdh7r4rltx/jtB8c41N6L2+VgxZwqVi2cyvK5VeTl6DqRcvZcOcc6G+6NhhYVJSW70Om7s8MYw/YjXfxu2zGe336Mg229OB3CR2eV85crZnPd+ZMpzstJ9zCVLGNScR6zq4p4o6GNuz42K93DSRoVJUUZAYwxfNDUyXPbjvHctmM0dfThcggfnV3BF6+cxbXzJ1NW6E73MJUs5/LZFazdcoh+n59c19jwwFWUlISoo5QcgYDhvcMneW7bMV7YfpwjJ/vIcQpLZ1fw5avquOa8SUxUIVJGkaWzK/jpWwd49+BJLps1NhLxqigpCdHzlGLjDxjeOdgREqLjXR7cTgdX1FXwtWvO5erzJlFSoFNzSnpYMrMMp0N4a1+ripIytlEhio3H6+fNhlZe2nGCV3adoK1nALfLwZXnVnLPBXNZcV6VrhEpGUFxXg4XVpfwRkMrX792TrqHkxQqSkpU/AEVpXA6e72srz/BSztO8NqeFnoH/BTlulg+t4pr5k1ixdwqijTnnJKBLJ1dwUMbGujyeMfElyX9K1Ki4vWrKB092cfLO0/w0s7jvN3Yjj9gqJqQy6cumsa18ydz6cyyMbN4rIxfls6u4D/WN/D2vjaunT853cNJiIqSEpUBXyB0PV5m8owx7Dlxipd2HOelnSfYdqQTgFmVhdz1sZlcO28SF1aX4nCMnd3xinLROaXk5zh5s6E1e0RJRFYC/wY4gUeMMd+NuJ8L/Ay4BGgDbjHGHLDvfQu4E/ADXzbGvBjPpojUAmuBMuBd4HZjzECqfYjIdLv+ZCAArDHG/FuqD2iss2F3MzMrC5lRXphSuwF/mChlcfzdgC/AlgPtrN/dzCu7TnCwrReAi88p5Z7r53LNvEnM0nOJlDFMrsvJ4toy3tzXlu6hJEVCURIRJ/AQcA3QBGwRkXXGmJ1h1e4EOowxs0VkNXA/cIuIzANWA/OBqcArInKu3SaWzfuBB40xa0XkYdv2j4fQhw/4ujHmXRGZALwjIi9HjDvr+dxPt+ByCA3//PGU2nnDRCnbONHl4dX6ZtbvbuaNva30DPhxOx1cNqucP//YLK4+r4qq4rx0D1NRho3LZpXz3ed303qqn4qi3HQPJy7JeEqLgQZjTCOAiKwFVgHhH+6rgP9rXz8N/FCsDICrgLXGmH5gv4g02PaIZlNEdgErgD+x6zxu2/1xqn0YYzYCxwCMMd227WkR485qgsLiG0LQQn8WTd/5A1aOuQ27m9lQ38yOo10ATCnJY9VF01g+p4qls8vH3LkzipIsi2vLANi8v52PXzAlzaOJTzJ/hdOAw2Gvm4AlseoYY3wi0gmU2+VvR7SdZl9Hs1kOnDTG+KLUH0ofAIhIDXARsCnaGxSRu4C7AM4555xoVcYkpzy+xJVi0NM/9LaZwMneAV7b08KG3c28tqeFjl4vDoFLZkzkGyvnsHxOFXMnTxhT2ZMVZahcMK2E/BwnmxrbskKUov3VRn53jlUnVnm01Mfx6g+lD6uRSBHwK+ArxpiuKHUxxqwB1gAsWrRojPsFg5w6C2EJbzsWPCVjDLuOdbOhvpkNu5t591AHAQNlhW6Wz6li2dwqrqyr1I2syrgkx+ngkhkT2bS/Pd1DSUgyotQETA97XQ0cjVGnSURcQAnQnqBttPJWoFREXLa3FF4/5T5EJAdLkH5pjPl1Eu81q+gd8A+5bbiXlama1NPv482GVluIWjje5QHg/GnF3L18NsvnVrGguhSnRsspCktqy3jglT2c7B2gtCBz010lI0pbgDo7Ku4IVlDBn0TUWQfcAWwEbgLWG2OMiKwD/ltEHsAKQqgDNmN5N2fYtNtssG2stW0+M5Q+7PWmR4FdxpgHUn0w2UDPwPB4SpnE/tae0NrQpsZ2BvwBinJdXFFXYXlEcyo1SEFRorC4tgxjYMuBDq7J4NNoE4qSvX5zN/AiVvj2Y8aYHSLyHWCrMWYd1of/z+0gg3YskcGu9xRWcIEP+JIxxg8Qzabd5TeBtSJyL/CebZtU+xCRy4HbgW0i8r5t42+MMc8N7VGNPYLrQkPxFLo83tB1OlMO9fv8bN5vhWy/Wt/C/tYeAGZXFXHHR2ewfG4Vi2aU4XbpYXiKEo8Lp5fidjnY1Ng2tkUJwP4gfy6i7Nth1x7g5hht7wPuS8amXd7IYIReeHlKfRhj3iD6etO4oaffmr5zD+H00hP2VFg6ONbZx4bdLWyob+bNhlZ6B/zkuqyQ7c8trWHZuVWcU16QtvEpylgkL8fJwumlbD6Q2etKGgObxQSn4HJzhiJK/aHrkfaTfP4A7x8+yfrdzWyob2HXMSseZVppPp++uJrlcyu5bGYF+W5N6aMoZ8OltWX8cEMD3R4vEzI0D56KUhbTesoSlpL81H/5jnT0UTUhl+bu/hFRpfaeAV7f08J6O2S7s8+L0yEsmjGRb10/l+Vzq6irKtKQbUUZRhbXlhNY38A7BztYNqcq3cOJiopSFtPabYlSTorTd8YYdh/v4vxpJZYoDQOWzW7W77YyKbxnh2xXFLm5Zt4kls+p4vK6iiEJqKIoyXHxjFJcDmHT/nYVJWX0OWavC/lSTBl0sK2Xjl4v508r4Q97W4ec+87j9bNxXxu/332C9buaOdppjWdBdQl/uaKOFXOruGBaiSY4VZRRosDtYkF1CZsaMzcPnopSFtNw4hSQ+jEUL+44DsCycyv58av7Umrb3jPASzuO8/LOE7y5rxWPN0CB28kVdRV85epzNWRbUdLM4tpyHvlDI30D/oxcp1VRylI6e700tARFKXlPqXfAx0/fOsDimjJqKqzM4oly57V09/PijuM8v/1Y6Nyh6on5rP7IOayYW8USPXdIUTKGJTPLePi1fbx7qIOlsyvSPZwzUFHKUn6x6SD+gOHC6hIOtfcm1cbnD/DNX23jeJeHH9yyMLQWFX62UhCvP8D63c08teUwG+qbCRiYWVHIF6+cxfUXTGbelGINUlCUDGTRjIk4BDY1tqkoKSOLMYY3G9p49I1GNtS3cP35k5lckkdjS0/CtjuPdnHv73by1r42vrFyDktmlodCysNFyeP18+SWwzz82j6OdXqonJDLn185i1ULpzJnkiY4VZRMZ0JeDvOnlvB2hubBU1HKAgZ8AZ798ChrXm9k9/FuKorcfOXqOr5w5SweeHnPaQf2RXK4vZd/fameZz44SnFeDv/vjy/g1sVWpnR3hKe0eX8733j6Aw609bJoxkT+8ZPzWTG3CtcQNucqipI+ltSW8bO3D+Lx+snLyaypdRWlMcypfh9rNx/ikT/s53iXh7qqIr736QWsumhqaA0nxylR14TaTvXzH+sb+OWmgzgdwhevnMWfXznrtJDsHKfl9Qz4Azzz/hG+/tQHTC3N52efX8wVdRXqFSnKGGXJzHIeeWM/Hxw+yZKZ5ekezmmoKI1BTnR5+MXbB3n8rQN0eXxcOrOM//fpC1h2buUZQuFyOPAHDIGAweEQegd8PPKH/ax5vZE+r5/PLJrOV66uY1KUiDgRwe1y8M7BDn7yWiMXz5jII3csojhDd4IripIci2vKEIFN+9tVlJSh4fMH2FDfwpNbDrF+dzMGuHbeJL5w5SwuOmdizHbBRKXeQIBXd7Xwf9ft4Finh5XzJ/N/rpvD7KqiuP3mOh28ta+Nygm5rLn9EhUkRckCSgpymDNpApv2t2EdrJA5qChlMF5/gLcb23hh+3Fe2nmClu5+Kifk8oUrZ/GZRdNDIdvxcNkbUx99Yz/fe6GeeVOK+Y9bL2JRTVlSY+i2gx3+YtmsjD6DRVGU1Lh0ZjlrtxxiwBfIqCz7KkoZRt+An9f2tPDSjuO8susEXR4fBW4ny+ZUsmrhNFbMrUopbVCw7vdeqOfKcytZ89lLUtozdPV5k9hQ3xwKflAUJTtYUlvGT986wLYjnVwyI/Zsy2ijopQBdPZ5Wb/7BC9sP85re1rweAOUFuRw7fzJXDd/MlfUVQw5QiYYrABw743np7yJ9QerF+LzBzIuQkdRlLNjca01W7Jpf5uKkmJlQXhp53Fe2H6cjfva8AUMk4pz+cyi6Vw3fzKLa8tSTqQajell1rlDn19aG7pOhaJc/RVRlGykvCiX2VVFbGps5y+WpXs0g+gnzihy9GQfL2y3hGjLwXaMgZryAu68opaV8ydzYXXpsCcnvfLcSn5x5xKWzExuDUlRlPHDktoynnn/KD5/IGP2G6oojTAer5/ntx/jqS1NbLQz886ZNIEvr6jj+gsmj3gWBBHh8rrMSyWiKEr6WTKznF9uOsTOY10sqC5N93AAFaURo6NngP966wCPv3WAzj4v55QV8PVrzuUTC6YwszJ+GLaiKMpocKm9rrRxX5uKUrbiDxh+8fZBvv9iPd39Pq6dN4nPLa1lSW2ZnhukKEpGUVWcx7wpxby88wR/fuWsdA8HgKQmEUVkpYjUi0iDiNwT5X6uiDxp398kIjVh975ll9eLyHWJbIpIrW1jr23TPdx9jBQnewe447HN/MO6HSw8p5QXv/Ix1nx2EZfNKldBUhQlI1l5/mTeOdRBc7cn3UMBkhAlEXECDwHXA/OAW0VkXkS1O4EOY8xs4EHgfrvtPGA1MB9YCfxIRJwJbN4PPGiMqQM6bNvD3cew09EzwM0Pb2Tz/na++8cX8LPPL2bO5Akj1Z2iKMqwcN38yRgDz35wLN1DAZLzlBYDDcaYRmPMALAWWBVRZxXwuH39NHCVWKv3q4C1xph+Y8x+oMG2F9Wm3WaFbQPb5o3D2UdyjyU1BnwBPvfTLRxs7+Wnn/8Iqxefo8lKFUUZE5w7qYiP1Ezk39fv5ejJvnQPJ6k1pWnA4bDXTcCSWHWMMT4R6QTK7fK3I9pOs6+j2SwHThpjfFHqD1cfZyAidwF32S9PiUgb0BqtbiKW3jeUVhlLBUN8DlmIPotB9FlYZN1zmPYPQ25aAcwYjjEkI0rRvvJHnoUQq06s8mgeWrz6w9nHmYXGrAHWBF+LyFZjzKJodccT+hwG0WcxiD4LC30Og9jPomY4bCUzfdcETA97XQ0cjVVHRFxACdAep22s8lag1LYR2ddw9aEoiqJkKMmI0hagzo6Kc2MFFayLqLMOuMO+vglYb4wxdvlqO3KuFitH+uZYNu02G2wb2DafGc4+knssiqIoSjpIOH1nr9/cDbwIOIHHjDE7ROQ7wFZjzDrgUeDnItKA5b2sttvuEJGngJ2AD/iSMcYPEM2m3eU3gbUici/wnm2bYe4jEWsSVxkX6HMYRJ/FIPosLPQ5DDJsz0IsZ0NRFEVR0k9mZOBTFEVRFFSUFEVRlAxCRSmM0U5LlA5E5DERaRaR7WFlZSLysp3a6WURmWiXi4j8u/08PhSRi8PaeI7xQgAAA3BJREFU3GHX3ysid0TrK5MRkekiskFEdonIDhH5K7t8PD6LPBHZLCIf2M/iH+3yYUv5NZawM8K8JyLP2q/H63M4ICLbROR9Edlql43834cxRn+sdTUnsA+YCbiBD4B56R7XCLzPjwEXA9vDyr4H3GNf3wPcb19/HHgeay/YpcAmu7wMaLT/nWhfT0z3e0vxOUwBLravJwB7sNJRjcdnIUCRfZ0DbLLf41PAarv8YeCL9vVfAA/b16uBJ+3refbfTS5Qa/89OdP9/obwPL4G/DfwrP16vD6HA0BFRNmI/32opzTIqKUlSifGmNexohfDCU/hFJna6WfG4m2sPWRTgOuAl40x7caYDuBlrLyDYwZjzDFjzLv2dTewCysTyHh8FsYYc8p+mWP/GIYv5deYQUSqgU8Aj9ivhzP1WTYw4n8fKkqDREunNC1G3WxjkjHmGFgf1kCVXR7rmWTVs7KnXS7C8hDG5bOwp6zeB5qxPjj2kWTKLyA85ddYfxY/AL4BBOzXSac+I7ueA1hfTF4SkXfESsUGo/D3oecpDZJMOqXxRqqpncYcIlIE/Ar4ijGmS2In0s3qZ2GsvX0LRaQU+A1wXrRq9r9Z+SxE5Aag2RjzjogsCxZHqZrVzyGMpcaYoyJSBbwsIrvj1B22Z6Ge0iDjOS3RCdvVxv632S7P6hROIpKDJUi/NMb82i4el88iiDHmJPAq1rrAcKX8GissBT4pIgewpu9XYHlO4+05AGCMOWr/24z1RWUxo/D3oaI0yHhOSxSewikytdNn7ciaS4FO22V/EbhWRCba0TfX2mVjBnvu/1FglzHmgbBb4/FZVNoeEiKSD1yNtcY2XCm/xgTGmG8ZY6qNlVh0Ndb7+l+Ms+cAICKFIjIheI31e72d0fj7SHeERyb9YEWQ7MGaT//bdI9nhN7jE8AxwIv1LeZOrHnw3wN77X/L7LqCdVDiPmAbsCjMzuexFnAbgM+l+30N4TlcjjWN8CHwvv3z8XH6LBZgpfT60P7g+bZdPhPrw7QB+B8g1y7Ps1832Pdnhtn6W/sZ1QPXp/u9ncUzWcZg9N24ew72e/7A/tkR/Dwcjb8PTTOkKIqiZAw6facoiqJkDCpKiqIoSsagoqQoiqJkDCpKiqIoSsagoqQoiqJkDCpKiqIoSsagoqQoiqJkDP8/oHrNcmlroVQAAAAASUVORK5CYII=\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": 14,
"metadata": {},
"outputs": [],
"source": [
"\n",
"\n",
"# probs = mixture.prob(test_q)\n",
"# probs_np = zfit.run(probs)\n",
"# probs_np *= np.max(calcs_test) / np.max(probs_np)\n",
"# plt.figure()\n",
"# plt.semilogy(test_q, probs_np,label=\"importance sampling\")\n",
"# plt.semilogy(test_q, calcs_test, label = 'pdf')\n"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"# 0.213/(0.00133+0.213+0.015)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Adjust scaling of different parts"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"total_f.update_integration_options(draws_per_dim=2000000, mc_sampler=None)\n",
"# inte = total_f.integrate(limits = (950., 1050.), 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": 17,
"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][0]*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": 18,
"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": 19,
"metadata": {},
"outputs": [],
"source": [
"total_f._sample_and_weights = UniformSampleAndWeights"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [],
"source": [
"# 0.00133/(0.00133+0.213+0.015)*(x_max-3750)/(x_max-x_min)"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"# zfit.settings.set_verbosity(10)"
]
},
{
"cell_type": "code",
"execution_count": 22,
"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": 23,
"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": 24,
"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": 25,
"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": 26,
"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": 27,
"metadata": {},
"outputs": [],
"source": [
"# jpsi_width"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [],
"source": [
"# plt.hist(sample, weights=1 / prob(sample))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Fitting"
]
},
{
"cell_type": "code",
"execution_count": 29,
"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": 30,
"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": 31,
"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": 32,
"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": 33,
"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": 34,
"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": 35,
"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": 36,
"metadata": {},
"outputs": [],
"source": [
"# 0.15**2*4.2/1000\n",
"# result.hesse()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Constraints"
]
},
{
"cell_type": "code",
"execution_count": 37,
"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), ztf.constant(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",
"# 7. Constraint on Ctt with higher limits\n",
"\n",
"# Ctt_abs = tf.pow(tf.pow(Ctt, 2.), 0.5)\n",
"\n",
"# constraint7 = tf.cond(tf.greater_equal(Ctt_abs, 0.5), lambda: 100000., lambda: 0.)\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, #constraint4_0, constraint4_1, constraint4_2,\n",
" constraint6_0, constraint6_1]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Reset params"
]
},
{
"cell_type": "code",
"execution_count": 38,
"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": 39,
"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 = 3.516E+05 | Ncalls=867 (867 total) |\n",
"| EDM = 6.79E-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 | False | False | True |\n",
"------------------------------------------------------------------\n",
"Function minimum: 351616.34467140574\n",
"----------------------------------------------------------------------------------------------\n",
"| | Name | Value | Hesse Err | Minos Err- | Minos Err+ | Limit- | Limit+ | Fixed |\n",
"----------------------------------------------------------------------------------------------\n",
"| 0 | DDstar_p | 1.95 | 0.29 | | |-6.28319 | 6.28319 | |\n",
"| 1 | p3770_s | 3.27 | 0.21 | | |0.918861 | 4.08114 | |\n",
"| 2 | bplus_0 | 0.479 | 0.018 | | | -2 | 2 | |\n",
"| 3 | Ctt | -0.44 | 0.19 | | | -1 | 1 | |\n",
"| 4 | bplus_2 | -0.23 | 0.08 | | | -2 | 2 | |\n",
"| 5 | Dbar_p | 5.30 | 0.26 | | |-6.28319 | 6.28319 | |\n",
"| 6 | p4040_p | 3.79 | 0.17 | | |-6.28319 | 6.28319 | |\n",
"| 7 | psi2s_p | 1.903 | 0.028 | | |-6.28319 | 6.28319 | |\n",
"| 8 | bplus_1 | -0.89 | 0.04 | | | -2 | 2 | |\n",
"| 9 | p4415_s | 1.09 | 0.18 | | |0.126447 | 2.35355 | |\n",
"| 10| p3770_p | -2.60 | 0.09 | | |-6.28319 | 6.28319 | |\n",
"| 11| DDstar_s | -0.300 | 0.016 | | | -0.3 | 0.3 | |\n",
"| 12| p4040_s | 1.02 | 0.16 | | |0.00501244| 2.01499 | |\n",
"| 13| p4160_p | -2.08 | 0.10 | | |-6.28319 | 6.28319 | |\n",
"| 14| p4415_p | 4.22 | 0.18 | | |-6.28319 | 6.28319 | |\n",
"| 15| Dbar_s | 0.300 | 0.013 | | | -0.3 | 0.3 | |\n",
"| 16| jpsi_p | 4.640 | 0.023 | | |-6.28319 | 6.28319 | |\n",
"| 17| p4160_s | 2.15 | 0.16 | | | 0.71676 | 3.68324 | |\n",
"----------------------------------------------------------------------------------------------\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"| | DDstar_p p3770_s bplus_0 Ctt bplus_2 Dbar_p p4040_p psi2s_p bplus_1 p4415_s p3770_p DDstar_s p4040_s p4160_p p4415_p Dbar_s jpsi_p p4160_s |\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"| DDstar_p | 1.000 0.173 0.000 -0.171 -0.324 -0.080 0.106 0.003 0.389 -0.061 0.262 0.029 -0.140 0.217 -0.024 0.003 0.173 -0.100 |\n",
"| p3770_s | 0.173 1.000 0.045 -0.234 -0.148 0.058 -0.027 -0.423 0.095 0.000 -0.171 0.024 0.080 0.020 -0.021 0.025 -0.006 -0.011 |\n",
"| bplus_0 | 0.000 0.045 1.000 -0.008 -0.011 0.019 0.022 -0.007 -0.832 0.017 0.025 0.000 0.018 0.014 0.018 0.001 -0.064 0.035 |\n",
"| Ctt | -0.171 -0.234 -0.008 1.000 0.689 -0.326 -0.291 0.166 -0.184 0.221 -0.263 -0.004 0.368 -0.425 -0.073 0.009 0.130 0.258 |\n",
"| bplus_2 | -0.324 -0.148 -0.011 0.689 1.000 -0.134 -0.069 -0.013 -0.337 -0.054 -0.134 0.005 0.099 -0.085 0.177 0.004 0.052 0.123 |\n",
"| Dbar_p | -0.080 0.058 0.019 -0.326 -0.134 1.000 0.011 0.052 0.180 -0.008 0.366 0.002 -0.089 0.105 -0.044 0.015 0.302 -0.091 |\n",
"| p4040_p | 0.106 -0.027 0.022 -0.291 -0.069 0.011 1.000 -0.228 0.020 0.031 0.180 0.029 -0.241 0.163 0.099 0.022 -0.071 0.295 |\n",
"| psi2s_p | 0.003 -0.423 -0.007 0.166 -0.013 0.052 -0.228 1.000 0.051 0.010 0.058 0.024 0.009 -0.131 -0.105 0.024 0.004 -0.083 |\n",
"| bplus_1 | 0.389 0.095 -0.832 -0.184 -0.337 0.180 0.020 0.051 1.000 0.100 0.128 -0.005 0.010 0.019 -0.100 -0.005 0.105 0.001 |\n",
"| p4415_s | -0.061 0.000 0.017 0.221 -0.054 -0.008 0.031 0.010 0.100 1.000 -0.081 -0.000 0.154 -0.055 -0.131 -0.001 -0.039 0.309 |\n",
"| p3770_p | 0.262 -0.171 0.025 -0.263 -0.134 0.366 0.180 0.058 0.128 -0.081 1.000 0.019 -0.177 0.252 0.072 0.022 0.115 -0.082 |\n",
"| DDstar_s | 0.029 0.024 0.000 -0.004 0.005 0.002 0.029 0.024 -0.005 -0.000 0.019 1.000 0.003 0.038 0.026 -0.001 0.054 0.007 |\n",
"| p4040_s | -0.140 0.080 0.018 0.368 0.099 -0.089 -0.241 0.009 0.010 0.154 -0.177 0.003 1.000 -0.562 -0.246 -0.001 -0.036 0.024 |\n",
"| p4160_p | 0.217 0.020 0.014 -0.425 -0.085 0.105 0.163 -0.131 0.019 -0.055 0.252 0.038 -0.562 1.000 0.282 0.024 0.016 -0.187 |\n",
"| p4415_p | -0.024 -0.021 0.018 -0.073 0.177 -0.044 0.099 -0.105 -0.100 -0.131 0.072 0.026 -0.246 0.282 1.000 0.014 -0.019 -0.216 |\n",
"| Dbar_s | 0.003 0.025 0.001 0.009 0.004 0.015 0.022 0.024 -0.005 -0.001 0.022 -0.001 -0.001 0.024 0.014 1.000 0.040 0.004 |\n",
"| jpsi_p | 0.173 -0.006 -0.064 0.130 0.052 0.302 -0.071 0.004 0.105 -0.039 0.115 0.054 -0.036 0.016 -0.019 0.040 1.000 -0.068 |\n",
"| p4160_s | -0.100 -0.011 0.035 0.258 0.123 -0.091 0.295 -0.083 0.001 0.309 -0.082 0.007 0.024 -0.187 -0.216 0.004 -0.068 1.000 |\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"Hesse errors: OrderedDict([(<zfit.Parameter 'DDstar_p' floating=True>, {'error': 0.2850165410517804}), (<zfit.Parameter 'p3770_s' floating=True>, {'error': 0.2124961626399453}), (<zfit.Parameter 'bplus_0' floating=True>, {'error': 0.01829244045185119}), (<zfit.Parameter 'Ctt' floating=True>, {'error': 0.19189240704670774}), (<zfit.Parameter 'bplus_2' floating=True>, {'error': 0.07723151076900314}), (<zfit.Parameter 'Dbar_p' floating=True>, {'error': 0.26113738425875166}), (<zfit.Parameter 'p4040_p' floating=True>, {'error': 0.16593345938513693}), (<zfit.Parameter 'psi2s_p' floating=True>, {'error': 0.02822480404964267}), (<zfit.Parameter 'bplus_1' floating=True>, {'error': 0.037964767888795214}), (<zfit.Parameter 'p4415_s' floating=True>, {'error': 0.17890980886132019}), (<zfit.Parameter 'p3770_p' floating=True>, {'error': 0.09336821683842733}), (<zfit.Parameter 'DDstar_s' floating=True>, {'error': 0.016393727908621925}), (<zfit.Parameter 'p4040_s' floating=True>, {'error': 0.1615214697208751}), (<zfit.Parameter 'p4160_p' floating=True>, {'error': 0.0976483880480612}), (<zfit.Parameter 'p4415_p' floating=True>, {'error': 0.17860121118891303}), (<zfit.Parameter 'Dbar_s' floating=True>, {'error': 0.01277429819793291}), (<zfit.Parameter 'jpsi_p' floating=True>, {'error': 0.022652863795177502}), (<zfit.Parameter 'p4160_s' floating=True>, {'error': 0.15625965574324152})])\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\ipykernel_launcher.py:166: UserWarning: Creating legend with loc=\"best\" can be slow with large amounts of data.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Toy 1/2\n",
"Time taken: 4 min, 18 s\n",
"Projected time left: 4 min, 18 s\n",
"Toy 1: Generating data...\n",
"Toy 1: Data generation finished\n",
"Toy 1: Loading data...\n",
"Toy 1: Loading data finished\n",
"Toy 1: Fitting pdf...\n",
"------------------------------------------------------------------\n",
"| FCN = 7.032E+05 | Ncalls=914 (914 total) |\n",
"| EDM = 0.000618 (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: 703225.2607413743\n",
"----------------------------------------------------------------------------------------------\n",
"| | Name | Value | Hesse Err | Minos Err- | Minos Err+ | Limit- | Limit+ | Fixed |\n",
"----------------------------------------------------------------------------------------------\n",
"| 0 | DDstar_p | 1.96 | 0.22 | | |-6.28319 | 6.28319 | |\n",
"| 1 | p3770_s | 3.08 | 0.16 | | |0.918861 | 4.08114 | |\n",
"| 2 | bplus_0 | 0.475 | 0.013 | | | -2 | 2 | |\n",
"| 3 | Ctt | -0.39 | 0.14 | | | -1 | 1 | |\n",
"| 4 | bplus_2 | -0.24 | 0.05 | | | -2 | 2 | |\n",
"| 5 | Dbar_p | -4.08 | 0.22 | | |-6.28319 | 6.28319 | |\n",
"| 6 | p4040_p | 3.71 | 0.12 | | |-6.28319 | 6.28319 | |\n",
"| 7 | psi2s_p | 1.961 | 0.025 | | |-6.28319 | 6.28319 | |\n",
"| 8 | bplus_1 | -0.875 | 0.027 | | | -2 | 2 | |\n",
"| 9 | p4415_s | 1.09 | 0.13 | | |0.126447 | 2.35355 | |\n",
"| 10| p3770_p | 3.69 | 0.07 | | |-6.28319 | 6.28319 | |\n",
"| 11| DDstar_s | -0.300 | 0.011 | | | -0.3 | 0.3 | |\n",
"| 12| p4040_s | 1.02 | 0.12 | | |0.00501244| 2.01499 | |\n",
"| 13| p4160_p | -2.10 | 0.07 | | |-6.28319 | 6.28319 | |\n",
"| 14| p4415_p | 4.18 | 0.13 | | |-6.28319 | 6.28319 | |\n",
"| 15| Dbar_s | -0.300 | 0.008 | | | -0.3 | 0.3 | |\n",
"| 16| jpsi_p | -1.642 | 0.017 | | |-6.28319 | 6.28319 | |\n",
"| 17| p4160_s | 2.12 | 0.11 | | | 0.71676 | 3.68324 | |\n",
"----------------------------------------------------------------------------------------------\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"| | DDstar_p p3770_s bplus_0 Ctt bplus_2 Dbar_p p4040_p psi2s_p bplus_1 p4415_s p3770_p DDstar_s p4040_s p4160_p p4415_p Dbar_s jpsi_p p4160_s |\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"| DDstar_p | 1.000 0.165 -0.005 -0.159 -0.333 -0.114 0.101 -0.005 0.405 -0.058 0.248 0.040 -0.141 0.222 -0.026 0.004 0.172 -0.100 |\n",
"| p3770_s | 0.165 1.000 0.043 -0.254 -0.143 0.044 0.025 -0.489 0.087 0.000 -0.162 0.023 0.074 0.047 -0.001 0.024 0.005 0.006 |\n",
"| bplus_0 | -0.005 0.043 1.000 -0.011 -0.013 0.020 0.025 -0.011 -0.821 0.016 0.023 0.000 0.015 0.017 0.020 0.001 -0.062 0.035 |\n",
"| Ctt | -0.159 -0.254 -0.011 1.000 0.685 -0.354 -0.282 0.181 -0.196 0.215 -0.292 -0.004 0.377 -0.430 -0.067 0.012 0.077 0.257 |\n",
"| bplus_2 | -0.333 -0.143 -0.013 0.685 1.000 -0.142 -0.063 -0.025 -0.347 -0.058 -0.155 0.004 0.105 -0.093 0.179 0.004 0.025 0.130 |\n",
"| Dbar_p | -0.114 0.044 0.020 -0.354 -0.142 1.000 0.005 0.093 0.195 -0.001 0.407 0.002 -0.085 0.119 -0.048 0.022 0.369 -0.099 |\n",
"| p4040_p | 0.101 0.025 0.025 -0.282 -0.063 0.005 1.000 -0.275 0.020 0.039 0.168 0.034 -0.245 0.148 0.095 0.026 -0.063 0.304 |\n",
"| psi2s_p | -0.005 -0.489 -0.011 0.181 -0.025 0.093 -0.275 1.000 0.070 0.015 0.075 0.038 0.031 -0.155 -0.129 0.038 0.015 -0.104 |\n",
"| bplus_1 | 0.405 0.087 -0.821 -0.196 -0.347 0.195 0.020 0.070 1.000 0.100 0.156 -0.004 0.001 0.034 -0.101 -0.005 0.138 -0.010 |\n",
"| p4415_s | -0.058 0.000 0.016 0.215 -0.058 -0.001 0.039 0.015 0.100 1.000 -0.078 -0.001 0.151 -0.055 -0.135 -0.001 -0.037 0.312 |\n",
"| p3770_p | 0.248 -0.162 0.023 -0.292 -0.155 0.407 0.168 0.075 0.156 -0.078 1.000 0.025 -0.185 0.260 0.061 0.029 0.164 -0.093 |\n",
"| DDstar_s | 0.040 0.023 0.000 -0.004 0.004 0.002 0.034 0.038 -0.004 -0.001 0.025 1.000 0.002 0.049 0.032 -0.002 0.069 0.008 |\n",
"| p4040_s | -0.141 0.074 0.015 0.377 0.105 -0.085 -0.245 0.031 0.001 0.151 -0.185 0.002 1.000 -0.559 -0.243 -0.003 -0.039 0.007 |\n",
"| p4160_p | 0.222 0.047 0.017 -0.430 -0.093 0.119 0.148 -0.155 0.034 -0.055 0.260 0.049 -0.559 1.000 0.277 0.030 0.040 -0.183 |\n",
"| p4415_p | -0.026 -0.001 0.020 -0.067 0.179 -0.048 0.095 -0.129 -0.101 -0.135 0.061 0.032 -0.243 0.277 1.000 0.017 -0.023 -0.205 |\n",
"| Dbar_s | 0.004 0.024 0.001 0.012 0.004 0.022 0.026 0.038 -0.005 -0.001 0.029 -0.002 -0.003 0.030 0.017 1.000 0.051 0.003 |\n",
"| jpsi_p | 0.172 0.005 -0.062 0.077 0.025 0.369 -0.063 0.015 0.138 -0.037 0.164 0.069 -0.039 0.040 -0.023 0.051 1.000 -0.078 |\n",
"| p4160_s | -0.100 0.006 0.035 0.257 0.130 -0.099 0.304 -0.104 -0.010 0.312 -0.093 0.008 0.007 -0.183 -0.205 0.003 -0.078 1.000 |\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"Hesse errors: OrderedDict([(<zfit.Parameter 'DDstar_p' floating=True>, {'error': 0.21976002920084792}), (<zfit.Parameter 'p3770_s' floating=True>, {'error': 0.15829775347665453}), (<zfit.Parameter 'bplus_0' floating=True>, {'error': 0.013027835592936077}), (<zfit.Parameter 'Ctt' floating=True>, {'error': 0.14012275276008618}), (<zfit.Parameter 'bplus_2' floating=True>, {'error': 0.0548143444334922}), (<zfit.Parameter 'Dbar_p' floating=True>, {'error': 0.21679058977592525}), (<zfit.Parameter 'p4040_p' floating=True>, {'error': 0.11774283300825683}), (<zfit.Parameter 'psi2s_p' floating=True>, {'error': 0.02483027542093197}), (<zfit.Parameter 'bplus_1' floating=True>, {'error': 0.027345870305282904}), (<zfit.Parameter 'p4415_s' floating=True>, {'error': 0.12735002863519618}), (<zfit.Parameter 'p3770_p' floating=True>, {'error': 0.07047893274838923}), (<zfit.Parameter 'DDstar_s' floating=True>, {'error': 0.010590941888212607}), (<zfit.Parameter 'p4040_s' floating=True>, {'error': 0.1157783313417895}), (<zfit.Parameter 'p4160_p' floating=True>, {'error': 0.07019659457672}), (<zfit.Parameter 'p4415_p' floating=True>, {'error': 0.12646052173735356}), (<zfit.Parameter 'Dbar_s' floating=True>, {'error': 0.008140502524761783}), (<zfit.Parameter 'jpsi_p' floating=True>, {'error': 0.016511357460833764}), (<zfit.Parameter 'p4160_s' floating=True>, {'error': 0.11145758454315047})])\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Toy 2/2\n",
"Time taken: 9 min, 6 s\n",
"Projected time left: \n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAaUAAAD4CAYAAABMtfkzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO29eXzc1XX3/z4zo9XWbsmrbAks7xgbjIGYBAcINiFgkkLjJBCehD60CWmztE3glzbpQ0Of0KaQtAkQ/wItJYuhhBSHGBwCGEICxgZveJd3WV60WdauWe7zx/fOaDye0czIGs1IOu/Xyy/N3Lnfc+98Lc1nzrnnnivGGBRFURQlE3ClewKKoiiKEkRFSVEURckYVJQURVGUjEFFSVEURckYVJQURVGUjMGT7glkGuPGjTNVVVXpnoaiZDzvH2tl3NgcJhTlntXe7Q2w71QbU0vzKcrLStPsHMLn6PUH2H2ijcnFeZSOyU7rvEYi7777bqMxpvx87agoRVBVVcWmTZvSPQ1FyXhmfPNFPn9VNffeMOus9tpT7Vz30Ov866cWcvPFk9I0O4fwOZ48083l//QK//jxeXzm8mlpnddIREQOD4YdDd8pijIg/Mbgcck57W7b5g8EhnpK52AwBKfoEudBIKB7MzMZFSVFUZLGGIM/YHBFESVPSJSGelbnEjBgtShMLFWUMhkVJUVRkib4ue6Wc0XJlUmekjEIznyCc/WrJmU0uqakKErSBL0Nd5SvtZnkKRkIhe/c7rPDdy+9f5z//OMhfv5nV0T1+AYLr9dLXV0d3d3dKRtjKMnNzWXKlClkZaUmiUVFSVGUpAnYmpnRPsyDazeZ4SkRit8FPSWfFaUvr95Cjy/A6S5vSrPx6urqKCgooKqqConiWQ4njDE0NTVRV1dHdXV1SsbQ8J2iKEkT8pSifMh6MmTtJlhsOjhDl/20Cwpqtsdp6PL6UzqP7u5uysrKhr0gAYgIZWVlKfX6VJQURUkavwmG72KvKfnSLkrOT1eEpxQUy2wbe+zq9aV8LiNBkIKk+r2oKCmKkjTBdRlXP55SIM3H4gTHj5V9F3ze7U1/mFHpQ0VJUZSkCXpBHnfsfUpp95Tsz6AzJyKInCuW6Z5nprB+/Xo+9rGPAdDT08N1113HggULePrpp4d0HprooChK0vTnKYU8kjTnXvd5Sn1z9Lgk5CkF5+7LhDTBDGPz5s14vV62bNky5GMn5CmJyHIR2SMitSJyb5TXc0Tkafv6BhGpCnvtPtu+R0SWxbMpItXWxj5rM7u/MUSkTEReE5F2EflhjPmvEZH3E7sliqLEo781pb79QJmxphSOS/pEKahV3lGwcenQoUPMmjWLO++8k/nz53PrrbfS2dnJSy+9xKxZs7jqqqt47rnnADh16hS33347W7ZsYcGCBezfv39I5xrXUxIRN/Aj4CNAHbBRRNYYY3aGdbsLaDHGTBeRlcCDwCdFZA6wEpgLTAJ+JyIz7DWxbD4IPGyMWS0ij1nbj8YaA+gG/h6YZ/9Fzv8TQHtSd0VRlH7pL/vO5XLCZOnOvgsS7s25wzylIL4hTF3/P7/ewc76M4Nqc86kQr5909y4/fbs2cPjjz/OkiVL+PznP89DDz3Ej3/8Y1599VWmT5/OJz/5SQAqKir4yU9+wve+9z1eeOGFQZ1rIiTiKS0Gao0xB4wxvcBqYEVEnxXAk/bxs8C14vjMK4DVxpgeY8xBoNbai2rTXnONtYG1eUt/YxhjOowxb+KI01mIyFjga8B3EnifiqIkSPBzPNamU0+UD/+hJjLRARwR9UekivtGgacEUFlZyZIlSwC4/fbb2bRpE9XV1dTU1CAi3H777WmeoUMia0qTgaNhz+uAy2P1Mcb4RKQVKLPtb0dcO9k+jmazDDhtjPFF6R9rjMZ+5v6PwL8Cnf29QRG5G7gbYOrUqf11VRSF8PBd9NfDw2TpIhi+C5dNl0tC62HBtSbvEK4pJeLRpIrIVO7W1taMTFVPxFOKNuvI37ZYfQarPdF59E1IZAEw3Rjzq1h9QkaMWWWMWWSMWVReft7HgSjKiCcyWSCSTPCU+rLvIsJ3EYtNo2FNCeDIkSO89dZbAPziF7/guuuu4+DBg6E1o1/84hfpnF6IRESpDqgMez4FqI/VR0Q8QBHQ3M+1sdobgWJrI3KsWGPE4krgUhE5BLwJzBCR9f2+U0VREiIYGvO4on+EuFyS9lTrqOE7l5xTk28o15TSyezZs3nyySeZP38+zc3NfPWrX2XVqlXceOONXHXVVUyblhlnTCUSvtsI1IhINXAMJ3Hh0xF91gB3Am8BtwKvGmOMiKwBfi4iD+EkOtQA7+B4PefYtNe8Zm2stjaf72+MWJM2xjyKkyCBzdR7wRizNIH3qyhKHILrMLHCdx6XpH3zbCh8F+4piYRq8o2m7DsAl8vFY489dlbb8uXL2b179zl9ly5dytKlS4doZmcTV5Ts+s2XgHWAG3jCGLNDRO4HNhlj1gCPA0+JSC2O97LSXrtDRJ4BdgI+4B5jjB8gmk075DeA1SLyHWCztU2sMaytQ0AhkC0itwDXR2QHKooyiIQKssYI37kzwFOKrH0HZ3tKwanrPqXMIqHNs8aYtcDaiLZvhT3uBm6Lce0DwAOJ2LTtB3Cy8yLb+xujKs78DxElXVxRlIERWaonEndYQkG66POU+tpcrnMrOngzJHU9lVRVVfH++8Njq6aWGVIUJWn8/RxdAU6YLO2ekv3pOid8F/Sghq6iQz8rDcOOVL8XFSVFUZIm0M/mWXAO1Et39l20RAdXWPZdX/gutfPMzc2lqalpRAhT8Dyl3NzclI2hte8URUmaeOE7j8uVdlGKtk/JEyWs6E1x9t2UKVOoq6ujoaEhpeMMFcGTZ1OFipKiKEnTX+07cCpzp12UOLcgqyssrDhUFR2ysrJSdkrrSETDd4qiJM2w8pQi9ikFPaXg7DT7LrNQUVIUJWniVXTIhM2zkSfPwtkVHYIe0mjIvhtOqCgpipI0gTjhu0zYPBuIsk8pvCZf8Kd6SpmFipKiKEkT/ByPlX2XEZ6S/XlO+C7oKdn5jZaKDsMFFSVFUZImFL7rr8xQukUpysmz4ecpBWvejZbad8MFFSVFUZImXvjO2Tyb3g/7aCnh4Ztng8e1pzshQzkbFSVFUZKmv5NnIZjlNpQzOpdoBVmzPK5QuE7Dd5mJipKiKEmTSO27tHtKBDME+9qyXBI61C8UvtNEh4xCRUlRlKRJRJTSHRYLhDylvrYstyuUCh7ylDR8l1GoKCmKkjT+OEdXZLkl7WExE2WOHrfjKQUCJhTeU08ps1BRUhQlaQJxPKUstyvt4btoDlCW24U3EDgrXT3VZYaU5FBRUhQlaeLVvvO4XWn3lIhS+y7LLfj85izB1PBdZqGipChK0gTilBnKcvclFKSLvjJDfW2OWEZ6Shq+yyRUlBRFSZp4iQ5ZLlfaRSmU6BC2UynbenD+MC8u3ZUnlLNRUVIUJWmCn+kxRckjaV+r6Tu6oq/NY1PCw89QUk8ps1BRUhQlafz2Q72/oyt6MzR85/Obs9LV1VPKLBISJRFZLiJ7RKRWRO6N8nqOiDxtX98gIlVhr91n2/eIyLJ4NkWk2trYZ21m9zeGiJSJyGsi0i4iPwyzky8ivxGR3SKyQ0S+m/ztURQlGsEkBk/M7Lv0e0p9VcrDw3dCrz9w1tzSn5ChhBNXlETEDfwIuAGYA3xKROZEdLsLaDHGTAceBh60184BVgJzgeXAIyLijmPzQeBhY0wN0GJtxxwD6Ab+HvibKNP/njFmFrAQWCIiN8R7v4qixCfoacQWpfSvKcXylAB6fBq+y1QS8ZQWA7XGmAPGmF5gNbAios8K4En7+FngWnHyMFcAq40xPcaYg0CttRfVpr3mGmsDa/OW/sYwxnQYY97EEacQxphOY8xr9nEv8B6QuoPlFWUUEfwg7y8l3BcwoQ2s6SBq7TsrSt1ef6hNw3eZRSKiNBk4Gva8zrZF7WOM8QGtQFk/18ZqLwNOWxuRY8UaIy4iUgzcBLwS4/W7RWSTiGxqaGhIxKSijGq8AUOWW876wA8n2+20pzM0Fkp0CGvLsvPqsqKUl+VOu0ennE0iohTtty7yNy1Wn8FqT3Qe5yAiHuAXwL8ZYw5E62OMWWWMWWSMWVReXh7PpKKMevwBgyfWYUr0hcnSWdUhFL4Lm2bQU+rqtaKU7U772pdyNomIUh1QGfZ8ClAfq48VgSKguZ9rY7U3AsXWRuRYscaIxypgnzHm+wn0VRQlAbz+QMz1JOj78E+np9R3HPrZte+gz1PK9aS/HJJyNomI0kagxmbFZeMkLqyJ6LMGuNM+vhV41TjB5DXASps5Vw3UAO/Esmmvec3awNp8Ps4YMRGR7+CI11cSeJ+KoiSIz29CH/DRyAqF79L3gR+1Srjr7DWl3Cy3rillGJ54HYwxPhH5ErAOcANPGGN2iMj9wCZjzBrgceApEanF8V5W2mt3iMgzwE7AB9xjjPEDRLNph/wGsNoKymZrm1hjWFuHgEIgW0RuAa4HzgDfBHYD79nY9w+NMT9J/jYpihKOL2BCIbpoBD2ldIbGTJT6fFke5/FZoqThu4wirigBGGPWAmsj2r4V9rgbuC3GtQ8ADyRi07YfwMnOi2zvb4yqGFOP/VVOUZQB44sTvgu+lk5PyR+lPl9wHSy4ppSblf7UdeVstKKDoihJ43hKsUUp2xNcU0p/+M4VLSXc7lPS8F3moaKkxGTX8TN881fbQxWhFSWI1x8Irc9EI+iRZEKiw1nHoQcTHXr7wnf+NO+nUs5GRUmJyef/cyM/23CEE2e643dWRhX+gIm5cRb6stzS6ymdu6bkcUcmOqRfPJWzUVFSFCVpvP7+Ex2y3ZkTvos85A/CU8LdQHr3Uylno6KkKErS+AKB0Ad8NIKeUjrXa/oOIuxri9w8m5vtiJJ6SpmDipKiKEkTL3wX2jzry6zwXdCD64z0lDQDL2NQUVLiot8hlUjiJTqENs+m0VOKlhKeY9eQOnqc8prBNSW/JvNkDCpKiqIkTfyKDpngKTk/w0Up6Bm1dzuilJdlw3cqShmDipISF92BrEQSr6JDMCU8nQkEoZTwsGnmWhFqD3lKGr7LNFSUlLjod0glEl+g/4oO2Z70H10RWlMK95SC4bteK0qa6JBxqCgpipI0Pr+JU2YoM1PCg55RR09flXDQlPBMQkVJUZSk8QVMaN0oGlmZUGYoSkp4jp3XueE79ZQyBRUlRVGSxucP9JsSHky97s2wlHARIcfjCs0rKEpalDVzUFFSFCVpvHGy74Kp1z1pFKVoKeHQJ0QifWtMQ7HJ94+1jdzzs/fSKtTDARUlJSaadafEwhfof59SMEyWTlHqOw797N/k4NxyPK4hXft6YO0ufrP9OJuPtKR8rOGMipISE42yK7HwBwzu/o6ucLsQgR5bOSEd+KNUCYc+TynH4w5t8h2KzbP7TrUDsL+hI+VjDWdUlBRFSRqv35DVz5pScO0mnZ5StJRw6AvZ5XhcfZt8h8BTyrfp5ydau1I+1nBGRUmJiYbvlFj4/IF+N8+C44mkV5ScnxJjTSkny9W39uVN7TyNMaHSRg3tPSkda7ijoqQoStJ445w8C1hPKX3hu2gp4dBXaijH4w497k7xPLu9gdAG3Ya23pSONdxRUVIUJSmMMfT6AuTE85SyXCn3QPojWko4wNhcD+CIZtBr6k7xPM90e0OP1VPqn4RESUSWi8geEakVkXujvJ4jIk/b1zeISFXYa/fZ9j0isiyeTRGptjb2WZvZ/Y0hImUi8pqItIvIDyPmdamIbLfX/JtE+vGKoiRN8Bt/tiezw3fB5IXIP/uivCzAEaVgJl53ihMyznQ5ouR2Cac71VPqj7iiJCJu4EfADcAc4FMiMiei211AizFmOvAw8KC9dg6wEpgLLAceERF3HJsPAg8bY2qAFms75hhAN/D3wN9Emf6jwN1Ajf23PN77VRSlf3ptUkB8UUpv+C6YEh7pKQVFKS/bPeSeUmVJHqc7vXF6j24S8ZQWA7XGmAPGmF5gNbAios8K4En7+FngWuuVrABWG2N6jDEHgVprL6pNe8011gbW5i39jWGM6TDGvIkjTiFEZCJQaIx5yxhjgP8Ks6UoygAJbv7sr8wQOAkFafWUYqSEB0VpTLZnCD0lJ8mhsjSfM93e0HqXci6JiNJk4GjY8zrbFrWPMcYHtAJl/Vwbq70MOG1tRI4Va4z+5l0XZ94AiMjdIrJJRDY1NDT0Y1JRFG8SnlKqP+z7I3R0RUT4rjg/K/TY5RKyPa6UJzqEPKXSfIyBtm5fnCtGL4mIUrR1mEiZj9VnsNoTnUciczq30ZhVxphFxphF5eXl/ZgcnRij3+qUPoKeUnbclPD07lMKVXSIEKXSMdlAX7WJXE/qEzKCa0pTS/MBaO3SEF4sEhGlOqAy7PkUoD5WHxHxAEVAcz/XxmpvBIqtjcixYo3R37ynxJm3oihJEvwwTyjRIY3Zd/4YKeGLq0vJzXKx8jLnIygny53yta8z1jOqLHFE6XSXJjvEIhFR2gjU2Ky4bJzEhTURfdYAd9rHtwKv2nWcNcBKmzlXjZNs8E4sm/aa16wNrM3n44wRFWPMcaBNRK6wa1WfDbOlKMoACXpKOfFEKSvN+5RipIRPLMpj27eXccNFEwGnwkPKEx26vOR4XFQU5gDqKfWHJ14HY4xPRL4ErAPcwBPGmB0icj+wyRizBngceEpEanG8l5X22h0i8gywE/AB9xhj/ADRbNohvwGsFpHvAJutbWKNYW0dAgqBbBG5BbjeGLMT+ALwn0Ae8KL9pyjKeZBc9l36z1OKthMkfO65HnfqEx26vRTmZYWSLDQDLzZxRQnAGLMWWBvR9q2wx93AbTGufQB4IBGbtv0ATnZeZHt/Y1TFaN8EzIv2mqIoA6NvTcndb79071MKmHO9pGjkZg2BKHX5KMz1UGxFST2l2GhFByUumueghNOXEp5AmaE0Z98loEk2SzD1+5QK87IoVFGKi4qSoihJkWhKeCbsU0qkiEvuUCQ6dHkpzM0iN8tNbpZLRakfVJQURUmKxLPvXPgCBl+ajho35txjK6KRm+WiK+Weki/kJRXlZdGqa0oxUVFS4qLhOyWcYKJDvOy7PHt+UFeaQnj+QGLhu/xsD129qd3M2trlpdAWgi3Oy6ZF69/FREVJUZSkSDTRIT/b+RDu7E2fKMU78wmgINeT0goLxhgnfGc9peL8LM2+6wcVJSUuRg9GV8LoTTB8FzxpNV2i5PUH8CTgKo3N9dDWkzpR6vL68QUMhbmOKJXkq6fUHypKiqIkRa9NCkhUlDpS+IHfH/4EDiIEKMjx0OsLpCzZIViMtTDP8RxLxqgo9YeKkhIXXVNSwgmuKcVLCR+Tk97wnS9g8Ljif8SNtfPs6EmRKNlirH2ekhO+05qS0VFRUhQlKRI95C/kKaU4iSAW/oBJaPPsWCsW7SlaVwoWYy2wiQ4l+dn4AialIcPhjIqSEhf9PqeE05NglfBQokOKPJB4OJ5SAqJkPaW2ntQkHwT3JAVLDJXYKuUtHRrCi4aKkhIXDTMo4fR4/eRmueJuTE2/pxRIyFMKejAp85S6I0TJnufUohl4UVFRUhQlKbq8fvKy+k8Hh741pa50rSn5Ewzf2Xm2pyicFtwoWxRKCVdPqT9UlJS4qJ+khNPV6yc3AVFKv6dk4h7ZDmGeUqpEKZR954hS8JBBzcCLjoqSoihJkainlONx4ZL0rSl5E0x0KEzxcRKtXV7ys90hgdTwXf+oKClx0SUlJZxubyAhT0lEGJPtSeuaUiKJDiX52YhAU4rCaWe6vaHQHTip4S7R8F0sVJSUBFBVUvro9vpDde3ikZ/jTl/2XYJrSm6XUJKfTVN7T0rm0dp1tii5XEKxVnWIiYqSoihJ0WWz7xIhvZ5SYhUdAMrGZNOcIs+lNazuXZCS/KyUjTfcUVFS4qLhOyWcrt7E1pTA8ZTSVWbIFzC4E6joAE7yQVN7isJ3EZ4SwPjCXE6e6U7JeMMdFSVFUZKi25dY9h1AQU5WSitw94c/YMhK5OwKYNzYHJo6UhO+O93pDR2DHsQRpdSMN9xRUVLioo6SEk53Ep5SUV5W2k5Z9SWYfQfWU0pBOC0QMDS29zCuIOes9vGFuZxq69aN6VFISJREZLmI7BGRWhG5N8rrOSLytH19g4hUhb12n23fIyLL4tkUkWprY5+1mX0eY3xVRHaIyPsi8gsRyU3u9iiKEklXEokO6RQlfyCQ8JrSuLE5nO70Dnql8NYuL76AoXxspCjl4PUbXVeKQlxREhE38CPgBmAO8CkRmRPR7S6gxRgzHXgYeNBeOwdYCcwFlgOPiIg7js0HgYeNMTVAi7U9kDEmA38FLDLGzAPctp+SJPplTgkn0ZRwcI5rCJbZGWqc7LvEgkGTip3vq8dPD2ydx+sPRPV6Gm1GX3kUTwnQEF4UEvkfWwzUGmMOGGN6gdXAiog+K4An7eNngWvFKYy1AlhtjOkxxhwEaq29qDbtNddYG1ibtwxwDAAPkCciHiAfqE/g/SoR6CF/ShBjjM2+S9xT6vam7qyi/ki0ICvAlJJ8AI6d7kp6nEONHVz6jy9z33Pbz3mtoc0RnXHneEpWlNo02SGSRERpMnA07HmdbYvaxxjjA1qBsn6ujdVeBpy2NiLHSmoMY8wx4HvAEeA40GqM+W20Nygid4vIJhHZ1NDQEPNGKMpoJ1ghPJk1JSAtIbxEj64AmFKSB8CxluRF6bnNxzjT7WP1xqMcbe4867WGmJ6S8/yUZuCdQyKiFO1/NfKrc6w+g9We9BgiUoLjRVUDk4AxInJ7lL4YY1YZYxYZYxaVl5dH6zKq0fCdEiRYXDXRfUrB/Tln0iBKvkAg7kGEQSYU5eISqGvpjN85gg0Hmiiz9exe3nnyrNeOtzqiExShIBUFuYhA/QDDhSOZRH6z6oDKsOdTODcMFupjQ2VFQHM/18ZqbwSKrY3IsZId4zrgoDGmwRjjBZ4DPpDA+1UUJQbBjbDBCuDx6POUhj4tPNGKDgBZbhcTCnOpG4CndKipgw/PquCC8jGs33t2pOVIcyelY7IpyD07JTzb42JSUR5HmpMXwZFOIqK0EaixWXHZOMkCayL6rAHutI9vBV41zqrfGmClzZyrBmqAd2LZtNe8Zm1gbT4/wDGOAFeISL5de7oW2JXYbVHCUU9JCRKspF2QoCil01Pq9QXIdicWZgSoLh/D/ob2pMbo6vVz8kwPVWX5LJ1RwdsHms46quNocyeVpflRr51ams/hpo6kxhsNxBUlu37zJWAdzof6M8aYHSJyv4jcbLs9DpSJSC3wNeBee+0O4BlgJ/AScI8xxh/LprX1DeBr1laZtT2QMTbgJES8B2y373XVAO6RoiiW4EF4yXtKQy9KPf4AWZ7EPCWAGeML2HuynUAg8W9hQU9nWtkYPjyrnF5fgDdrG0OvH27qpNKuV0UyrSyfw03qKUWS0G+WMWYtsDai7Vthj7uB22Jc+wDwQCI2bfsB+rLnwtsHMsa3gW9Hu0ZJHM2+U4IEPaWxuYmJUnHoWIih3Y9jjKHXFyAngfOUgswcX0CX109dSxdTy6J7N5Ecsp7OtLJ8Zk0opCDHw8s7T/CROeM50+3lSHMnt106Jeq108rG0NTRS1u395zw3mhGKzoocdHwnRIk2fBdcX42rhQeCxELr9/5pc32JCFKEwoA2H3iTMLXBMNv00rHkO1xsXRWBa/sOoU/YNhxzLEzb0pR1GunWeFTb+lsVJQURUmYZMN3bpdQOiYntIl0qOj1O6nryYjSjPEFuF3C9mOtCV9zuKmTkvwsiuzBfR+dN4Gmjl5e3X2KP+5vxCWwsLI46rVVZWMAONCo60rhqCgpipIwyYbvAMaNzQ5tIh0qvHY/VXYS4bsxOR7mTCxk46HmhK853NTJVCsuAB+ZM57JxXn880u7eWbTURZXl1Kcnx312gsrxuBxCbuPJ+6ZjQZUlBRFSZigKI3JTlyUygtyaEjRsRCx6POUEs++A7isqpQtR0/Ta0UtHoeaOqgKW3/yuF1866Y51Da009DWw19eUxPz2hyPmwvLx7JLReksVJSUuOiakhKkvdtHfrY74f0/AOVjc2gcYk8pKCrJhO8ALqsqodsbYGvd6YTGqD/dxbSIlO9lcyfw4pc/yG/+6oMsmT6uXxuzJxaw63hbUnMc6agoKTFRLVIiae/xMTbB9aQgjqfUM6THNATLISVa0SHIkppxZLnlnMoM0Tja0knAQNW4Mee8NmtCIbMnFsa1MXtiISfOdNOi1cJDqCgpcdGUcCVIe48vqfUkcIqR9voCtA3hCbRBTyknSU+pMDeLKy8cx7odJ+KKaCjzruxcUUqUeZOdzLxEPLPRgoqSEhcN3ylB2rp9CaeDBxlX4Cz0nxrCYxoGkn0XZNnc8Rxu6oybhXeo0UnlrkpwT1M0FlQW43YJmw61DNjGSENFSYmJipESyenOXopiZJPFYmKRU9HgROvQFR8NrSklUWYoyE0XTyIvy81P3z7cb7/DTR0U5HgoHZPc/QhnTI6HeZMKeSdOxt+pM93c/+udo+JQQBUlJS6qTUqQlk4vJfnJVR8IHgsxkArcA2WgiQ7ghPBuWTiJNVvraepnf9Xek+1cUDEWp7TmwFlkM/76O3PqkfX7eeIPB3n8zQPnNdZwQEVJiYmuJSmRtHT2UpKkpzShMBe3SwZUgXug9PqdD/iBiBLAXVddQK8vwCPr90d93RjD+/WtzJ0UP5khHpdXl9LrC/BuPyG84JrTaAjzqSgpcRnKrCklc/H5A7R1+yhO0lPyuF1MLModUk+pqze5wwgjmV4xlj+5ZApPvXWYQ1EqLhxt7qKt28e8SdFLCCXDkunjyPa4+N2uU1FfDwQMe044aeO1p5KrYj4cUVFSYiL2/ESVJAXgtK30naynBE4Ibyg9pU577lN+9sBECeCvr59JTpaLv/nvrfgjKoe/dcCpBH7JtOglhJJhTI6HD1xYxiu7T0b9AnikuZPOXj/TK8bS1NFLa+fQV1wfSlSUlJho+E4JJ1jpO1lPCWBKSStHYvkAACAASURBVP6QilKX1wnf5Z2HKE0oyuUfbprLpsMtfPfFs49ie2XXKSYU5jJzfMF5zTPItbOdjL99UTyhYNsN8yYAsL9xZHtLKkpKXDR6p4CT5AAD85SmleZz4kw3HUO0V6nTHrR3Pp4SwCcumcxnr5zG///7g3xv3R4CAcOhxg5e2X2KmxdMOu8khyDL5o7H7RJ+tfnYOa8FQ3bXzh4PEDWcOJJIbsOBMqpQMVLCCVYdGIgo1YwfCzgfsBfHqJo9mARFKTfJ2neRiAjfvmku3V4/P3ytlpd2nAiVWvr8kurBmCoAFQW5LJ1RznPv1fE31888q4zT/oZ2KgpymD2xAJfAoRF+1IV6SkoCqDophPbIlIxJPnw3w4a59p4cmjpvnT2OcLiSqNEXC7dLePBP5vODlQsozc+malw+//m5xUwoyh2EmfZx66VTOHmmhzf2NZzVvr+hnQvLx5LjcTOpOG/EH6GunpISE5UiJZxTtqhqeUFO0tdOK3MOwYu2ZpIKOr3+8w7dhSMirFgwmRULJg+azUiumV1BRUEOP/n9AT48swJwMu9qT7azYuEkwDmDST0lZdSjYTwF4OSZbkrys8gZQEjM7RIuLB8bSm1ONV29/vNKckgHOR43d11VzR9qm9h61NmXdKCxnbYeHxdPcUKe08ryR7ynpKKkxEU1SQHHUxpfOPCQ1eyJBeyobx2SfW+dvT7ys4ZfIOjTl0+lKC+L7764G2MMG+1m2YVTHVGqKhvD6U5vKBNyJJKQKInIchHZIyK1InJvlNdzRORp+/oGEakKe+0+275HRJbFsyki1dbGPmsz+zzGKBaRZ0Vkt4jsEpErk7s9iqIEOXWme0ChuyCXTC2hsb2Xo82pTw3v6PGTnzO8PCWAgtws/mbZTN460MTP3znCmi31TC3N58JyJ1Fkmi3+OpJDeHFFSUTcwI+AG4A5wKdEZE5Et7uAFmPMdOBh4EF77RxgJTAXWA48IiLuODYfBB42xtQALdZ20mPYa34AvGSMmQVcDJy92UBJCA3fKQAnz5yfp3TptBIA3j2S+HHjA6W1y0txXvIJGZnApxdP5YM14/jmr97nrQNN3HHFtFDqebU9u2kkh/AS8ZQWA7XGmAPGmF5gNbAios8K4En7+FngWnHu4gpgtTGmxxhzEKi19qLatNdcY21gbd4ykDFEpBD4EPA4gDGm1xijh5YkgYqREiQQMDS091BxHp7SjPEFjM3xDEn9ttNdvRQNU1Fyu4Qf33Epf7tsJn9342w+f1Vf6nllaT4ifcdmjEQSCbpOBo6GPa8DLo/VxxjjE5FWoMy2vx1xbTB9JZrNMuC0McYXpX+yY3QBDcB/iMjFwLvAl40x53zFEJG7gbsBpk6dGus+jFq09p3S0N6DP2DOKw3a7RIuqyrhzdpGjDGDtvE0Gq2dXooHsJ8qU8jP9nDPh6ef056b5WZiYe6o95Si/eZEfkrF6jNY7QMZwwNcAjxqjFkIdADnrIcBGGNWGWMWGWMWlZeXR+sySlExUhyONDvfzKeWDvxAO4BrbDmdAymsSuAPGM50+ygcpp5SPKaVjeHQKBelOqAy7PkUoD5WHxHxAEVAcz/XxmpvBIqtjcixBjJGnTFmg21/FkeklCRRaVKONA2SKM1y9t+8suvkec8pFmds4djhuqYUj6px+RwezYkOwEagxmbFZeMkFayJ6LMGuNM+vhV41TgxnzXASps5Vw3UAO/Esmmvec3awNp8fiBjGGNOAEdFZKa95lpgZwLvV4lAo3fKkeZORJzCqufD5OI85k0u5H82R36vHTyC1cyH65pSPKaVjaGpo5cz3SOzWnhcUbLrO18C1uFkrz1jjNkhIveLyM222+NAmYjUAl/DhsmMMTuAZ3DE4CXgHmOMP5ZNa+sbwNesrTJrO+kx7DV/CfxMRLYBC4B/SvYGKYoCR5s7mVSUN+BD88L500WV7Dx+hvePtQ7CzM7l1Bnn2PXzSV/PZKqCaeEjtDBrQrvLjDFrgbURbd8Ke9wN3Bbj2geABxKxadsP4GTnRbYPZIwtwKJo1yiJo0dYKIebO6kszRsUWysunswDv9nFU28d5sFb5w+KzXBOWFGaOMi16TKFYB3B3cfbmD8l9cVthxqt6KDERMN2CjjZl3tPtlFTMThnBxXlZ7Hyskp++V4dR5sHf23kpBWl8SNUlKrKxjA2x8P2FHma6UZFSYmPitOopr61m7ZuHzMnDI4oAXxh6XRcLuGf1+0ZNJtBTrT2kJ/tpiBn+JUZSgSXS5gzqVBFSRl9mIifyuhkz4kzAMwaRFGaUJTLF5deyK+31vO7nYObiXe0pZNJxXkp3QeVbi6aXMSu42fw+QPpnsqgo6KkKEq/7DruVPaeMYiiBPDFpdOZNaGAv/7vrRwcxEX7/afamW5rxY1ULppcRI8vQG3DyDsaXUVJiYuuLY1utte1MrU0n8LcwU2xzva4WHXHIlwCdz7xzqCsL/X6Ahxu7mR6xcgWpeDpve8dHnmV01SUFEWJiTGGTYdbQsVUB5upZfk88b8uo7XLy62P/ZF3Dp5fsda9J9vwB8yge3WZRlVZPuUFObxzsCndUxl0VJSUuGhK+OjlaHMXje09KRMlgIVTS3jmz68kN8vNylVv8Z0XdtLaNbCNoUFRu6wqdfPNBESExdWlbDjYPOJqU6ooKTEZab/sSvJsPOR8yKdSlABmTijgN3/1Qf50USWP/+EgS//lNf7tlX00tfckZefN2kYqS/OYWDQ4e6oymcurSzne2k1dS+rPpxpKVJSUuKg2jV7W721g3NgcZo5PfThsbI6H7/7JfF74y6tYUFnMQy/v5QPffZVvPLstoeoPzR29vLG3gY/Om5jyuWYCi6tLAc475JlpjMxEfmVQUU0anfj8AV7fc4plcyfgcg1devXcSUX8x+cWs+9kG0/84SC/2nyMpzcd5ZKpxXz2yipuuGgCOZ5zT5V94s2D+AKGP7l0ypDNNZ3MqCigJD+LP+5vGlHvWT0lJSYqRqObjYdaONPt49rZFWkZv2Z8Af/3E/PZcN91/N2Ns2nu6OUrT29hyXdf5V/W7ab+dF/YasOBJn78xn5uunhSqAzPSMflEpZMH8fv9zWMqFC7ekpKXEbSL7ySOM+9V8eYbDcfmpHeM8aK8rP4sw9ewOeXVPP72kaeeusQj6zfz6Pr93Pd7PF8+vKpfOXpLUwtzec7K+alda5DzYdmlPPCtuPsPtHG7ImF6Z7OoKCipMREtWj00tHj4zfbj3PT/EnkZ2fGx4TLJVw9o5yrZ5RztLmTn204wn/84SC/3XmSbI+L5++5jKL8kXlcRSw+VON8YXhjb8OIESUN3ylxUW0afTz3Xh2dvX5uW5SZaxWVpfnce8MsfnKncwjAv9w6n2llY9I8q6FnQlEuM8cX8Ma+hnRPZdDIjK9AiqJkDF5/gMdeP8AlU4tTngp+vnywppyt375+xB7olwgfmjGOJ/94mM5eX8Z4teeDekpKfNRVGlU8s+kox0538cWl04dFUdPRLEjgrCv1+gNsODAyUsNVlJS4aEWH0UNLRy//sm4Pl1eXpi3rTkmOy6pKyc1y8frekRHCU1FSYqJZd6MLYwzfXrODtm4f/3Dz3GHhJSmQm+XmygvKeENFSRktqDaNDp7eeJQ1W+v56nU1IyaTa7Rw9YxyDjR2cKRp8E/yHWpUlJSYqBaNHt7Y28DfP/8+V00fxxeWTk/3dJQkuXqmE2p9fQRk4SUkSiKyXET2iEitiNwb5fUcEXnavr5BRKrCXrvPtu8RkWXxbIpItbWxz9rMHugY9jW3iGwWkRcSvy0KEFIl9ZRGNm8faOIvfvou0ysKeOT2S3APYUkhZXCoKstnamk+r+85le6pnDdxRUlE3MCPgBuAOcCnRGRORLe7gBZjzHTgYeBBe+0cYCUwF1gOPGJFoj+bDwIPG2NqgBZrO+kxwub2ZWBXYrdDCSegajTieXH7cT77xDtMKs7jyc9dNugH+SlDg4izsfiP+5vo8fnTPZ3zIhFPaTFQa4w5YIzpBVYDKyL6rACetI+fBa4VZ5V0BbDaGNNjjDkI1Fp7UW3aa66xNrA2bxngGIjIFOBG4CeJ3Q4lnEDQU0rvNJQU4PUH+O6Lu/nCz95j3qRC/vvPr6SiMDfd01LOg6tnlNPZ6+fdQy3pnsp5kYgoTQaOhj2vs21R+xhjfEArUNbPtbHay4DT1kbkWMmOAfB94OtAoL83KCJ3i8gmEdnU0DD8Y7KDRTAVXLPwRhb7G9pZueptHnt9P59aPJWf/+8rKBmTne5pKefJlReWkeWWYZ8anogoRQswR35KxeozWO1JjyEiHwNOGWPejfL62Z2NWWWMWWSMWVRent7ik5lEQLVoRNHj8/P93+3lhu//nn0n2/jBygX8309cRG7WucdAKMOPMTkeLqsqHfailEhNijqgMuz5FKA+Rp86EfEARUBznGujtTcCxSLisd5QeP9kx7gZuFlEPgrkAoUi8lNjzO0JvGeFPg9JtWl4EwgYfr2tnode3svhpk5WLJjE3904h/KCnHRPTRlkls4s55/W7uZEazcTioZnODYRT2kjUGOz4rJxkgrWRPRZA9xpH98KvGqcT7Q1wEqbOVcN1ADvxLJpr3nN2sDafH4gYxhj7jPGTDHGVFn7r6ogJYdG7YY3xhjW7znFx/79Tb68egv52R6eumsxP1i5UAVphHL1DCc1fDhvpI3rKRljfCLyJWAd4AaeMMbsEJH7gU3GmDXA48BTIlKL472stNfuEJFngJ2AD7jHGOMHiGbTDvkNYLWIfAfYbG0zkDGU8yOYfafiNLzwBwwvvX+Cx17fz/ZjrVSW5vH9Ty7g5osnDekJssrQM2P8WCYU5rJ+7yn+9LLK+BdkIAmVlDXGrAXWRrR9K+xxN3BbjGsfAB5IxKZtP4DNnotoT3qMsNfXA+tjva5Ep29NSVVpONDt9fPce8dY9cZ+DjV1UlWWzz99/CJuvXQK2R7dJz8aCKaGr33/OD5/AI97+P2/D/8654oyyqk/3cXPNxxh9cYjNLb3Mn9KEY985hKWzZ2gG2FHIVfPLOfpTUfZcvQ0i6pK0z2dpFFRUuKi4bvMwxjDH2qb+K+3DvG7XScxwLWzKvjckmo+cGGZFlMdxSyZPg63y0kNV1FSFCWlnOn28st363jq7cMcaOigJD+Luz90IZ+5fCqVpfnpnp6SARTlZbGwspjX9zbw19fPTPd0kkZFSYmLOkrpJRAwbDjYzH9vOsra94/T7Q2woLKYh/70Yj560UTdZ6Scw9UzyvnXl/fS2N7DuLHDK9NSRUlRMpTjrV388t06ntlUx5HmTgpyPPzJJVNYedlULppSlO7pKRnM1TMdUfr9vgY+vnBKuqeTFCpKSlx0TWno6PH5eWXXKZ7eeJTf72sgYODKC8r46kdqWD53InnZ6hUp8Zk3qYiS/Cze3NekoqSMPPQ49NRijGHn8TM8+24d/7P5GC2dXiYW5XLPh6dz26WVTC3TtSIlOVwu4YoLynj7QBPGmGGV+KKipChpov50F89vqedXm+vYe7KdLLdw/ZwJ3LZoCh+sKdd0buW8uPLCMl58/wRHm7uG1RcbFSUlLhq+GzzOdHt5afsJfrX5GG8fbMIYuHRaCf94yzw+dtFErdatDBpXXlAGwFsHGplaNjXNs0kcFSVFSTFef4DX9zTwqy3H+N3Ok/T4AlSV5fOVa2dwy8JJTCsbk+4pKiOQ6RVjGTc2h7f2N/HJy1SUlBGEOkrJY4xhy9HT/GrzMV7Ydpzmjl5Kx2Sz8rJKblk4mQWVxcMqzq8MP0SEKy4o5a1htq6koqQog8jek238ems9v95az6GmTnI8Lq6bM55PLJzMh2aUkzUMa5Epw5crLyzjhW3HOdjYwQXlY9M9nYRQUVLioifP9s/Bxg5e2FrPr7fVs/dkOy6BD1w4ji8unc7yiyZQmJuV7ikqo5S+daUmFSVFGckcO93FC1vreWHbcbYfawXgsqoS7l8xlxvmTdTzipSMoHrcGMYXOutKn7l8WrqnkxAqSkpUAnoW+jmcautm7bbj/Hrbcd493ALAxVOK+LsbZ/PRiyYyqTgvzTNUlLMREa68oIw3a4fPupKKkhIVbyAQejyao3fNHb289P4Jfr21PpTCPWtCAX+7bCYfmz9RM+eUjOeKC8r4ny317G9oZ3pFQbqnExcVJSUqXv/oVaL60138dscJ1u04yTuHmvEHDBeMG8NfXlPDTfMnUjM+8/+wFSXIZdXO8RWbDrWoKCnDF68vzFMaBUnhtafaWbfjBL/dcYKtdc4a0fSKsfzF1Rdww7yJzJ1UOCxCH4oSyQXjxlA6JptNh1tYuTjz9yupKI1wHl2/n0unlbC4OrnDvsLDdyMRYwzb6lpZt+ME63acYH9DB+CsEX19+UyWzZ3AhcMkW0lR+kNEuHRaCZsONad7KgmhojTCefCl3QAc+u6NSV0XHr4bKWtKPn+Adw4189sdJ/ntjhPUt3bjdgmXV5fy2SuruH7ueCYWabKCMvK4rKqEl3eepKGtJ+MzQxMSJRFZDvwAcAM/McZ8N+L1HOC/gEuBJuCTxphD9rX7gLsAP/BXxph1/dkUkWpgNVAKvAfcYYzpTXYMEam0/ScAAWCVMeYHyd6g4YzXP3Bvp9vrDz0ezqLU7fXz+32NrNtxgld2naSl00uOx8WHZpTztetncu2sCq03p4x4gseiv3u4meXzJqZ5Nv0TV5RExA38CPgIUAdsFJE1xpidYd3uAlqMMdNFZCXwIPBJEZkDrATmApOA34nIDHtNLJsPAg8bY1aLyGPW9qMDGMMH/LUx5j0RKQDeFZGXI+Y9omnv9g342o6egV+bbs50e3lt9ynW7TjB+j0NdPb6Kcj1cO2sCpbNncDVM8vJz9YggTJ6mDepiByPi42HWoa/KAGLgVpjzAEAEVkNrADCP9xXAP9gHz8L/FCcVeEVwGpjTA9wUERqrT2i2RSRXcA1wKdtnyet3UeTHcMY8xZwHMAY02ZtT46Y94im/TyEpaMnzFMajMmkmFNt3by88yTrdpzkrf2NeP2G8oIcPr5wMsvmTuCKC8rI9miJH2V0ku1xcXFl8bBYV0pElCYDR8Oe1wGXx+pjjPGJSCtQZtvfjrh2sn0czWYZcNoY44vSfyBjACAiVcBCYEO0NygidwN3A0ydmvnZKYnSFRaCS5bO3sz3lI40dYYSFd490oIxMK0sn88tqWbZ3PEsrCzBpWcSKQrgrCv9+PUDdPb6MjpSkMjMov1VR355jtUnVnu0r6z99R/IGM5FImOBXwJfMcacidIXY8wqYBXAokWLhoNjkBDnE4Lr6A1fU8qMW2KMYdfxtpAQ7T7RBsCciYV85doZLJs3npnjCzR1W1GisKiqlB+9tp8tR0/zgQvHpXs6MUlElOqAyrDnU4D6GH3qRMQDFAHNca6N1t4IFIuIx3pL4f2THkNEsnAE6WfGmOcSeK8jik4rLAM5wfRMlzf0OJ2SFAgY3jvSYoXoJEeaOxGBRdNK+LsbZ7Ns7gQqS4fPqZqKki4umVqCiLOJdriL0kagxmbFHcNJKvh0RJ81wJ3AW8CtwKvGGCMia4Cfi8hDOEkINcA7ON7NOTbtNa9ZG6utzecHMoZdb3oc2GWMeSjZGzMSCHpK2QM4LqGxvWewp5Mwvb4Af9zfyLodJ3l550ka23vIcgtLpo/jC0sv5LrZ4zM+rVVRMo2ivCxmji9gY4avK8UVJbt+8yVgHU769hPGmB0icj+wyRizBufD/ymbZNCMIzLYfs/gJBf4gHuMMX6AaDbtkN8AVovId4DN1jbJjiEiVwF3ANtFZIu18f8ZY9YO7FYNP4Ke0kAW+BvawkRpCFyljh4fr+9t4KX3T/Da7lO09fjIz3bz4ZkVXD93PB+eVaFHQCjKeXLptBKe31KPP2AGFEEZChJa7bIf5Gsj2r4V9rgbuC3GtQ8ADyRi07YfoC9DL7w9qTGMMW8Sfb1p1NBqQ3B5We6krz3V1kNBjoe2FKaGN3f08rtdzkbWN/Y10usLUJKfxQ0XTWDZ3AksmT6O3AHMXVGU6CyqKuFnG46w50QbcyYVpns6UcncFAzlvAmG4Apyk/9v3t/QTnX5GLbVtQ5q7btgsdOXdpzgnYPNBAxMKsrl04unsmzuBC6rKsGjp7MqSkpYNK1vE62KkjLkNLb3AuBPMnuu2+vncFMntyyYzDZbnPR8ONzUwYvvn+DF90+w9ehpAGoqxvLFpdNZNncC8yZrsVNFGQqmlORRUZDDpsMt3HFlVbqnExUVpRHMsdNdAPiSPIZioz2uYeHUYn75Xt2AygzVnmrjxe0nWPv+CXYddzLxL5pcxN8um8nyeVrsVFHSgYiwqKqETYda0j2VmKgojWD2nXT28fiSrIG3dvtxcjwurrywDEg8z+Focyf/s/kYz2+tp/ZUO+AsrGrqtqJkDpdOK2Xt9hOcaO1mQlFuuqdzDipKI5QjTZ0cb+0GwJvE0eaHmzp47r1jfOKSyaG1KF8/17d1e1mztZ7/2XyMjfbb1+LqUu5fMZdlcycwvjDzfukVZTSzaFoJAJsON/Ox+ZPSPJtzUVEagRhjeGR9LS6Bq2eUs9mu48TjTLeXL/18M9luF391bQ05bifzLfzAvyAHGtr5zz8e4pfv1tHR66emYixfXz6TFQsmM7lYj39QlExlzqRC8rLcbDrUoqKkpBafP8Bvth/n0fX72X2ijf/9wWp8ARPyYPpj38k2vvTzzRxobOfHd1zKxKK8UP273rDw3/HWLv71t3t57r06PC4XH7t4IndeWcX8KUWarKAow4Ast4uLK4t493BmriupKI0AfP4Az2+p54ev1XKwsYOairH8620X8/GFk3nwpd39nqvU1u3lkfX7efzNgxTkeHjif13GB2vKgb5KEEFP6YVt9dz33HZ6fAHuuqqaP7/6QsaN1coKijLcWDStlEdf309Hj48xOZklA5k1GyUpfP4Av9p8jB++Vsvhpk7mTCzksdsv4fo5E0LVsT1uibom5PUHWP3OEb7/u300dfTy8YWTue+GWVSErQG5XYKI4yn94p0j3Pfcdi6ZWszDn1zAtLIxQ/Y+FUUZXC6tKsH/mmHr0dN8YHpm1cFTURqGeINi9GotR5o7mTupkFV3XMpH5ow/J4TmcbnwBwzGGEQEYwxrt5/ge7/dw8HGDhZXl/IfN85m/pTic8YREbLdLjYfOc2j6/ezdGY5P77jUnI8WmVBUYYzl1QGkx1aVJSUgdPU3sPqjUf5+YYjHDvdxUWTi/jJZxdx7eyKmOs5WW6n3es31J/u5OvPbuOdQ83MGD827rXg1M17s7aRcWOz+cEnF6ogKcoIoCg/ixnjx7IpA9eVVJQyHH/AsOFgE89uquOFbcfp9QdYMr2M+1fM5ZpZ/QsKECrZs+FgE1/86Xsg8N1PXMRtiyoTKsjYZo9U//xV1RTla0FURRkpXDqtlBe2Zl5xVhWlDMQYw9a6VtZsqeeFbfWcauthTLabTy2u5I4rpzG9oiBhWx77y/bFn71H2dhsfvpnlzOlJPFNrAsqi9ly9DSfWTwt6fehKErmsmhaCb945wh7T7Yxe2Lm1MFTUcog9p5sY82Wen69rZ7DTZ1ku11cPbOcmy+exLWzKwZ0hHGW9ZTaun388NOXJCVIAP+2ciGdXp96SYoywlhU1beupKKkhDja3MmarfX8ems9u0+04RL4wIXjuGfpdJbNm0BR3vmJQY49S2lKSR4fqkl+QXNqmZYGUpSRyNTSfMaNzeHdQ83ccUXmREJUlNJAQ1sPa7cf5/ktx3jviFNt4ZKpxfzDTXP46PyJVBQMXmmeJdPHUVMxln+57WLd3KooSggRYdG0koxLdlBRGiK6ev2OEG2t5w+1jfgDhlkTCvj68pncNH9SyoqVVpbm8/LXrk6JbUVRhjeLqkp4accJTp7pzpg6lSpKKWZ/Qzs/ffswz75bR1u3jyklefzF1Rdw88WTmTkh8YQFRVGUwebSYHHWQy3cOH9immfjoKKUImpPtfODV/bxwrZ6PC7hhnkT+czlU1lcXaphNEVRMoK5k4rI8bjYdLhZRWmk0uPz88NXa3l0/X6yPS6+cPWFfG5JNeUFWiNOUZTMItvjYkFlMRsONKd7KiFciXQSkeUiskdEakXk3iiv54jI0/b1DSJSFfbafbZ9j4gsi2dTRKqtjX3WZvZgj5EqGtp6WLnqbf791VpuXjCJN77+Yb6+fJYKkqIoGcs1syrYefwMR5s70z0VIAFREhE38CPgBmAO8CkRmRPR7S6gxRgzHXgYeNBeOwdYCcwFlgOPiIg7js0HgYeNMTVAi7U92GMMOk3tPXzyx2+x6/gZHv3MJTz0pwu0graiKBnPDfOcsN1vth9P80wcEvGUFgO1xpgDxpheYDWwIqLPCuBJ+/hZ4FpxFk5WAKuNMT3GmINArbUX1aa95hprA2vzlsEcI7Hbkhy9vgCff3ITx0538dRdl3PDRZkRm1UURYnH1LJ8FleX8u+v7OOEPa06nSSypjQZOBr2vA64PFYfY4xPRFqBMtv+dsS1k+3jaDbLgNPGGF+U/oM1xjmIyN3A3fZpu4g0AY3R+sZj8QMDuSpjGccA78MIRO9FH3ovHEbcfZj4jwO+dBwwKDtwExGlaKlikQf0xOoTqz2ah9Zf/8Ec49xGY1YBq4LPRWSTMWZRtL6jCb0Pfei96EPvhYPehz7svagaDFuJhO/qgMqw51OA+lh9RMQDFAHN/Vwbq70RKLY2IscarDEURVGUDCURUdoI1NisuGycpII1EX3WAHfax7cCrxpjjG1faTPnqoEa4J1YNu01r1kbWJvPD+YYid0WRVEUJR3EDd/Z9ZsvAesAN/CEMWaHiNwPbDLGrAEeB54SkVoc72WlvXaHiDwD7AR8wD3GGD9ANJt2yG8ATTpK0wAAA9BJREFUq0XkO8Bma5tBHiMeq+J3GRXofehD70Ufei8c9D70MWj3QhxnQ1EURVHST0KbZxVFURRlKFBRUhRFUTIGFaUwhrosUToQkSdE5JSIvB/WVioiL9vSTi+LSIltFxH5N3s/tonIJWHX3Gn77xORO6ONlcmISKWIvCYiu0Rkh4h82baPxnuRKyLviMhWey/+j20ftJJfwwlbEWaziLxgn4/W+3BIRLaLyBYR2WTbUv/3YYzRf866mhvYD1wAZANbgTnpnlcK3ueHgEuA98Pa/hm41z6+F3jQPv4o8CLOXrArgA22vRQ4YH+W2Mcl6X5vSd6HicAl9nEBsBenHNVovBcCjLWPs4AN9j0+A6y07Y8BX7CPvwg8Zh+vBJ62j+fYv5scoNr+PbnT/f4GcD++BvwceME+H6334RAwLqIt5X8f6in1MWRlidKJMeYNnOzFcMJLOEWWdvov4/A2zh6yicAy4GVjTLMxpgV4Gafu4LDBGHPcGPOefdwG7MKpBDIa74UxxrTbp1n2n2HwSn4NG0RkCnAj8BP7fDBLn40EUv73oaLUR7RySpNj9B1pjDfGHAfnwxqosO2x7smIulc27LIQx0MYlffChqy2AKdwPjj2k2DJLyC85NdwvxffB74OBOzzhEufMbLuAzhfTH4rIu+KU4oNhuDvQ89T6iORckqjjWRLOw07RGQs8EvgK8aYMxL7AMYRfS+Ms7dvgYgUA78CZkfrZn+OyHshIh8DThlj3hWRpcHmKF1H9H0IY4kxpl5EKoCXRWR3P30H7V6op9THaC5LdNK62tifp2z7iC7hJCJZOIL0M2PMc7Z5VN6LIMaY08B6nHWBwSr5NVxYAtwsIodwwvfX4HhOo+0+AGCMqbc/T+F8UVnMEPx9qCj1MZrLEoWXcIos7fRZm1lzBdBqXfZ1wPUiUmKzb663bcMGG/t/HNhljHko7KXReC/KrYeEiOQB1+GssQ1Wya9hgTHmPmPMFOMUFl2J874+wyi7DwAiMkZECoKPcX6v32co/j7SneGRSf9wMkj24sTTv5nu+aToPf4COA54cb7F3IUTB38F2Gd/ltq+gnNQ4n5gO7AozM7ncRZwa4HPpft9DeA+XIUTRtgGbLH/PjpK78V8nJJe2+wHz7ds+wU4H6a1wH8DObY91z6vta9fEGbrm/Ye7QFuSPd7O497spS+7LtRdx/se95q/+0Ifh4Oxd+HlhlSFEVRMgYN3ymKoigZg4qSoiiKkjGoKCmKoigZg4qSoiiKkjGoKCmKoigZg4qSoiiKkjGoKCmKoigZw/8DzKwcax3qENcAAAAASUVORK5CYII=\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 = 25\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",
" _1 = np.where((total_samp >= x_min) & (total_samp <= (jpsi_mass - 50.)))\n",
" \n",
" tot_sam_1 = total_samp[_1]\n",
" \n",
" _2 = np.where((total_samp >= (jpsi_mass + 50.)) & (total_samp <= (psi2s_mass - 50.)))\n",
" \n",
" tot_sam_2 = total_samp[_2]\n",
"\n",
" _3 = np.where((total_samp >= (psi2s_mass + 50.)) & (total_samp <= x_max))\n",
" \n",
" tot_sam_3 = total_samp[_3]\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, 1.5e-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": 40,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"2/2 fits converged\n",
"Mean Ctt value = -0.41593044149928\n",
"Mean Ctt error = 0.16600757990339696\n",
"Sensitivy = 0.00011574576965860746\n"
]
}
],
"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)))\n",
"print('Sensitivy = {}'.format(np.mean(Ctt_error_list)**2*4.2/1000))"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAR/0lEQVR4nO3dbYxcV33H8e+/dh5QodhJNlFkm24QVkVelBCtUkupqjYhkAeE8yJIQaixqCVLJZVAVKJGSJWQ+iLpC0IjVSBDUJ2Kh9AAipVAqeUkQpWagENCSOqCN27arGzFhiQGVEEb+PfFnI2H9ezO7DzPud+PNJp7zz0zc+7Zmd+cPXPnTmQmkqS6/NakGyBJGj7DXZIqZLhLUoUMd0mqkOEuSRXaOOkGAFx00UU5Pz8/6WZI0kx54oknfpyZc522TUW4z8/Pc/jw4Uk3Q5JmSkT812rbnJaRpAoZ7pJUIcNdkipkuEtShQx3SaqQ4S5JFTLcJalChrskVchwl6QKGe5SReb3PjTpJmhKGO6SVCHDXZIqZLhLUoUMd0mqkOEuSRUy3CWpQoa7JFXIcJekChnuklQhw12SKmS4S1KFDHdJY+c5cEbPcJekChnuklQhw12SKmS4S1KFDHdJqpDhLkkV6incI+L5iPhBRDwVEYdL2QURcTAijpbrzaU8IuLuiFiMiKcj4spR7oAk6WzrGbn/SWZekZkLZX0vcCgztwOHyjrADcD2ctkDfHpYjZUk9WaQaZmdwP6yvB+4ua383mx5DNgUEZcO8DiSpHXqNdwT+JeIeCIi9pSySzLzBEC5vriUbwFeaLvtUimTJI3Jxh7rXZ2ZxyPiYuBgRPzHGnWjQ1meVan1JrEH4E1velOPzZAk9aKnkXtmHi/XJ4GvA1cBLy5Pt5Trk6X6ErCt7eZbgeMd7nNfZi5k5sLc3Fz/eyBJOkvXcI+I346INywvA+8EngEOALtKtV3AA2X5AHBbOWpmB3B6efpGkjQevUzLXAJ8PSKW638xM/85Ir4LfCUidgP/Dby31P8GcCOwCPwP8IGht1qStKau4Z6Zx4C3dSj/CXBth/IEbh9K6yRJffEbqpJUIcNdkipkuEtShQx3SaqQ4S5JFTLcJalChruksZrf+9Ckm9AIhrskVchwl6QKGe6SVCHDXZIqZLhLUoUMd0mqkOEuSRUy3CWpQoa7JFXIcJekChnukibC0xCMluEuSRUy3CWpQoa7VJlpnu6Y5rbVxnCXpAoZ7pJUIcNdkipkuEtShQx3SaqQ4S5JFTLcJalCPYd7RGyIiCcj4sGyfllEPB4RRyPivog4t5SfV9YXy/b50TRdkrSa9YzcPwQcaVu/E7grM7cDLwO7S/lu4OXMfAtwV6knSRqjnsI9IrYCNwGfK+sBXAPcX6rsB24uyzvLOmX7taW+JI1dU78V2+vI/VPAR4Ffl/ULgVcy89WyvgRsKctbgBcAyvbTpf5viIg9EXE4Ig6fOnWqz+ZLkjrpGu4R8W7gZGY+0V7coWr2sO1MQea+zFzIzIW5ubmeGitJ6s3GHupcDbwnIm4Ezgd+h9ZIflNEbCyj863A8VJ/CdgGLEXERuCNwEtDb7kkaVVdR+6Z+bHM3JqZ88CtwMOZ+X7gEeCWUm0X8EBZPlDWKdsfzsyzRu6SNE5Nm3sf5Dj3vwI+EhGLtObU7ynl9wAXlvKPAHsHa6Ikab16mZZ5TWY+Cjxalo8BV3Wo8wvgvUNomySpT35DVVL1mjYlA4a7JFXJcJdmSBNHoOqP4S5VyDcBGe6SVCHDXZIqZLhLUoUMd0mqkOEuSRUy3CU1RpOOIjLcJalChrukiWnSSHrcDHepUgZnsxnuklQhw13SRPkfxmgY7pKmgiE/XIa7JFXIcJdUpab/J2C4S1KFDHdJqpDhLmnimj6FMgqGu6RGacobieEuSRUy3CVVqymj9E4Md0mqkOEuSRUy3CXNrCZPu3RjuEtShbqGe0ScHxHfiYjvR8SzEfGJUn5ZRDweEUcj4r6IOLeUn1fWF8v2+dHugqSm6DRSd/TeWS8j918C12Tm24ArgOsjYgdwJ3BXZm4HXgZ2l/q7gZcz8y3AXaWeJGmMuoZ7tvy8rJ5TLglcA9xfyvcDN5flnWWdsv3aiIihtVhqKEeoWo+e5twjYkNEPAWcBA4CzwGvZOarpcoSsKUsbwFeACjbTwMXdrjPPRFxOCIOnzp1arC9kKSilzfBJrxR9hTumfmrzLwC2ApcBby1U7Vy3WmUnmcVZO7LzIXMXJibm+u1vZI0dDWG/bqOlsnMV4BHgR3ApojYWDZtBY6X5SVgG0DZ/kbgpWE0VpJWU2NAD6KXo2XmImJTWX4d8A7gCPAIcEuptgt4oCwfKOuU7Q9n5lkjd0lajUE9uI3dq3ApsD8iNtB6M/hKZj4YEf8OfDki/gZ4Erin1L8H+MeIWKQ1Yr91BO2WJK2ha7hn5tPA2zuUH6M1/76y/BfAe4fSOklnmd/7EM/fcdOkmzF2a43m+x3p19yXfkNVkipkuEuaGtM81z7NbevEcJc0kwYN21kL6/Uy3CWpQoa7JPVolkb7hrtUsVkKo2Xd2tzLUTOzuN/DZrhLUoUMd2kGOBKdjFnud8NdqtwsB9SySe/DpB+/H4a7JFXIcJekDmZxtN7OcJekChnuklQhw13SVJqWaZFpacd6Ge6Sptqow3VWw7sbw12SKmS4S5pakxpV1zCaN9wlqUKGu6SZU8PIetQMd6mBDMf6Ge6SVCHDXZIqZLhLEvVNVRnuUgP4C0XNY7hPIV+Aajeu54PPu97MSj8Z7pJUIcNdaohZGXFqOAx3SapQ13CPiG0R8UhEHImIZyPiQ6X8gog4GBFHy/XmUh4RcXdELEbE0xFx5ah3QtL6DWsk738E06mXkfurwF9m5luBHcDtEXE5sBc4lJnbgUNlHeAGYHu57AE+PfRWS9IUmOY3tq7hnpknMvN7ZflnwBFgC7AT2F+q7QduLss7gXuz5TFgU0RcOvSWSxqpaQ6uUalpn9c15x4R88DbgceBSzLzBLTeAICLS7UtwAttN1sqZSvva09EHI6Iw6dOnVp/y6XK1RQ0Gr+ewz0iXg98FfhwZv50raodyvKsgsx9mbmQmQtzc3O9NkOS1IOewj0izqEV7F/IzK+V4heXp1vK9clSvgRsa7v5VuD4cJorNcM4R+0rH2u1x24v97+K6dfL0TIB3AMcycxPtm06AOwqy7uAB9rKbytHzewATi9P30ham6E5ebX8DXoZuV8N/ClwTUQ8VS43AncA10XEUeC6sg7wDeAYsAh8Fvjg8JstadJqCcFabexWITP/lc7z6ADXdqifwO0DtkvSGA3zmPfn77hpKPelwfgNVWmKjWJ0PMjcea/z800w7fvedeQuaXTaA2IaRryDBNa0h13TOHKXpAoZ7tIYzcrodtDpmtrNwv4a7tKUmIXA0Oww3CVpQNP4xmy4S2PWLQimLSimrT3qjeEujUhNoVjTvjSF4S5JFTLcpQlwJKxRM9ylHqwnjPs5y6Jm37T9PQ13qQ/T9kKWVjLcJalChrsaq9/Rdz/TLLM60u/U7lndl6Yx3KUuDDPNIsNd6pOhr7VM+vlhuEvrMMg0xaRf7BqNaf1tWcNdjTdNL0jNvml5PhnumlnTeI6W5cdceS2Nm+EuSRUy3KUVHG2rBoa7RG8flBr6miWGuxrBYFbTGO7SEPjmoWljuKuRej2FgKGtfkzD88ZwVzWm4QUlTQvDXY3mG4JqZbirMeb3PtTzmRs9pYBmXddwj4jPR8TJiHimreyCiDgYEUfL9eZSHhFxd0QsRsTTEXHlKBuv0Zu28Jq29kjTqpeR+z8A168o2wscysztwKGyDnADsL1c9gCfHk4zVYPVjiUfJLBX+wDUr/+r6bqGe2Z+G3hpRfFOYH9Z3g/c3FZ+b7Y8BmyKiEuH1VhpmaEtra3fOfdLMvMEQLm+uJRvAV5oq7dUys4SEXsi4nBEHD516lSfzdCsGHUYT+tpV6VJGfYHqtGhLDtVzMx9mbmQmQtzc3NDboamRa/Hja+cnqnxJ+vUTJN6vvYb7i8uT7eU65OlfAnY1lZvK3C8/+ap6Qadk1/tPqVxmORzrd9wPwDsKsu7gAfaym8rR83sAE4vT99Iksanl0MhvwT8G/B7EbEUEbuBO4DrIuIocF1ZB/gGcAxYBD4LfHAkrVa11pqa6XcU5EhdTbSxW4XMfN8qm67tUDeB2wdtlMZnfu9DPH/HTWO9z2H/gpLhLZ3Nb6hqYIarNH0Md63K0JZml+GuofF0udL0MNzVVS8hbZBL08Vw128YNKQ9p4s0HQx3GcRShQz3KTOpoO321f/Vvik6rikb34Ck9THc1RfDVppuhrskVchwr0S3KZNOy2sdurjW2Rslrd+4XzuGe2U8WkUSGO5TYdg/NTesut3uwzcQqTeTeK0Y7n3qNKUxrGPEB60zjNtImm2G+wwyrCV1Y7j3YZinrB3FaN/wl2S4D9nK+ei1fht03CHsj0hLkzXO153hvg7jmFMf9v0Z4tJ0Gddr0nBfoZ9jw6clQKelHZImz3DvoNeQHGWY9vLGIkmrMdzHZNBDGA12SevRqHAf9g8zD/J40/DfgaTJGMfrulHh3q9JfBA6jvuQNDmjfg03OtzXmvYwPCXNskaHuyTVqvpwX21E7py3pEkbZb5UH+7LhvHhpiTNiurCvZeRuYcYSqpddeHezuCW1FQbR3GnEXE98HfABuBzmXnHKB5nWbcjXQx5SU0z9JF7RGwA/h64AbgceF9EXD7sx1lmcEvS2UYxLXMVsJiZxzLzf4EvAztH8DiSpFWMYlpmC/BC2/oS8AcrK0XEHmBPWf15RPwE+PEI2jOLLsK+WGZfnPFaX8SdE27JGsbQtqqeEwP21++utmEU4R4dyvKsgsx9wL7XbhRxODMXRtCemWNfnGFfnGFftNgPvRnFtMwSsK1tfStwfASPI0laxSjC/bvA9oi4LCLOBW4FDozgcSRJqxj6tExmvhoRfwF8i9ahkJ/PzGd7uOm+7lUaw744w744w75osR96EJlnTYdLkmZc1d9QlaSmMtwlqUITD/eIuD4ifhgRixGxd9LtGYWI+HxEnIyIZ9rKLoiIgxFxtFxvLuUREXeX/ng6Iq5su82uUv9oROyaxL4MKiK2RcQjEXEkIp6NiA+V8sb1R0ScHxHfiYjvl774RCm/LCIeL/t1XzkwgYg4r6wvlu3zbff1sVL+w4h412T2aDARsSEinoyIB8t6I/thaDJzYhdaH7g+B7wZOBf4PnD5JNs0ov38I+BK4Jm2sr8F9pblvcCdZflG4Ju0vi+wA3i8lF8AHCvXm8vy5knvWx99cSlwZVl+A/AjWqepaFx/lH16fVk+B3i87ONXgFtL+WeAPy/LHwQ+U5ZvBe4ry5eX1855wGXlNbVh0vvXR398BPgi8GBZb2Q/DOsy6ZF7I05VkJnfBl5aUbwT2F+W9wM3t5Xfmy2PAZsi4lLgXcDBzHwpM18GDgLXj771w5WZJzLze2X5Z8ARWt9qblx/lH36eVk9p1wSuAa4v5Sv7IvlProfuDYiopR/OTN/mZn/CSzSem3NjIjYCtwEfK6sBw3sh2GadLh3OlXBlgm1ZdwuycwT0Ao84OJSvlqfVNdX5d/pt9MasTayP8pUxFPASVpvUM8Br2Tmq6VK+369ts9l+2ngQuroi08BHwV+XdYvpJn9MDSTDveeTlXQMKv1SVV9FRGvB74KfDgzf7pW1Q5l1fRHZv4qM6+g9U3uq4C3dqpWrqvsi4h4N3AyM59oL+5Qtep+GLZJh3uTT1XwYpleoFyfLOWr9Uk1fRUR59AK9i9k5tdKcWP7AyAzXwEepTXnvikilr9g2L5fr+1z2f5GWtN9s94XVwPviYjnaU3NXkNrJN+0fhiqSYd7k09VcABYPsJjF/BAW/lt5SiRHcDpMk3xLeCdEbG5HEnyzlI2U8rc6D3Akcz8ZNumxvVHRMxFxKay/DrgHbQ+g3gEuKVUW9kXy310C/Bwtj5JPADcWo4iuQzYDnxnPHsxuMz8WGZuzcx5WhnwcGa+n4b1w9BN+hNdWkdD/IjWXOPHJ92eEe3jl4ATwP/RGl3spjVHeAg4Wq4vKHWD1o+dPAf8AFhou58/o/Uh0SLwgUnvV5998Ye0/lV+GniqXG5sYn8Avw88WfriGeCvS/mbaYXSIvBPwHml/Pyyvli2v7ntvj5e+uiHwA2T3rcB+uSPOXO0TGP7YRgXTz8gSRWa9LSMJGkEDHdJqpDhLkkVMtwlqUKGuyRVyHCXpAoZ7pJUof8HCol+hyjh9nkAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.hist(tot_sam, bins = int((x_max-x_min)/7.))\n",
"\n",
"plt.show()\n",
"# _ = np.where((total_samp >= x_min) & (total_samp <= (jpsi_mass - 50.)))\n",
"\n",
"# total_samp[_]"
]
},
{
"cell_type": "code",
"execution_count": 42,
"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
}
