{
"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_const1 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_cut(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",
" #Cut out jpsi and psi2s\n",
" \n",
" tot = tf.where(tf.math.logical_or(x < ztf.constant(jpsi_mass-60.), x > ztf.constant(jpsi_mass+70.)), tot, 0.0*tot)\n",
" \n",
" tot = tf.where(tf.math.logical_or(x < ztf.constant(psi2s_mass-50.), x > ztf.constant(psi2s_mass+50.)), tot, 0.0*tot)\n",
" \n",
" return tot\n",
" \n",
"class total_pdf_full(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",
" #Cut out jpsi and psi2s\n",
" \n",
"# tot = tf.where(tf.math.logical_or(x < ztf.constant(jpsi_mass-60.), x > ztf.constant(jpsi_mass+70.)), tot, 0.0*tot)\n",
" \n",
"# tot = tf.where(tf.math.logical_or(x < ztf.constant(psi2s_mass-50.), x > ztf.constant(psi2s_mass+50.)), tot, 0.0*tot)\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), lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"rho_s = zfit.Parameter(\"rho_s\", ztf.constant(rho_scale), 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), lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"omega_s = zfit.Parameter(\"omega_s\", ztf.constant(omega_scale), 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), lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"phi_s = zfit.Parameter(\"phi_s\", ztf.constant(phi_scale), 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.7, upper_limit=1.7)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Load data"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"x_min = 2*pdg['muon_M']\n",
"x_max = (pdg[\"Bplus_M\"]-pdg[\"Ks_M\"]-0.1)\n",
"\n",
"# # Full spectrum\n",
"\n",
"obs_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, jpsi_mass - 60.))\n",
"obs2 = zfit.Space('q', limits = (jpsi_mass + 70., psi2s_mass - 50.))\n",
"obs3 = zfit.Space('q', limits = (psi2s_mass + 50., x_max))\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_cut(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_full(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/Sum:0' shape=() dtype=float64>"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"total_f_fit.normalization(obs_fit)"
]
},
{
"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",
"\n",
"\n",
"test_q = np.linspace(x_min, x_max, int(2e6))\n",
"\n",
"probs = total_f_fit.pdf(test_q, norm_range=False)\n",
"\n",
"calcs_test = zfit.run(probs)\n",
"\n",
"Ctt.set_value(0.5)\n",
"\n",
"calcs_test1 = zfit.run(probs)\n",
"\n",
"# Ctt.set_value(4.9)\n",
"\n",
"# calcs_test2 = 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:14: UserWarning: Creating legend with loc=\"best\" can be slow with large amounts of data.\n",
" \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+17YcXAAAgAElEQVR4nOy9eZycVZX//z61955OL9k6SSd0h5CEJCQhIAFkCRANEtSoQXBQEfw5qKMzzogzo8ww4EsYBRdEvggqbgSMOIQdZFHWhGyEJGRpsnb2pfel1vP743mqu9JdW3e6u6q77/v1alJ1n/uce6tInk+fc889V1QVg8FgMBiyAUemJ2AwGAwGQxQjSgaDwWDIGowoGQwGgyFrMKJkMBgMhqzBiJLBYDAYsgZXpieQbZSWlmplZWWmp2EwZD2hAxvxuwrIK590UntLwzHyWvYRGHk6Hl9uhmZnETnwLu3uEeSWTaSx4QSFLXsIj6zG6cvP6LyGImvXrj2mqmWnaseIUhcqKytZs2ZNpqdhMGQ99f9VwZaRl3He1399UvvqZ37D/NX/xN5P/5oJ087J0Owsmm8dxdaxH2fel+/npacf49J3buTEp3/ByGkXZ3ReQxER2dMXdkz4zmAw9AoHEXB0f4SIwwmARsIDPaVuCAqI9dqeVySc+XkZEmNEyWAw9ApBUYknSlYAJhIODvSUuiEAYolSVECzQSwNiTGiZDAYeoWDCBLnESLO7PBIVNXylGzhzCYPzpAYs6aUBsFgkNraWtrb2zM9FUMKfD4fFRUVuN3uTE9lyOMgEtdTctiekmZclMARE75DoqIUydykDCkxopQGtbW1FBQUUFlZiURDAYasQ1U5fvw4tbW1TJo0KfUNhlPCqREkzpoSTjt8F8ls+E6j/7X/zTqc1lwjkVDG5mRIjQnfpUF7ezslJSVGkLIcEaGkpMR4tAOE5Sk5u7dHw2QZ95T0pEQHxPbgjKeU1RhRShMjSIMD8/9p4HBJBHUkEaUMr91E9OREh6hXp2rWlLIZI0oGg6HnRL2NOJ6SOK31vExn3ylqpa3b615RscQkOmQ1RpSGEK+++ipXXnklAH6/n4ULFzJ79mweffTRbn2/8Y1v8Pe//x2wEjluueUWqqurmTFjBvPnz+fZZ58F4Pvf/37HPfX19dx33319Mle/389nPvMZqqqqOOecc9i9e3fcfs899xynn346VVVV/OAHP+hoX7ZsGTt27OiTuRh6gb0uEy98F82+y7SnpB2ekt0g0axAs6aUzRhRGqKsX7+eYDDIhg0b+MxnPnPStRMnTvD2229z4YUXAvDd736XgwcPsmnTJjZt2sSTTz5JU1MT0H+i9NBDD1FcXExNTQ3f/OY3+fa3v92tTzgc5uabb+bZZ59ly5YtPPLII2zZsgWAr3zlK9x11119MhdDL4iGwOKE77Il9VoVHKJEH3NRscQcbJrVpCVKIrJIRLaJSI2I3BLnuldEHrWvrxKRyphr37Hbt4nIFalsisgk28YO26Yn2RgiUiIir4hIs4jcm2D+K0VkU3pfSfaxe/dupk6dyvXXX8/MmTNZunQpra2tgOVJTJ06lfPPP5/HH38cgCNHjnDdddexYcMGZs+ezQcffHCSvRUrVrBo0SIAWltb+eUvf8nPfvYzvF4vAKNGjeLTn/40t9xyC21tbcyePZtrr72WW265hQ8++IDZs2fzr//6r6f0mZ544gmuv/56AJYuXcpLL71E11OQV69eTVVVFZMnT8bj8bBs2TKeeOIJAC644AL++te/EgqZ33ozQYe3EcdTctrhO81wlptqNMTYZU3JZN9lNSlTwkXECfwcuAyoBd4RkZWquiWm2w1AnapWicgy4E7gMyIyDVgGTAfGAn8VkSn2PYls3gnco6rLReR+2/YvEo0BtAPfBWbYP13n/wmguUffShL++8nNbDnQ2FfmAJg2tpBbPzY9aZ9t27bx0EMPsWDBAr74xS9y33338dWvfpUbb7yRl19+maqqqg6PqLy8nAcffJAf/vCHPPXUU91svfHGGyxduhSAmpoaJkyYQGFhYbd+P/jBD7j33nvZsGEDYInjpk2bOt535YILLujwsGL54Q9/yMKFC09q279/P+PHjwfA5XJRVFTE8ePHKS0tjdsHoKKiglWrVgHgcDioqqri3XffZe7cuYm/OEO/EAmHrN9o4yU6ZMnm2Uj0l5yoKEl2ZAUakpOOpzQfqFHVnaoaAJYDS7r0WQI8bL9eAVwqVhrUEmC5qvpVdRdQY9uLa9O+5xLbBrbNq5ONoaotqvo6ljidhIjkA/8M3J7G58xqxo8fz4IFCwC47rrreP3119m6dSuTJk2iuroaEeG6665Ly9bBgwcpKzvlYr7deO2119iwYUO3n66CBHTziqB75lyqPuXl5Rw4cKAPZm7oKeEknlK0zBCZ9pTsZAyJ1r6z9091eFCGrCSdzbPjgH0x72uBrqV/O/qoakhEGoASu/3tLveOs1/Hs1kC1KtqKE7/RGMcSzL3/wF+BLQm+4AichNwE8CECROSdU3p0fQXXR/Y0vHbX89ToHNycjr28lRVVbF3716ampooKCg4pTn2xFOqqKhg3759VFRUEAqFaGhoYOTIkXH7RKmtrWXs2LEd79vb28nJyTmlORt6R3S9KG5KuDM7KjpEPSWNbp6NbvQ12XdZTTqeUrynXtdfYRP16av2dOfROSGR2UCVqv4lUZ8OI6oPqOo8VZ3XHx5EX7B3717eeustAB555BHOP/98pk6dyq5duzrWjB555JG0bJ1xxhnU1NQAkJubyw033MDXv/51AoEAYHlSv//97wFwu90Eg1Zqb0FBQVzRidITT+mqq67i4Yctx3fFihVccskl3QT27LPPZseOHezatYtAIMDy5cu56qqrOq5v376d6dMz80vCcCfZmlKHKGV6P1AkGr47ufZdxGyezWrSEaVaYHzM+wqga8yko4+IuIAi4ESSexO1HwNG2Da6jpVojER8CJgrIruB14EpIvJq0k+axZxxxhk8/PDDzJw5kxMnTvCVr3wFn8/HAw88wOLFizn//POZOHFiWrYWL17Mq6++2vH+9ttvp6ysjGnTpjFjxgyuvvrqjvDeTTfdxMyZM7n22mspKSlhwYIFzJgx45QTHW644QaOHz9OVVUVd999d0e694EDB/joRz8KWGtN9957L1dccQVnnHEGn/70pztE6PDhw+Tk5DBmzJhTmoehd4RCluBIkjWlTCcUREVRuqwpkWmxNCRHVZP+YIX4dgKTAA/wLjC9S5+bgfvt18uAx+zX0+3+Xvv+nYAzmU3gT8Ay+/X9wD8mGyNmDp8H7k3wGSqBTak+q6oyd+5c7cqWLVu6tQ0ku3bt0unTp/epzQULFmhdXV2f2hxI7r77bn3wwQfjXsv0/6/hQP3+Haq3Fuprj97d7dq+PTtVby3UtSt+mIGZdXK8vkH11kJd87v/UFXVNZu2qt5aqDVPdZ+z4dQB1mgaz9hUPynXlNRav/kq8LwtKL9S1c0icps9iZXAQ8DvRKQGy3tZZt+7WUQeA7YAIeBmtX99iWfTHvLbwHIRuR1Yb9sm0Ri2rd1AIeARkauBy/Xk7EBDF370ox+xd+9eRowYkemp9IoRI0bwuc99LtPTGLaE7PBdR1JDDFFPKdNrNx2JDh0FWV3RC5makiEN0qoSrqrPAM90aftezOt24FMJ7r0DuCMdm3b7TqzsvK7tycaoTDH/3cRJFx8sVFZWsmlT326zOueczB5Tfap84QtfyPQUhjXRRId44TtnlqwpdSY6mPOUBhOmooPBYOgxHSnhzjiekiuafZfhTaoJNs9m2oMzJMeIksFg6DEaiobvElcJz/jD3/aUovuUHGaf0qDAiJLBYOgx4XBiUXJ2rN1kOHzXJSU8a8TSkBQjSgaDocdESwhFj6mIxemy2zKdEk6i2ndGlLIZI0pDiKF4dEVlZSVnnnkms2fPZt68eR3t3/rWt3j55Zf7ZC6GnhM9KymupxRdU8r0IX8m+25QYkRpiDIUjq6I8sorr7BhwwbWrFnT0fa1r33tpPOVDANL1FNyZPGaUufaUbTMUHbMy5AcI0qDgOF6dEUyJk6cyPHjxzl06NApzcPQOyJ2aK7jjKIYnA4hqE4k05UTuq0pOYioZDxV3ZCctPYpGWJ49hY49F7f2hx9Jnwk+W/9w/HoCrBCL5dffjkiwpe//GVuuummjmtz5szhjTfe4JOf/GTC783QP0RDcw5HnDUlEYI4Mu6RdJ6nZIuSQBgHYsJ3WY0RpUFC16MrfvrTn7Jw4cKOoyui7Q888EBKW/15dEW6xPOK4lU8f+ONNxg7dixHjhzhsssuY+rUqR1hR3N0ReboWFOK4yk5HEIYR8az77qLkhBBMr7WZUiOEaWeksKj6S+G49EVQMdRFeXl5Xz84x9n9erVHaJkjq7IHNFjKeIlOoDlkWQ8+y5aDdzOunM6hAgOk+iQ5Zg1pUHCcDy6oqWlpWO8lpYWXnjhBWbM6KwWtX379pPeGwaOsC04jji174Ds8JQ6PCK7zJBABDGilOUYURokDMejKw4fPsz555/PrFmzmD9/PosXL+5I0AgGg9TU1JyUJm4YQMJ2urUrvihFcCAZPrcoXvgunAVrXYbkmPDdIMHhcHD//fd3a1+0aBFbt27t1n7RRRdx0UUXxbV1wQUX8J3vfIf6+npGjBiBx+Phrrvu4q677urW98477+TOO+/seP/HP/6x9x8iBp/Px5/+9Kdu7WPHjuWZZ6w6vZMnT+bdd9+Ne/9TTz3F0qVLcSV4KBr6l+iaUryUcIAwTtAsqX1nwneDCuMpDVOiR1cMVkKhEP/yL/+S6WkMW6Jp1Y44BVkBIpIFnlLH5tmu4TvjKWUz5tfMQYA5uqI7n/pU3FNMDANENNHBGSf7DixPSTLsKXUkOtii5DThu0GB8ZTSpCcbOw2Zw/x/GhiiR51LnH1KAKFsECU92VOyUsJN+C7bMaKUBj6fj+PHj5sHXpajqhw/fhyfz5fpqQx5OjwlVwJPSVxItqSER2vf2Z6SGE8pqzHhuzSoqKigtraWo0ePZnoqhhT4fD4qKioyPY0hT6enFP8REsKFIxIcyCl1I7rupR1lhiCsmffgDMkxopQGbrebSZMmZXoaBkPWEH3gOxNkP4bFhSPDnlLHIX8xnlKQLKjJZ0iKCd8ZDIYeEz3qPHH2nQvRDHtKHQVZrRCjFb5zZjysaEhOWqIkIotEZJuI1IjILXGue0XkUfv6KhGpjLn2Hbt9m4hckcqmiEyybeywbXqSjSEiJSLyiog0i8i9MXZyReRpEdkqIptFxJxzYDD0FZEU2XdZ4SlFSyFFj66wEzAyPS9DUlKKkog4gZ8DHwGmAdeIyLQu3W4A6lS1CrgHuNO+dxqwDJgOLALuExFnCpt3AveoajVQZ9tOOAbQDnwX+Fac6f9QVacCZwELROQjqT6vwWBITUeV8CSekiNLsu+iKeEuh4MQDiNKWU46ntJ8oEZVd6pqAFgOLOnSZwnwsP16BXCpWIHcJcByVfWr6i6gxrYX16Z9zyW2DWybVycbQ1VbVPV1LHHqQFVbVfUV+3UAWAeYFXCDoS/o8JQSryk5My1KHdl3ljfndAghXCbRIctJR5TGAfti3tfabXH7qGoIaABKktybqL0EqLdtdB0r0RgpEZERwMeAlxJcv0lE1ojIGpNhZzCkRu3MOocz/j6liCPzogQnH4cOmDWlQUA6ohTvbISuG3YS9emr9nTn0Q0RcQGPAD9V1Z3x+qjqA6o6T1Xn9cc5QwbDkCPqKbkTiFI2hO+i+5Gkc90rLCb7LttJR5RqgfEx7yuArierdfSxRaAIOJHk3kTtx4ARto2uYyUaIxUPADtU9cdp9DUYDGkgdkFWp8sT97o63Jn3lDqy7zp/n43gRDK8f8qQnHRE6R2g2s6K82AlLqzs0mclcL39einwslrlD1YCy+zMuUlANbA6kU37nldsG9g2n0gxRkJE5HYs8fpGGp/TYDCkSyREWCVh9l1EXDjJrChFOo5s75xjSFw4op5SOAQ1f+3Yz2TIDlKKkr1+81XgeeB94DFV3Swit4nIVXa3h4ASEakB/hm4xb53M/AYsAV4DrhZVcOJbNq2vg38s22rxLadcAwAEdkN3A18XkRqRWSaiFQA/4GV3bdORDaIyJd6/hUZDIauSCRICBdOR/yTj9XhwpnpQ/6ite8cnY+5SGxNvlX3w+8/CTVxl5oNGSKtig6q+gzwTJe278W8bgfilm1W1TuAO9KxabfvxMrO69qebIzKBFPv+VnhBoMhJRIJEcRJrsT/JxZxuHFlfE3p5IKsYK0pdeyfOvq+9WfDvq63GjKIqehgMBh6jOUpOXEk9JTcGQ/fdZz5FBO+i4izM3znzrP+DLV3vdWQQYwoGQyGnhMJEyL+ehIADhcuMhu+i4RPPnkWumQFRovJmhTxrMKIksFg6DGiQevI8wSow4Urw55S9NwkRzdRssUy6kEZUcoqjCgZDIaeEw4RksRL0upwW55SBjPbIpGTq4RDAk8pbEQpmzCiZDAYeoxEkntKOD040IwePa4dBVlP3jzb4SlFxSrsH+ipGZJgRMlgMPQYiYQIJUve7VivydxG1Wj2XWyigzpi9k/ZG4AJtg301AxJMKJkSMiRA3t5a4UphGHojkNDhJOE77Br4mk4MEAz6o6mWlOKriUZUcoqjCgZElL/0Cf40KZbOXpgT6anYsgyJBIiLMmy7yxRCgcz7ynFbp5VcXaWP4oKpkkJzyqMKBkSUhS2SguGw6ZWmOFkRENEkqwpie0phUKZ9JS6rylFHC6c0VR1E77LSowoGQyGHuOIhIgkCd+Jfc5SKJhBUYp0D98hLpxErKzAaPjOeEpZhRElg8HQY1KvKVnVw8OBzGW2dYiSnOwpWS9CnetdRpSyCiNKBoOhxzg01PmAj0N2hO+6V3ToOFspHOTgiSYAmpqbB3hmhmQYUTKkxpT2N3TBoSnCdy470SGUwfVIW5SczpjsOzsBg0iI442WGAXbWwZ8aobEGFEyGAw9xplKlKLhuyxYU5KY8J3EhO/c9nHpjrAJ32UTRpQMqUlwPIFh+OLQMJokJVzsE2nDGVyvibtPqaO0UBAX9um5pqJDVmFEyZAaE74zdCHVmpLD7QMgksFEB6Ip4THhu9jK4NFDCI2nlF0YUTIYDD3GqWE0SfguKkrhQAY9pY7su1hRsteUwoGOKubOSOZCjIbuGFEypMaE7wxdcBFCow/4ODg9UVHK3MZUtT38kw75c1hhRcJBHBETvstGjCgZUmPCd4YuODWEJgnfuTy5AIQzWS1Bu5cZithrXYTaEbvckFv95u94FmFEyWAw9Bgn4aSekttrrykFM+iF2MdmOGNr33V4SgEcsYf7mQ20WUNaoiQii0Rkm4jUiMgtca57ReRR+/oqEamMufYdu32biFyRyqaITLJt7LBtepKNISIlIvKKiDSLyL1d5jVXRN6z7/mpiIlD9QrztRm64CTcmTQQB5cnB4BIJj2lSDTRoVM8I06v9SLkx6Exe6hM/busIaUoiZXk/3PgI8A04BoRmdal2w1AnapWAfcAd9r3TgOWAdOBRcB9IuJMYfNO4B5VrQbqbNsJxwDage8C34oz/V8ANwHV9s+iVJ/XEAcT2jB0wZ1iTcntjYpS5rPvHM6Y1HVnZ/juZE/JrCtlC+l4SvOBGlXdqaoBYDmwpEufJcDD9usVwKW2V7IEWK6qflXdBdTY9uLatO+5xLaBbfPqZGOoaouqvo4lTh2IyBigUFXfUmvF87cxtgxpoBgPyRCHSAQPISLRB3wcouE7DWYwLBaJilKMeLqseREO4NQQAbUFK2Q8pWwhHVEaB+yLeV9rt8Xto6ohoAEoSXJvovYSoN620XWsRGMkm3dtinkDICI3icgaEVlz9OjRJCaHF4LxkAxxsAuZajQUFgeP10p00Eyu1UTXlGI8JXFZc44E/Tg0RDOWR0cmxdNwEumIUrxfl7s+rRL16av2dOeRzpy6N6o+oKrzVHVeWVlZEpPDE+MvGU7CTqHWJJ6S1+smoM7Mekodm2c7176i+6dCwXZcGqJJLfE0nlL2kI4o1QLjY95XAAcS9RERF1AEnEhyb6L2Y8AI20bXsRKNkWzeFSnmbUhCNHynxmMyxBBdJ0omSh6nAz+eztNdM0A05VtiEjKcHstTCgfacGA8pWwkHVF6B6i2s+I8WIkLK7v0WQlcb79eCrxsr+OsBJbZmXOTsJINVieyad/zim0D2+YTKcaIi6oeBJpE5Fx7reofYmwZ0kCSOqmG4UowuiHWlTh853U58ONGsiB8R0yNPoers9KEKzZ8Z1LCs4Ykp3RZqGpIRL4KPA84gV+p6mYRuQ1Yo6orgYeA34lIDZb3ssy+d7OIPAZsAULAzWqfURzPpj3kt4HlInI7sN62TaIxbFu7gULAIyJXA5er6hbgK8BvgBzgWfvH0EOMJBliCQX8eKEzky0OIkIQF5LBagnScZ5SpyhFK01EQn7chGjSARSl5qPw3mNw9o3gSvzdDXdSihKAqj4DPNOl7Xsxr9uBTyW49w7gjnRs2u07sbLzurYnG6MyQfsaYEa8awaDoXeE7Hp2jiSeEkAAD5LBVGvREGEVnDH77JzRQrF+6wyllo7w3QCsKb35U+uncBxMN4nAiTAVHQwpMStKhliCfluUbK8jYT/x4IhktqJDpMsjzmWvKUXa7VNnB9JTOrzJ+vPQxv4faxBjRMmQko5jpQ0GOteUUnlKQXHjyGiiQ5gQJ5/55Ha58KsLApan1IRdoy/Q2v8Tajps/Vm/L3m/YY4RJYPB0CNCtig5U3pKXhwZPBZCNExETn7EeV0OArjBbx2FHvWUQgNwxEZ74xEAGo/s6fexBjNGlAwpMVWGDLGE7YP7UnlKIYcHZwbDdxoJE+7iKXlcDgK4IGCJUjT7LuTvZ09JFVd7HQCRhtoUnYc3RpQMCTEVHQzxiHpKLru+XSLC4saVaU+pyyPO7bRT1YNW+C7k9BFSR/+Lkr+p41BBnz/Z9kqDESVDaoyrZIghbG+ejWayJezn8ODKskQHj8tBQN04gpan5HJ58OMm0t+HEbYeB2BnZDQ+bTObdZNgRMmQkI6KDkaUDDFERSladDURAWcO3kjmHr6iYcLSJXxne0pOO3zncntpxUfEXmPqN1ot7+gDtctv2iJl6I4RJUNCTPjOEI9omSFXCk8p5MzFo5mt6BDPU2rFhyvYCIDb46FRc6G9sX/nYovQBzoWAG091r/jDWKMKBnSwKSEGzqJ2Ht6PCk8pbArB18GRclaU+ruKbWqF3fAEiGX20szOYi/qV/nEmmrB+CAYwwALXXmNIJEGFEyGAw9QqOekid5okPElYeXQGcNuoEm0j0lPOopRbMCPV4vTZqDBPrXU/K3WJl3wRGVALTUH+7X8QYzRpQMKdGICeMZOlG7dJDHl9xTwm0fC2FvVB1o4npKLgctdKayezxemsjFGehfT8nfbHlK7tLTrPcNRpQSYUTJkBJzdIUhlqgoeVOE79STZ73ImChFiHRNdHA5aNXOeXu9Ppo1B2ewfxMdgi31BNVJ8agJRFQINJk1pUQYUTIYDD0j1E67uvG4k9dzdngtUQq293NmWwJEQ6h03acktBArSpan5A717xxDbY00kUNleSH15BFpNqKUCCNKhtSY2neGGBzBNtrw4nEmf3yINx+A9tZ+zmxLNL5GumffOa01pSi+DlFqgUj//T2PtDXQpLlMKs3nhBYiJiU8IUaUDAZDj5BQG+14cKUQJaftKflb+3e9JhEODaNdwnciQtDRmaDh9uVbiQ4o9Oe6kr+RJnIpzfdQL0W42k1Vh0QYUTKkxGyeNcTiDLfTTvK6dwBOXwEAgQyJksQRJYCQq1OUxJNLm8Ne++rHtHCHv5EmzaUox02zawS+gBGlRBhRMqTEaJIhFke4Db+kFiW3z3rYhzK2phTulugAEHLldvbx5OJ3WWFG2hv6bS7OYBPN5JDncdHmLiYvVNdvYw12jCgZDIYe4Qy1pydKOYUAhNoyI0pOQoTF3a09HBUhAE8ura5i63VL/yUfuIPNtDvzcDiEoK+E/Ehj5vZvZTlGlAypMYkOhhic4TaCjhR7lABPrvXwD/Z3XbkEuDRExNFdlFo9IzteO1w+Qj77fUv/VVnwhFsI2GIYySnBgXbUwzOcjBElg8HQI9yR9rREyZtreUqR9szsU3JqkIh0T1v3e0s7XrucDsI59vv+8pRU8UVaCLmtNTZnfhkAkWZTaigeaYmSiCwSkW0iUiMit8S57hWRR+3rq0SkMubad+z2bSJyRSqbIjLJtrHDtuk5hTG+KSKbRWSTiDwiIqn/JRm6YxaVDDG4In7CaYhSbm4BYRXwZyYl3Klhwo7uohTwlXS89rqdOHJHEsbRf55SoAUnEcIeS5TcheUANJ042D/jDXJSipKIOIGfAx8BpgHXiMi0Lt1uAOpUtQq4B7jTvncasAyYDiwC7hMRZwqbdwL3qGo1UGfb7s0Y44CvA/NUdQbgtPsZeoip6GCIxRNpI+xKLUo5XhdN5CL+/ksgSIaLIBGHp1t7NCtwv5bgdTkoyPVQTyH0V+VuW5TVa3mOvuJRALScMKWG4pGOpzQfqFHVnaoaAJYDS7r0WQI8bL9eAVwqImK3L1dVv6ruAmpse3Ft2vdcYtvAtnl1L8cAcAE5IuICcoEDaXxeg8GQBK/6CTlzU/bL8zpp1FwkQ56SS0NonPBdrtvJRf4f8XH/bXhdDgp9bo5rYf+F7+xjMcRriWF+8WgA2hoO9c94g5x0RGkcsC/mfa3dFrePqoaABqAkyb2J2kuAettG17F6NIaq7gd+COwFDgINqvpCvA8oIjeJyBoRWXP0qInzdsXsUzLE4lE/mo6n5HbSSD6uDHlKbuInOuR4nOzWMRyhGK/LSVGOmyORArT5SP9MxN7/5MwpAmBEyWgiKgQb+2m8QU46oiRx2ro+pRL16av2Ho8hIsVYXtQkYCyQJyLXxemLqj6gqvNUdV5ZWVm8LgaDASASxksQdSc/tgKs6gktjjxcwcxsnnUSQhOIUhSv22+Q3FwAACAASURBVEFhjpsjFKON/RNICdjHVrhyRwBQVphLHfmoSXSISzqiVAuMj3lfQfcwWEcfO1RWBJxIcm+i9mPACNtG17F6OsZCYJeqHlXVIPA4cF4an9fQFZMSbogSbLP+dKUO3wG0OQvwBDMTvnMnSAnPixUll4OiHDe1Woo0HYRwqFv/U6W10apz58639kMV5rg4QZGpf5eAdETpHaDazorzYCULrOzSZyVwvf16KfCyWjGflcAyO3NuElANrE5k077nFdsGts0nejnGXuBcEcm1154uBd5P72sxxGKid4YoGmwFQFIc8Bcl4CrAF86Mp+QmBM7uiQ5FuZ1tXpeTQp+L/VqGaBga9/f5PPz2MRXeAiv1XERodIzA4zeiFI/kteex1m9E5KvA81gZbL9S1c0ichuwRlVXAg8BvxORGizvZZl972YReQzYAoSAm1U1DBDPpj3kt4HlInI7sN62TS/GWCUiK4B1dvt64IHeflEGgwECrY14AbUX7VP2dxeSm4HNs6qKizA4u3tKI2NFye2gtMBLrdp7lRr2QfHEPp1LqMkSn9yizv1Rze4Sxge29ek4Q4WUogSgqs8Az3Rp+17M63bgUwnuvQO4Ix2bdvtOOrPnYtt7M8atwK3x7jGkj0l0METxN9fjBcSTniiFPYV48UPID67UpYn6ikAwiFciSBxRKs7rbMvzuCgv8FKr9lpy/d4+n0u45TiNmkNhXmfIs8U3iuLGN6zjMhymhkEs5tswJETM/iRDFwIt1rHeDruuXSrUZ2Wc9Wex03gEAtbpuBInfDcyr7PN6RBK870c1BIU6RdR0tYT1Gs+I3I7xTCQN9YKL/bX3qhBjBElQ2pMooPBJmAf2JeuKOGzMs4GWpSCfluUXHFEKffkNp/bSU5OLnWeMXC070NqzvYT1FFAUU6nKDlHVADgP973IjjYMaJkSIkJ3xmihFotT8lt77lJhSPXyjgLNA3snpxAMLEoleZbYcQCX+fqRVmBl1rXBDi6tc/n4vLXUa/5FPg6RSmn1Fq3qj+0u8/HG+yktaZkMBgMAEHbU/LlpylKBVadt/a6Q3gm9du0uhGyw3eOOOE7h0P4443nnBTGKy/wsqNuPDOP/R+Eg3ETJHqLJ9BAs3MyTkfnlsqi0ZUAtBzd3WfjDBWMp2RIA+MpGSxCbVYYLqegOK3+zsIxAPgHuKRO0G/tpxJ3/OSK804rZerozhDk6EIfmwJjIRKE4x/03URUKQgepdFVclLz6NHjaFc3gRP7Etw4fDGiZDAY0ibS1kBAnRTk5afuDOQWjyKiQqB+YIuPhvw92081oSSXVa1WTToOvdfzAV+/x/rpGupuq8OtQVq85Sc1jx6Rw0EtQRpqez7WEMeIkiElZk3J0IHfOta7ICe98FapXVIn3DSwohQOWKLk8KRXeWJiSS7bIhVE3HlQu7png9Wuhb/+l/VT89LJ1+zSRcHc0Sc1u50ODjrHkteyp2djDQOMKBkMhvTxN9GsOeR701uOLsv3ckyLoL+KnSYgYntKzjRq9AFMLMkjjJP6kTNh79s9G2zfKvuFwNpfn3wtWk/PDmPGciJnPKX+faZkSheMKBlSY/7RGGwcwWZaJReXM71HR3Geh2NahKttYIuPdnhK3vQ8pcqSPAD25J4JhzdBT6pQHNsGOSPhQzfD9uegpbN8UKTeWjNyjhjf7bZg0WR82g5N5giLWIwoGVJiwneGKO5AIy2O9NaTwApTNTiL8Q1wnbdIwEp0cKYpSsW5bgp9LtbLGda+vD1vpD/Y0W1QNhVmLYNICDY/3nGp/dD7tKqXnJKKbre5y6cA0Hyg79PQBzNGlAwGQ9rkhOppcaa5cdam2VNKfvDYgHrcantKrjRFSUSYNraQpxsngzsPtj+f5kCKHt3KW00lXPdUK+GyafDu8o7L4cPb+EDHUFbYfR4jxp8BwPG9W9Iba5hgRMmQBsZTMljkhhpoc43o0T0tOWPxaABaBi6EF7GP2HB789K+Z8bYIjYdaiMy+WJLlNIR0ZZjSFsdLx4Zwes1x3gzbyHsXwPHakAVz7EtbNcKxo/svrZVMbGKVvXiP2gOL4jFiJIhIab2neEkVMmPNOJ3p7dxNoo/zw5d9UNduUSE7UQHX276ojR9XCH+UITDYy6Gxlo4sC71TXYFiB06jrkTi7lj3wxUHLBxORyvwes/xjuRqUws6T6P8SX5bGc83mObu10bzhhRMqREI6b2nQFob8BJhJAvvY2zUXTEBOvPuoFLf46uKfly0helM8dZHuCb7nPB6T0pDJeQY1atvLrcSXz90mq2tuRzuPxCWPUA/O0uIgibfHPjZis6HUKtbwqlzdtMMlEMRpQMCVH7pHk1HpMBIHpSas7IHt2WV27VF2o7uquvZ5SQ6JpSbl56R2wAnFaWR1mBl7/tDcIZV8J7f7KO3EjG0W20SQ6F5RO5oKqUSaV5/E/ws1bCw3uP8aZ3Qcfnj4e/ZDp52oLW7U57nkMdI0qGhJjwnSGWkJ3q7Mwv69F95aWlnNB82o4MnCg5gk20qBen05m6s42IcN5pJbz5wXF09rXQVgebHk96jx7dRk1kLFWjCnA4hOvOncjTB/LZsWQlwcU/4ebmLzF7QuI1uNzKuQAc2/FO2vMc6hhRMqTGhBYMQHOdtQHWU1CSoufJVBTnsE/L0RM7+2NacXEGmmmR9DLvYllQVcqxZj9bcuZA2Rnw5k+T/v2PHNnKtsg4qsutNPmlcyvIcTv55VY3bxR+lIawh7MnJvYsx58+l4A6afighxt2hzBGlAwGQ1q01VvZc76innlKFcU51Og4chp29Me04uIONdEs6e+ninLp1HKcDuHp9w7Dgn+CI1sSp4e31eFsOUxNZCyn2aJUlONm6dwK/rxuP997YjNFOW7Ory6Nfz8wZVwp7+lp+A70sLTREMaIkiE1xlMyAIF6q2ROzsjuJXOSUZTjZreMJ89/1AqJDQCeYBNtvfCUSvK9LKgq5cmNB9AZn4TiSqumXTjUvfMRK5V7q46nurxz7erfFp3OmeOKONjQxq0fm4bPnTiE6HE52JM/m9HNW8BeB0vE5qd/wfpXk4cThwJGlAypMaJkACKNB2nSHIpH9Cx8JyLU5Vdbb44MTPUCT7iZdmfPPSWAq2aNZd+JNt7a0wiX3QZH34f1v+ve8Yi16fWAZzKl+Z1nMxX43PzlH89j039fwSfmdK/k0I2J5+EiTMvOtxJ2qT+4i+nv3MJZr36B1rbk4jXYSUuURGSRiGwTkRoRuSXOda+IPGpfXyUilTHXvmO3bxORK1LZFJFJto0dtk3PKYwxQkRWiMhWEXlfRD7Us6/HYDBEcTQf4rAWU5zXiwPwyq3qBRwZmD05vnBLr0XpypljGJnn4Vev74YzroIJ58FLt0HXSueH3qNF8igsn4iInHRJRPC60kuyqJh5MWEVDm98KWGffe93rjltX/tquh9lUJJSlETECfwc+AgwDbhGRKZ16XYDUKeqVcA9wJ32vdOAZcB0YBFwn4g4U9i8E7hHVauBOtt2j8ew7/kJ8JyqTgVmAWbrdC8wKeEGAHfrYQ5TTHFu99NcU1E6bjL1mkeodn0/zKw7OZEWQu7eiZLP7eTacybw0tbD1Bxtho/9GAIt8OQ/nRw12PMm6/R0qkaln3Yej5lV49lIFd7dLyfs03648+DBxl1rT2m8bCcdT2k+UKOqO1U1ACwHlnTpswR42H69ArhUrF8dlgDLVdWvqruAGtteXJv2PZfYNrBtXt2bMUSkELgQeAhAVQOqWp/e12IwGLqS036EBmcJ7jQrhMdSPaqQdZFqQnsGIMssEqZQGwl4exZmjOXz51WS53Fx53PboOx0uPR7sP1ZePNnVoe6PXBsO38PTmXKKYqSz+2kZsT5jGt9H208GL9T3S6aNIc6inAe3nRK42U76fztGgfEntlba7fF7aOqIaABKElyb6L2EqDettF1rJ6OMRk4CvxaRNaLyIMiEnd7t4jcJCJrRGTN0aMDW2J/MGCqhBusY72P0eorT903DlNG5bMmMgVffQ20nujjyXWh9bhVeSI3cdZbKkryvXzlotN4ccth3vrgOJz7jzDtanjxe/D2L+C1H6IIT4fP5fTRpyZKAL4ZiwE4su7JuNdzmvZwyDmGQ7nVlLZsP+Xxspl0REnitHV9SiXq01ftvRnDBcwBfqGqZwEtQLf1MABVfUBV56nqvLKynqW7DguMKBlajuEiRDB3VK9uryzN412Zar3p6SF6PSTYaJ1PpLm9E9AoX1wwiQkjc7nl8Y20hiLw8fuh+nJ47hZY91veH/8ZDlB6yp4SwPz557NfS2h576m410f491OfU0F78elMjOyjsbX9lMfMVtIRpVog9oSqCuBAoj4i4gKKgBNJ7k3UfgwYYdvoOlZvxqhV1eixkCuwRMpgMPQUuwxOpGhCr253Ox0Ex8ylTXJgR5rHQvSS1hPWI8NR2DsBjZLjcfK/S2ey53grtz/9Prhz4JrlcO0K+MzveSD3JsoLvCdl3vWWUUU5rM9dwLjjb0F740nXNBxkVPgw/oKJeMediU+C7N7+3imPma2kI0rvANV2VpwHK6lgZZc+K4Hr7ddLgZfVivmsBJbZmXOTgGpgdSKb9j2v2DawbT7RmzFU9RCwT0ROt++5FDAHl/QCE74z+I/WAOAomdxrGzMnlvO38Ex023PQj0V+/Sf2A+As6tl+qnicM7mEL184mT+u2suj7+wFhwOqL0OnXsnbu+qZP2lkt8y73hI84xN4CXBk9Z9Oaq87uAu3hJGRkymrOguAEzs39MmY2UhKUbLXb74KPI+VvfaYqm4WkdtE5Cq720NAiYjUAP+MHSZT1c3AY1hi8Bxws6qGE9m0bX0b+GfbVoltu8dj2Pd8DfiDiGwEZgPf7+kXZDAYoOWQJUo55b0XpTkTi3kuNBdpPgS1/VfBIHy0hqA68ZRU9om9f73idC6oLuU//28TL71vpYVv2FfPocZ2Lqzuu3D/eR9exG4dRevak6uTH99nVSL3jaqitHImYYTQoaGb7NC9nnocVPUZ4Jkubd+Led0OfCrBvXcAd6Rj027fiZWd17W9N2NsAObFu8fQA4ynNOwJHPmAgzqSsaU9qxAey7zKYr4VmUfAmYdnza9hwrl9OMMYjtewT8sYNSL9YyuS4XI6uPeaOXzuV6v48u/WcvPFVbz1wXHyvS4WnTm6T8YAK4T3f0UL+VjDI0QaDuAoGgtAy0ErsWHk+KmIJ5fDrnHk1W/rs3GzDVPRwZASxZynNNyR+l3s1XIqS3teuidKeYGP08aN4iXPJbD5L9B8pA9n2ImvoYadOobRhb4+s1mU6+YPXzqHy6aN4icv7WDNnhP811XTKfT1YiNxEgrmX4uTCHtf/mVHW+RYDW3qYUyF5aU2FFQzxr+LUHho/rs0omQwGJKjSkHTLvbKOMryvadk6uKp5fxvw0VoJAR/u6uPJhhDWz3FLTt5j2qKcvpYMHxufnHdXP72rxfxxi2XsHRuGiWEesiCc87lLWZStOm3HfX2vI272e8Yg89jBbYiZdOYwGH2HBqa21eMKBlSY8J3w5umQ+SGGziaV33Ki/qXTxvFzsgYtld8Etb+Gg708YL97tesP3LP7LMEhK5MLMljTFFOv9j2uZ0cOv0fKA4f48g7fwagtLWGo76JHX0KJ87CIUrt9oGpjjHQGFEyGAzJOWzlILWPPOOUTc0YV8T0sYXc2nQ15JXDn284tc20J3bCDybAYTuxdsMfqZci6koH7+6PBR/5LHu1nPCrdxGs38+oyBFaSmd3XB9VbX22lr0bMzXFfsWIkiE1xlMa1gQPWHtivONm9Im9ZWeP5+1D8P6Cu6F+L/z+E92LnabL6l9CewP8ahE89++w7Rn+ELmCyvLEp71mO+Uj8nhj4j8ypr2Glgc/BkDOlIs6rntKJ9OOB+exoVnK04iSISVGk4Y3zbvXUKulVE3omzWUpXPHU5rv5baNxfDp31nHWfy/C2DjY+nvXwq2wZv3WqIE4HTD2z+nfcpV/MS/mNPKeleMNVu4/FNf4RXmMaL5A9boGZw1/8LOiw4nR3yTKG4euEMTB5K0UsINBsMwRRXvgdW8Hqli9pjCPjGZ43Fy88Wn8d9PbuGZD83hoze+BH/5/+DxG60D9c74GFScbR2wlzvS+q0o5IemA3BiF+xbBTUvQdsJqL4CPvEAePIh2MIbu9oJbFzD1D6oR5dJSgp8lHzxUe554UnmnnMRud6Tkzbaiqcy+cArHGtqp7Sg77IMswEjSoY0GJqpp4Y0aNhHrv8I7zmv5Mrivlvc/9y5E3l83X7+8/82cdbXz2fMTX+DrU/C+j/A2t/AqvsT35w/CqoWwpzPQeUFEE1ocBaxZs9B3E5h1vjBG76LMnNCKTO/9IW411wVZ1F68AlW17xP6VlnDfDM+hcjSoaEiDlHybDXKh3ZVDa3T7PZXE4Hd396Fh+/702+8Ot3WH7TuYyYtgSmLYFwEI5uhcYDVhKEiBWeKxgLxROhYEynEHXh79uPMrNiRNIjyIcC5dM+DO9A3dbXwIiSYdhhFpWGLYHtf6VV8zoyvvqS6lEF3HftHL708Bo+8Ys3eeBz86gqz7cEaPSZ1k8PqDnSxOYDjfzn4lPPEsx2CibOopUcHLWrUnceZJhEB0NC1D4VxBRkHaZEwrDjBV6NzGL+5FM7BiIRF04p43c3zOdES4CP/vQ17n5xO/WtgV7Z+tnLNfjcDpbM7nrc2xDE4eRgwQzGN28kEBpa4XUjSoaEmPDdMGf/Wjz+E/xN53DWhOJ+G+acySW88M0LuWzaKH760g7O+8HLfPPRDTy36RB1LakFSlX5zRu7eGLDAW68YDJlBadWdWKwEBl/DlPYx+Zd+1J3HkSY8J0hNcZTGp68t4IALlomXEKOp3/XaMoLfPz8s3P42iWN/Pr13Ty3+RB/WW8dQVFZkktVeQETRuYydoSPfK+LHI+T9mCYgw3tvLL1CO/WNnDp1HK+fml1v84zmxg142IcW+6l9t1XOav6+tQ3DBKMKBkMhu6E/ITffZQXwvM4/8zTBmzYqaMLuXPpTG7/+AzW7qlj3d463t1Xz57jrbz5wTFaA+Fu90wbU8jtV8/gmvkTcDr6p7RQNlJYvYAAbpy7/0bnUXODHyNKhoQ47PCdmjDe8GPrUzj99fwp/GF+MO3UTnDtDW6ng3Mnl3Du5JKONlWlsT1EayBEiz+Mz+1gZJ6HXM8wfYy5c6gtmMVpje/QHgwPmYxDs6ZkSIItRiZ8N7xQRV//MftkDMGJH+634qM9RUQoynEzpiiHqvJ8Kopzh68gRTntYk6XvazbvDXTM+kzjCgZEmISHYYpNX9FDm3kp4GPsezcSZmejSEJ4+Z8BIBDG57L8Ez6DiNKhoQ4OkRpaKWcGpIQDsLz/8Fh52he813CFdMHPnRnSB9vxVk0OQrJrX0901PpM4woGZJgPKVhx6r74dg2/qPtWm68eCpe19BYpxiyOBwcKzuX2cH17D3WkunZ9AlGlAwJEXstySwpDRMOvYe+dBur3PPZmPshrj1nQqZnZEiDvOmLGC11bFw3NLyltERJRBaJyDYRqRGRW+Jc94rIo/b1VSJSGXPtO3b7NhG5IpVNEZlk29hh2/T0dgz7mlNE1ovIU+l/LQaAjuRao0pDn7Y6+NMXaHUW8pWmL3LrVTOGTDbXUKfsrMUABLc8k+GZ9A0pRUlEnMDPgY8A04BrRGRal243AHWqWgXcA9xp3zsNWAZMBxYB99kikczmncA9qloN1Nm2ezxGzNz+CRiap2H1MybRYZgQbIdHPkukbg83tt7MWVOr+OiZozM9K0OaSMFoanPPYFLdG/hD3fdxDTbS8ZTmAzWqulNVA8ByYEmXPkuAh+3XK4BLxSopvARYrqp+Vd0F1Nj24tq077nEtoFt8+pejoGIVACLgQfT+zoMsYhJCR/6BFpg+Wdh75v8l/Nr7Mmfzf9+alafVgQ39D+ByQuZSQ3r3x/8B/+lI0rjgNjiSrV2W9w+qhoCGoCSJPcmai8B6m0bXcfq6RgAPwb+jRTpYyJyk4isEZE1R48eTdZ1WGE8pSFOy3H47dXozlf4Uc7XWOGfzy+um8PIPE+mZ2boIWPPvhqHKAfXPp3pqZwy6YhSvF+Zuj6tEvXpq/YejyEiVwJHVHVtnOsnd1Z9QFXnqeq8srKyVN2HDVFRUjUp4UOO2rXw/y5ED77LbTn/xgNNC3jw+nnMrBj8h+MNR3zj51DvKKao9uVMT+WUSUeUaoHxMe8rgAOJ+oiICygCTiS5N1H7MWCEbaPrWD0dYwFwlYjsxgoPXiIiv0/j8xpsTABnCBIJwxs/hV9dgT8C/yD/w4rWOfz682dz3mmlmZ6dobc4HBwZ/WHmBdex50h9pmdzSqQjSu8A1XZWnAcrqWBllz4r6awIuBR4Wa1DeFYCy+zMuUlANbA6kU37nldsG9g2n+jNGKr6HVWtUNVK2/7Lqnpdmt+LARO+G3Ic3gIPLoQXv8sHxQv40Ilb2e2p5vGvnMd5VUaQBjsjZl1JobTy/uq/Znoqp0TKwlGqGhKRrwLPA07gV6q6WURuA9ao6krgIeB3IlKD5b0ss+/dLCKPAVuAEHCzqoYB4tm0h/w2sFxEbgfW27bpzRiGU6MzfGfEaVDTegL+/r+w+peEPAXcU/Btfr5/Jktmj+N/rp5Boc+d6Rka+oDy2YsIPusisv15On+vH3yIeeCczLx583TNmjWZnkZWELq1GJdEWPehe5lzxecyPR1DTwkFYM1D8OoPUH8ja4sX8+WDi5G8Mr575RnD44TWYcbOH11KuPEQ4/9z44DvMxORtao671TtDPMSu4ZkRD0lMbXvBhchP6z/Pbz+Y2jYy56i+XyjdSnvHqzg+vMq+eZlU4x3NETR6supXvd93tq4gQ/NnZvp6fQKI0qGhEQTHYw3PUgItsG631pi1HSAQwUz+G/+nWcPT2fxmWO5a2E11aMKMj1LQz9Scc7HYd33ObH+KTCiZBhqOMSI0aCg+Sis/TWs/iW0HGFP3kxuC3+Bl45O49Kpo3j68ilMH1uU6VkaBgDvqCkcdo1lxIG/o6qDchO0ESVDaoynlJ0c3gxv34du/BMS9rMpdz53BG5kTXA6S2ZX8Oz5kzhjTGGmZ2kYYOrGXshZe/7CrkPHmTxm8GVVGlEyGAYTkTDseBHevg92/Y2Qw8ezzov5ceulHJNKrvvwBH7yoUrKC32ZnqkhQ5TMWkzu3uVsW/0Ck5d8NtPT6TFGlAwpMY5SFtB4ANb9Dl33W6SxlnpXGQ+Gr+H37Rcxcfx4vnz5BK6cNcYcD26g7MyFBJ50W7+8YETJMCQx2XcZIRKGmpdg7W/Q7c8hGmadaza/DCzlbZ3P4jkT+MM5E8x6keFkPLnUFs6humEVLf4Qed7B9ZgfXLM1GIYDjQdh/e/Rtb9BGmtpcIzgj6HFPBK6mBGlU7jm8gn8aNbYQfewMQwcUr2QqrV38NrGjVxw9pxMT6dHmL/VhtSY+F3/E2yDrU+j7z4CH7yMaIRVnMlvA59gY94Crjy7kgfnjGOKSek2pMG4s6+CtXdwYsPTYETJYDCkhSrsWwUb/khk019wBBo5IqU8FryKp+Rips2YzTVzxvGz00pxOgZfaq8hc3hGnc4x12iKDw6+1HAjSoaUqCnM2rfU74V3HyW0/g+46nfRjpenw2fz5/CHkcrz+fjcCfx5xmjyTXjO0FtEqBv7YebseYIdB08wZWxJpmeUNuZvvSE1RpNOnbZ62PoUwXWP4N73OgBrImewIvxl9oy6lIWzqvjfWWMZNyInwxM1DBVKZn+U/L2Psm31C0y5+ppMTydtjCgZ0sCoUq8ItsH25whueAzHBy/ijAQ5oOX8ObSUdcWXM/+sOfzjzDFMLsvP9EwNQ5CR0xcSWOmyU8ONKBmGEibRIX3CQdj5KsENjyLbnsYVaqVOR/BU+FJW5V1C1VkX8rHZ4/jmqIJBFec3DEK8+RwonE1VwxraAmFyPANbNby3GFEyGE6VSAT2voV/w2Ow+f/wButp1TyeDc/ndd9FjJm1kCtnj+cLFUVGiAwDik66kDPevZu3t+3g3DOnZno6aWFEyRAXjUQ6j0M3nlJ3VGH/Wlo3rEA3PU5e+2Ei6uGvkTm85r2IwjMXcfnMiXx6YjEOkzlnyBBjZl8B797N4XdfBCNKhsFMOBzq+MthJMkmEoHa1bSs/zP6/kry2w/hUievRWbyZu5nyTnzYyycdRp3Go/IkCX4JsyjVXLx1b4OfC3T00kLI0qGuISCgY6/HDKcZSkShr1v0bTuzzi2riQvcAyXuvh7ZBbr8q4hf+aVXDx7Cv85xqwRGbIQp4sDxfOYenwdje3BQXG4oxElQ1yCwQDROtPD7pC/cIjwrtc5seYxcj94lrzgCdzq5tXIbDYWfp6iWVdy6exqLis3WXOG7Md12kVMPPF3XntvIxecnf0H/xlRGuJs/OFigpMXMvcT3+zRfaGAv59mlKUE22nd/jLH3nmckn0vkBduIE+9vBKZzfaRl1I486NcNvs0FpXkZnqmBkOPGDvnCnjnNuo2vQhDRZREZBHwE8AJPKiqP+hy3Qv8FpgLHAc+o6q77WvfAW4AwsDXVfX5ZDZFZBKwHBgJrAM+p6qBno4hIuPt/qOxylw/oKo/6ekXNNiZ2fw6bHwdeihK/rbmzjdD1VNqOszRdU/Q8t7TjD7+FrnqZ6T6+JvM5eC4Kxh11mLOnzaBxbmeTM/UYOg1ntHTqXeMpODAG5meSlqkFCURcQI/By4DaoF3RGSlqm6J6XYDUKeqVSKyDLgT+IyITAOWAdOBscBfRWSKfU8im3cC96jqchG537b9i16MEQL+RVXXiUgBsFZEXuwy7yFN7LpQT/G3NqfuNNhQJbh/IwdXP46z5nnGtb5PGRDQEGAq5AAAEPZJREFUEl7wXkrbpMuZfPYiLp80CpfTkenZGgx9gwiHS89hxuE3OdbUTmlBdh8Amc4zaz5Qo6o7AURkObAEiH24LwH+y369ArhXrFXfJcByVfUDu0SkxrZHPJsi8j5wCZ0nUz1s2/1FT8dQ1beAgwCq2mTbHtdl3kOaluZGenvSTrA9VpQGsacUbKfh/Zc4tvYJRu5/meLQUSpU2Kin8Vbx5/FM+yhnzTufq0ryMj1Tg6Hf8E25mLIjz/K3DW/z4QsuyvR0kpKOKI0D9sW8rwXOSdRHVUMi0gCU2O1vd7l3nP06ns0SoF5VQ3H692YMAESkEjgLWBXvA4rITcBNABMmTIjXZVDS3lTXa1EKtLd0vB5siQ7adIgD7zyBf/PTjD3+NkX4camXdxyzODbuBkrnfIyzZ0xltil4ahgmjJuzCF7/N5refwmGgCjFy3Pt+pRK1CdRe7zYSLL+vRnDukkkH/gz8A1VbYzTF1V9AHgAYN68eYPrCZyE9ta4HzctQjGilPWo0l77LvvffhzPzhcY3/Y+44D9WsLLOQsJnHYF1fMXceH4crOR1TAscY2cyBHXGEYcWZ3pqaQkHVGqBcbHvK8ADiToUysiLqAIOJHi3njtx4ARIuKyvaXY/j0eQ0TcWIL0B1V9PI3POqTwtzb1+t5wW6ygZaFOB9s5tulFjq9bSemBVygJH2WSCu9xGutKvohvhhWW+2ihqbptMAAcLz2b6Qdf5kRzOyPzs3ddKR1RegeotrPi9mMlFXy2S5+VwPXAW8BS4GVVVRFZCfxRRO7GSkKoBlZjeTfdbNr3vGLbWG7bfKI3Y9jrTQ8B76vq3T39YoYCwZaG3t/bfKwPZ9I3hBsPsW/VXwhueYbxdasoxU+Oelnrmk3dxBsZNfcqzpo2hVmuwVF40mAYSLynXUDxoZW8sXE1C867MNPTSUhKUbLXb74KPI+Vvv0rVd0sIrcBa1R1JdbD/3d2ksEJLJHB7vcYVnJBCLhZVcMA8WzaQ34bWC4itwPrbdv0dAwROR/4HPCeiGywbfy7qj7Tu69q8BFoszyldnXT09+LIjGiJJlaU1Klac969q96nJxdLzKxfSuVwAEt4W95lxGpXsSUcz7CBWNKTDUFgyEF42YvhDegceurMJhFCcB+kD/Tpe17Ma/bgU8luPcO4I50bNrtO+nM0Itt79EYqvo68debhg3hdkuUAtJzUZK2Ex2vBzLRQYNtHNzwAnUbnmTUwVcpjRxligqbpIpnyr5E7plXcta8BVxh9g4ZDD3CWzqJY45SCg7FzffKGkz60RAm3FoPQICeP8C9LftpUR950t7X0+qGv/4Ae976C5GtzzKxYTVj8VOkXjZ4zuKdypsYO38JM6ZMYaZJUjAYeo8IR0bO4/Sjb9H8/7d370FR3ucCx78Py0XjDeIV8QKIclERAZGAuWo10bZJk5zGTHKSOc1p0pO2czo9c1LTZJqTzEknpmea00vajG06rU1bY9OkGtvGoLmZxLvIAgEUiaKCIsrFoNx//eP9KQRBF1zYhX0+Mzu8+9vfe3t099n3fX/7vI0tjPTTOnialIYw81kVAOel96Vxws8fpSI4iplth/D6QAdjOFO2l4qdbzDiyBZimoqZhXNa7uNRSyH+NhKuW072uAjvrlepABccu4jx1W+xu2AfC9K7/rLHP2hSGsKCzjlJqbdVvpsbzzG5rQL36Bvg7CGvbItpOc+RvZup3/8mkVXvM779FOFG+CQojrcnfZ0x875Ecmo2i/W3Q0r1m6iUJbDrSWqK3gNNSmqgDW9wRtO7Lv4W2TNl+98nQVppnXIdFL3T59p3jfXVlH38Ou1FfyOmbgfRNNJgwsgPS8Ud/Q2iMu5g9owZzNFBCkoNiBGRCdRIOMOPb/f1pvRIk9JQZQyRTWUAuGjr1az1e9bTaEKYOHcxFD3bq4EOtccPcuSj9Qwv20zs+XySpJ0qE8GeMUuQhBUkZa0gM7yvdSaUUldFhMrwVOLO7KeppZWwEP9LAf63RcorjhTtYTrOQIcQPD9SqiovYe6pTeSNuYXY8LEezVNfVU7Ze68wsnQDcc3FhAOlMo1tE+9nVMrtzEm/kRtD/fOiqlIBJ3oRkTXv4C4uJHnuPF9vzSU0KQ1BjQ311L35fRpNCIWjspn5mWelRRrqTnN27b2MQJh617PIhRH1pv2Svu2tLZR8sB52v0z8uX2kiOGAxPLutG8RmflV4hOTidPTckr5nah5iyH3GU4VbAVNSqo/nTtbQ8GGF5hR+luSqePjhMcJqj1M8Nkrn777dP/7hG58hGltVeQt+iXp02dSd/okAKa144Z/bS3NuDf9nMnuF0k01VQylm2Tv8bErPuIn5PKLE1ESvm1MdOSqZdRhB7zz+tKmpSGgPozVRT99XkSyv9IBg24w1KpvHkVWZnL2LHm2wRf5vTd0eLdnHnzKeY1fMQpIiha+gfSs28DIHSYM5T8QlI6UrgDXv8689vKKXQl8mnaU6Qsvocbw8L6fyeVUt4RFMTRUSlE1+XR1m5w+dnv/zQpDWKnK8sp3fgccyr+wkJpJHd4FsOXPEZy2s0dnYJCCJU2THs7EtRRnL2yrJDKDT8gpXYr4Qzjw2kPM/fOx0mOuPZin9AwW8y0tYmiD/9KTM6/c1ZGsnPhz1iw9H6C9EZ4Sg1K7dOymFqwjaKDJSTGJ/h6cz5Hk9IgdOygm4q3/o+U6r+RThu5o28mYtkq5s+59HcHxuUMMGhtbSEkNIzTlUf49LUnSanexBhC2B55Pwl3Pcmi8ZMumdcVHEyLcRFaW8rULWs57opi1MObWDhp6iV9lVKDx+SUL0DBak7k5WhSUn3T1NhAwZZXGO5+haRmNxNMMLljVxC14nukz5jd43zickoMtTQ3kvvaauaW/Ixk2tg1/k7i7nyK7MmXv6lhMyGk1W+lwQwj7F9fZYImJKUGvbGxadQzkpDybcC3fb05n6NJyY81NTZQ8tGbNOVvYFbt+6TRwHGZyPbobzJz2TdYGOnBXXJdzj9x0a8eIqMuh9xrMhl714/Jius5kXV2ofZdXuTdZMX41zcqpVQfBQVRPjqVmPp9fnddSZOSn2mor6HkwzcwRRtJqN9OsjRSb66hZEw2YQseYE7WCqJcnt8vSOzpu7S6HLZPvJfMR37xuWtLnopd8d1ez6OU8l/t0dcT5f6AkpIC4hPn+npzLtKk5AfqTp/gwAd/JvTAJhLO7SVVWjjNGArGLmVY8h0kXLecBWF9u4OqBDs3rWgyIaQ/9JNeJ6S8G9bQVH+KjKkz+rR+pZR/ikxZBu5nOel+W5OSguqKwxza9iojy/5OfKObBdLOCcaRO/FORs//CvELvsDC4Kv/55kw5yaa8n+Ie/4zLAjt/dDtebfcc9XboJTyP+Njkjkj4YSWfwT8l6835yJNSgPoRPkBDm9bR/jhfzCruYhxYigPimJ31AOMy7ibuORsJvXh1NrlRCek0vJkBQtC9KZ4SqlORCgfnUZs3T7a2tpx+clPPDQp9bPzDWcp3LKWawr/RFJzPpOAQ0Ex7Ix+hMmZX2V6YhoeDFe4KiGakJRS3TDR1zMhbysHinOZNTvN15sDaFLqN7WnKine8DxJx9aRzjmOSSTbox9l6qL7mBE3B71Co5TytajUWyHvf6ja/5YmpaGqrbWVPa/9iNlFPyFTzrNvxPWEZT9KUuatTPHyqTmllLoaE6YncjRoCqOPbAae8PXmAODRp6SI3CoiJSJSKiKrunk9TEReta/vFJHoTq89bttLRGTZlZYpIjF2GQftMkO9vY7+Unf6JJ/8aAkLi5+jbHgSh+/ZSup/b2J21vI+DcNWSqn+Vhm1lKSmfE6dPO7rTQE8SEoi4gJeBG4DkoB7RSSpS7eHgBpjTBzwArDazpsErARmA7cCvxAR1xWWuRp4wRgzE6ixy/b2OryutvoENS8uJr4xn11zn2buY1uITkzvr9UppZRXjMv4F4KlnYPvrPX1pgCeHSllAKXGmDJjTDOwDri9S5/bgd/Z6deAxSIitn2dMabJGPMpUGqX1+0y7Ty32GVgl3mHN9fhWVh6p7mpkaqXvkRk2wkOLP0tGXd9R4+MlFKDQszsTIpDkkgseZETR0t9vTkeXVOKAo52en4M6Fr582IfY0yriNQBY237ji7zRtnp7pY5Fqg1xrR2099b67iEiDwMPGyffiYip4Hq7vpe0TNf7tNsfmocfY3D0KOx6KCxcAy9ODw9s69zjgOme2MTPElK3RVFMh726am9u8OIy/X35joubTRmDbDmwnMR2WOMCfhzbxqHDhqLDhoLh8ahg41FtDeW5ck5pmNA59LQU4CKnvqISDAwBjhzmXl7aq8Gwu0yuq7LW+tQSinlpzxJSruBmXZUXCjOoIKNXfpsBB6003cD7xhjjG1faUfOxQAzgV09LdPO865dBnaZG7y5Ds/CopRSyheuePrOXr/5FrAZcAG/McYUisgzwB5jzEbgZeD3IlKKc/Sy0s5bKCLrgU+AVuCbxpg2gO6WaVf5PWCdiPwvkGuXjZfXcSVrrtwlIGgcOmgsOmgsHBqHDl6LhTgHG0oppZTv6bhlpZRSfkOTklJKKb+hSamTgS5L5Asi8hsRqRKRgk5t14pIji3tlCMiEbZdROSnNh5uEUntNM+Dtv9BEXmwu3X5MxGZKiLvikiRiBSKyH/a9kCMxTAR2SUieTYWT9t2r5X8GkxsRZhcEdlknwdqHA6LSL6I7BeRPbat/98fxhh9ONfVXMAhIBYIBfKAJF9vVz/s5w1AKlDQqe15YJWdXgWsttPLgX/g/BYsE9hp268FyuzfCDsd4et962UcIoFUOz0KOIBTjioQYyHASDsdAuy0+7geWGnbXwL+w04/Crxkp1cCr9rpJPu+CQNi7PvJ5ev960M8vgv8EdhknwdqHA4D47q09fv7Q4+UOgxYWSJfMsZ8gDN6sbPOJZy6lnZaaxw7cH5DFgksA3KMMWeMMTVADk7dwUHDGFNpjNlnp88CRTiVQAIxFsYY85l9GmIfBu+V/Bo0RGQKsAL4tX3uzdJnQ0G/vz80KXXorpxSVA99h5qJxphKcD6sgQm2vaeYDKlY2dMu83GOEAIyFvaU1X6gCueD4xAelvwCOpf8Guyx+H/gMaDdPve49BlDKw7gfDF5W0T2ilOKDQbg/aH3U+rgSTmlQNPb0k6DjoiMBP4CfMcYU+980e2+azdtQyYWxvltX4qIhANvAInddbN/h2QsROSLQJUxZq+I3HShuZuuQzoOnWQbYypEZAKQIyLFl+nrtVjokVKHQC5LdNIeamP/Vtn2IV3CSURCcBLSH4wxr9vmgIzFBcaYWuA9nOsC3ir5NVhkA18WkcM4p+9vwTlyCrQ4AGCMqbB/q3C+qGQwAO8PTUodArksUecSTl1LOz1gR9ZkAnX2kH0zsFREIuzom6W2bdCw5/5fBoqMMT/u9FIgxmK8PUJCRIYDS3CusXmr5NegYIx53BgzxTiFRVfi7Nd9BFgcAERkhIiMujCN8/+6gIF4f/h6hIc/PXBGkBzAOZ/+hK+3p5/28U9AJdCC8y3mIZzz4FuBg/bvtbav4Nwo8RCQD6R3Ws7XcC7glgL/5uv96kMcFuGcRnAD++1jeYDGIhmnpJfbfvD8wLbH4nyYlgJ/BsJs+zD7vNS+HttpWU/YGJUAt/l6364iJjfRMfou4OJg9znPPgovfB4OxPtDywwppZTyG3r6TimllN/QpKSUUspvaFJSSinlNzQpKaWU8hualJRSSvkNTUpKKaX8hiYlpZRSfuOf0QlWeWtOd2MAAAAASUVORK5CYII=\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 (Ctt = 0.0)')\n",
"plt.plot(test_q, calcs_test1, label = 'pdf (Ctt = 0.5)')\n",
"# plt.plot(test_q, calcs_test2, label = 'pdf (Ctt = 4.9)')\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",
"constraint7 = tf.cond(tf.greater_equal(Ctt*Ctt, 0.25), lambda: 100000., lambda: 0.)\n",
"\n",
"constraint7dtype = tf.float64\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]#, constraint7]"
]
},
{
"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": [],
"source": [
"# # zfit.run.numeric_checks = False \n",
"\n",
"# fitting_range = 'cut'\n",
"# total_BR = 1.7e-10 + 4.9e-10 + 2.5e-9 + 6.02e-5 + 4.97e-6 + 1.38e-9 + 4.2e-10 + 2.6e-9 + 6.1e-10 + 4.37e-7\n",
"# cut_BR = 1.0 - (6.02e-5 + 4.97e-6)/total_BR\n",
"\n",
"# Ctt_list = []\n",
"# Ctt_error_list = []\n",
"\n",
"# nr_of_toys = 1\n",
"# if fitting_range == 'cut':\n",
"# nevents = int(pdg[\"number_of_decays\"]*cut_BR)\n",
"# else:\n",
"# nevents = int(pdg[\"number_of_decays\"])\n",
"# # nevents = pdg[\"number_of_decays\"]\n",
"# event_stack = 1000000\n",
"# # nevents *= 41\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",
"# if fitting_range == 'cut':\n",
" \n",
"# sampler.resample(n=nevents)\n",
"# s = sampler.unstack_x()\n",
"# sam = zfit.run(s)\n",
"# calls = 0\n",
"# c = 1\n",
" \n",
"# else: \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",
"# if fitting_range == 'cut':\n",
" \n",
"# total_samp = sam\n",
" \n",
"# else:\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 - 60.)))\n",
" \n",
"# tot_sam_1 = total_samp[_1]\n",
" \n",
"# _2 = np.where((total_samp >= (jpsi_mass + 70.)) & (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",
"# plt.clf()\n",
"# plt.hist(tot_sam, bins = int((x_max-x_min)/7.), label = 'toy data')\n",
"# plt.savefig(plotdirName + '/toy_histo_cut_region{}.png'.format(toy))\n",
"\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.axvline(x=jpsi_mass-60.,color='red', linewidth=0.7, linestyle = 'dotted')\n",
"# plt.axvline(x=jpsi_mass+70.,color='red', linewidth=0.7, linestyle = 'dotted')\n",
"# plt.axvline(x=psi2s_mass-50.,color='red', linewidth=0.7, linestyle = 'dotted')\n",
"# plt.axvline(x=psi2s_mass+50.,color='red', linewidth=0.7, linestyle = 'dotted')\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)/(toy+1))*((nr_of_toys-toy-1)))))\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [],
"source": [
"# with open(\"data/results/Ctt_list.pkl\", \"wb\") as f:\n",
"# pkl.dump(Ctt_list, f, pkl.HIGHEST_PROTOCOL)\n",
"# with open(\"data/results/Ctt_error_list.pkl\", \"wb\") as f:\n",
"# pkl.dump(Ctt_error_list, f, pkl.HIGHEST_PROTOCOL)"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [],
"source": [
"# print('{0}/{1} fits converged'.format(len(Ctt_list), nr_of_toys))\n",
"# print('Mean Ctt value = {}'.format(np.mean(Ctt_list)))\n",
"# print('Mean Ctt error = {}'.format(np.mean(Ctt_error_list)))\n",
"# print('95 Sensitivy = {}'.format(((2*np.mean(Ctt_error_list))**2)*4.2/1000))"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [],
"source": [
"# plt.hist(tot_sam, bins = int((x_max-x_min)/7.))\n",
"\n",
"# plt.show()\n",
"\n",
"# # _ = np.where((total_samp >= x_min) & (total_samp <= (jpsi_mass - 50.)))\n",
"\n",
"# tot_sam.shape"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [],
"source": [
"# Ctt.floating = False"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {},
"outputs": [],
"source": [
"# zfit.run(nll.value())"
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {},
"outputs": [],
"source": [
"# result.fmin"
]
},
{
"cell_type": "code",
"execution_count": 46,
"metadata": {},
"outputs": [],
"source": [
"# BR_steps = np.linspace(0.0, 1e-3, 11)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# CLS Code"
]
},
{
"cell_type": "code",
"execution_count": 47,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[0.15430335 0.53452248 0.74001287 0.89973541 1.03509834 1.15470054\n",
" 1.26302735 1.36277029 1.45569489 1.5430335 ]\n",
"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",
"Step: 1/10\n",
"Current Ctt: 0.1543033499620919\n",
"Ctt floating: True\n",
"Toy 0/1 - Fit 0/5\n"
]
},
{
"ename": "RuntimeError",
"evalue": "exception was raised in user function\nUser function arguments:\n bplus_0 = -0.627540\n p4160_s = +3.373292\n DDstar_s = +0.299698\n p3770_s = +3.780920\n jpsi_p = -5.063578\n p3770_p = +2.196764\n bplus_1 = +1.062520\n DDstar_p = +4.628551\n omega_s = +5.708895\n phi_s = +16.739100\n Ctt = +3.195816\n p4415_s = +0.911410\n psi2s_p = +2.482877\n Dbar_p = +0.982372\n p4040_p = -2.632544\n p4160_p = +2.741969\n p4415_p = -2.983702\n phi_p = -4.941643\n rho_p = +1.587532\n rho_s = +0.149543\n bplus_2 = -0.691416\n Dbar_s = +0.150965\n omega_p = -0.692164\n p4040_s = +1.222875\nOriginal python exception in user function:\nKeyboardInterrupt: \n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\minimizers\\minimizer_minuit.py\", line 101, in func\n loss_evaluated = self.sess.run(loss_val)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\client\\session.py\", line 929, in run\n run_metadata_ptr)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\client\\session.py\", line 1152, in _run\n feed_dict_tensor, options, run_metadata)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\client\\session.py\", line 1328, in _do_run\n run_metadata)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\client\\session.py\", line 1334, in _do_call\n return fn(*args)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\client\\session.py\", line 1319, in _run_fn\n options, feed_dict, fetch_list, target_list, run_metadata)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\client\\session.py\", line 1407, in _call_tf_sessionrun\n run_metadata)\n",
"output_type": "error",
"traceback": [
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[1;31mRuntimeError\u001b[0m Traceback (most recent call last)",
"\u001b[1;32m<ipython-input-47-d417fd816fd8>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m 101\u001b[0m \u001b[0mminimizer\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mzfit\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mminimize\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mMinuitMinimizer\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mverbosity\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m5\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 102\u001b[0m \u001b[1;31m# minimizer._use_tfgrad = False\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 103\u001b[1;33m \u001b[0mresult\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mminimizer\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mminimize\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnll\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 104\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 105\u001b[0m \u001b[1;31m# print(\"Function minimum:\", result.fmin)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32m~\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\minimizers\\baseminimizer.py\u001b[0m in \u001b[0;36mminimize\u001b[1;34m(self, loss, params)\u001b[0m\n\u001b[0;32m 205\u001b[0m \u001b[0mtuple\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mstack\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0menter_context\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mparam\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mset_sess\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msess\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mparam\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mparams\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 206\u001b[0m \u001b[1;32mtry\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 207\u001b[1;33m \u001b[1;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_hook_minimize\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mloss\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mloss\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mparams\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mparams\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 208\u001b[0m \u001b[1;32mexcept\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mFailMinimizeNaN\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mRuntimeError\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mas\u001b[0m \u001b[0merror\u001b[0m\u001b[1;33m:\u001b[0m \u001b[1;31m# iminuit raises RuntimeError if user raises Error\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 209\u001b[0m \u001b[0mfail_result\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mstrategy\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfit_result\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32m~\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\minimizers\\baseminimizer.py\u001b[0m in \u001b[0;36m_hook_minimize\u001b[1;34m(self, loss, params)\u001b[0m\n\u001b[0;32m 214\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 215\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m_hook_minimize\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mloss\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mparams\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 216\u001b[1;33m \u001b[1;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_call_minimize\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mloss\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mloss\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mparams\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mparams\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 217\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 218\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m_call_minimize\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mloss\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mparams\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32m~\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\minimizers\\baseminimizer.py\u001b[0m in \u001b[0;36m_call_minimize\u001b[1;34m(self, loss, params)\u001b[0m\n\u001b[0;32m 218\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m_call_minimize\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mloss\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mparams\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 219\u001b[0m \u001b[1;32mtry\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 220\u001b[1;33m \u001b[1;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_minimize\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mloss\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mloss\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mparams\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mparams\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 221\u001b[0m \u001b[1;32mexcept\u001b[0m \u001b[0mNotImplementedError\u001b[0m \u001b[1;32mas\u001b[0m \u001b[0merror\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 222\u001b[0m \u001b[1;32mtry\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32m~\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\minimizers\\minimizer_minuit.py\u001b[0m in \u001b[0;36m_minimize\u001b[1;34m(self, loss, params)\u001b[0m\n\u001b[0;32m 136\u001b[0m minimizer_setter)\n\u001b[0;32m 137\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_minuit_minimizer\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mminimizer\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 138\u001b[1;33m \u001b[0mresult\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mminimizer\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmigrad\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m**\u001b[0m\u001b[0mminimize_options\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 139\u001b[0m \u001b[0mparams_result\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[0mp_dict\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mp_dict\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mresult\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 140\u001b[0m \u001b[0mresult_vals\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[0mres\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m\"value\"\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mres\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mparams_result\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32miminuit\\_libiminuit.pyx\u001b[0m in \u001b[0;36miminuit._libiminuit.Minuit.migrad\u001b[1;34m()\u001b[0m\n",
"\u001b[1;31mRuntimeError\u001b[0m: exception was raised in user function\nUser function arguments:\n bplus_0 = -0.627540\n p4160_s = +3.373292\n DDstar_s = +0.299698\n p3770_s = +3.780920\n jpsi_p = -5.063578\n p3770_p = +2.196764\n bplus_1 = +1.062520\n DDstar_p = +4.628551\n omega_s = +5.708895\n phi_s = +16.739100\n Ctt = +3.195816\n p4415_s = +0.911410\n psi2s_p = +2.482877\n Dbar_p = +0.982372\n p4040_p = -2.632544\n p4160_p = +2.741969\n p4415_p = -2.983702\n phi_p = -4.941643\n rho_p = +1.587532\n rho_s = +0.149543\n bplus_2 = -0.691416\n Dbar_s = +0.150965\n omega_p = -0.692164\n p4040_s = +1.222875\nOriginal python exception in user function:\nKeyboardInterrupt: \n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\minimizers\\minimizer_minuit.py\", line 101, in func\n loss_evaluated = self.sess.run(loss_val)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\client\\session.py\", line 929, in run\n run_metadata_ptr)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\client\\session.py\", line 1152, in _run\n feed_dict_tensor, options, run_metadata)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\client\\session.py\", line 1328, in _do_run\n run_metadata)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\client\\session.py\", line 1334, in _do_call\n return fn(*args)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\client\\session.py\", line 1319, in _run_fn\n options, feed_dict, fetch_list, target_list, run_metadata)\n File \"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\client\\session.py\", line 1407, in _call_tf_sessionrun\n run_metadata)\n"
]
}
],
"source": [
"# zfit.run.numeric_checks = False \n",
"\n",
"load = False\n",
"\n",
"bo5 = True\n",
"\n",
"bo5_set = 5\n",
"\n",
"fitting_range = 'cut'\n",
"total_BR = 1.7e-10 + 4.9e-10 + 2.5e-9 + 6.02e-5 + 4.97e-6 + 1.38e-9 + 4.2e-10 + 2.6e-9 + 6.1e-10 + 4.37e-7\n",
"cut_BR = 1.0 - (6.02e-5 + 4.97e-6)/total_BR\n",
"\n",
"Ctt_list = []\n",
"Ctt_error_list = []\n",
"\n",
"nr_of_toys = 1\n",
"nevents = int(pdg[\"number_of_decays\"]*cut_BR)\n",
"# nevents = pdg[\"number_of_decays\"]\n",
"event_stack = 1000000\n",
"# nevents *= 41\n",
"# zfit.settings.set_verbosity(10)\n",
"\n",
"mi = 1e-4\n",
"ma = 1e-2\n",
"ste = 10\n",
"\n",
"# mi = 0.0\n",
"# ma = 1e-3\n",
"# ste = 11\n",
"\n",
"BR_steps = np.linspace(mi, ma, ste)\n",
"\n",
"Ctt_steps = np.sqrt(BR_steps/4.2*1000)\n",
"\n",
"print(Ctt_steps)\n",
"\n",
"total_samp = []\n",
"\n",
"start = time.time()\n",
"\n",
"Nll_list = []\n",
"\n",
"sampler = total_f.create_sampler(n=nevents)\n",
"\n",
"__ = 0\n",
"\n",
"#-----------------------------------------------------\n",
"\n",
"if not load:\n",
"\n",
" for Ctt_step in Ctt_steps:\n",
" \n",
" __ += 1\n",
" \n",
" newset = True\n",
" \n",
" for floaty in [True, False]:\n",
"\n",
" Ctt.floating = floaty\n",
"\n",
" Nll_list.append([])\n",
" Ctt_list.append([])\n",
" Ctt_error_list.append([])\n",
" \n",
" if bo5:\n",
" \n",
" if __ >= 6:\n",
" \n",
" while len(Nll_list[-1])/bo5_set < nr_of_toys:\n",
"\n",
" print('Step: {0}/{1}'.format(__, ste))\n",
"\n",
" print('Current Ctt: {0}'.format(Ctt_step))\n",
" print('Ctt floating: {0}'.format(floaty))\n",
"\n",
" print('Toy {0}/{1} - Fit {2}/{3}'.format(int(len(Nll_list[-1])/bo5_set), nr_of_toys, len(Nll_list[-1])%bo5_set, bo5_set))\n",
"\n",
" reset_param_values()\n",
"\n",
" if floaty:\n",
" Ctt.set_value(Ctt_step)\n",
" else:\n",
" Ctt.set_value(0.0)\n",
" \n",
" if newset or len(Nll_list[-1])%bo5_set == 0:\n",
" sampler.resample(n=nevents)\n",
" s = sampler.unstack_x()\n",
" total_samp = zfit.run(s)\n",
" calls = 0\n",
" c = 1\n",
" newset = False\n",
"\n",
"\n",
" data = zfit.data.Data.from_numpy(array=total_samp[:int(nevents)], obs=obs_fit)\n",
"\n",
" ### Fit data\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",
"\n",
" if result.converged:\n",
" Nll_list[-1].append(result.fmin)\n",
" Ctt_list[-1].append(params[Ctt]['value'])\n",
" Ctt_error_list[-1].append(params[Ctt]['minuit_hesse']['error'])\n",
"\n",
" else:\n",
"\n",
" while len(Nll_list[-1]) < nr_of_toys:\n",
"\n",
" print('Step: {0}/{1}'.format(__, ste))\n",
"\n",
" print('Current Ctt: {0}'.format(Ctt_step))\n",
" print('Ctt floating: {0}'.format(floaty))\n",
"\n",
" print('Toy {0}/{1}'.format(len(Nll_list[-1]), nr_of_toys))\n",
"\n",
" reset_param_values()\n",
"\n",
" if floaty:\n",
" Ctt.set_value(Ctt_step)\n",
" else:\n",
" Ctt.set_value(0.0)\n",
"\n",
" if floaty:\n",
" sampler.resample(n=nevents)\n",
" s = sampler.unstack_x()\n",
" total_samp = zfit.run(s)\n",
" calls = 0\n",
" c = 1\n",
"\n",
"\n",
" data = zfit.data.Data.from_numpy(array=total_samp[:int(nevents)], obs=obs_fit)\n",
"\n",
" ### Fit data\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",
" Nll_list[-1].append(result.fmin)\n",
"\n",
" _t = int(time.time()-start)\n",
"\n",
" print('Time Taken: {}'.format(display_time(int(_t))))\n",
"\n",
" print('Predicted time left: {}'.format(display_time(int((_t/(__+1)*(ste-__-1))))))\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"if load:\n",
" Nll_list = []\n",
" CLs_values = []\n",
"\n",
" _dir = 'data/CLs/finished/f1d1'\n",
" \n",
" jobs = os.listdir(_dir)\n",
" \n",
" for s in range(ste):\n",
" CLs_values.append([])\n",
" \n",
" for s in range(2*ste):\n",
" Nll_list.append([])\n",
" \n",
" for job in jobs:\n",
" if not os.path.exists(\"{}/{}/data/CLs/{}-{}_{}s--CLs_Nll_list.pkl\".format(_dir, job, mi,ma,ste)):\n",
" print(job)\n",
" continue\n",
" \n",
" with open(r\"{}/{}/data/CLs/{}-{}_{}s--CLs_Nll_list.pkl\".format(_dir, job, mi,ma,ste), \"rb\") as input_file:\n",
" _Nll_list = pkl.load(input_file)\n",
" \n",
"# print(_Nll_list)\n",
" \n",
" if bo5: \n",
" for s in range(2*ste):\n",
" if _Nll_list[s]: \n",
" Nll_list[s].append(np.min(_Nll_list[s]))\n",
" print(np.min(_Nll_list[s]))\n",
" else:\n",
" print('empty in {}'.format(s))\n",
" else:\n",
" for s in range(2*ste):\n",
" Nll_list[s].extend(_Nll_list[s])\n",
" \n",
"# with open(r\"{}/{}/data/CLs/{}-{}_{}s--CLs_list.pkl\".format(_dir, job, mi,ma,ste), \"rb\") as input_file:\n",
"# _CLs_values = pkl.load(input_file)\n",
" \n",
"# for s in range(ste):\n",
"# CLs_values[s].extend(_CLs_values[s])\n",
" \n",
" print(np.shape(Nll_list))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dirName = 'data/CLs'\n",
"\n",
"# if bo5 and not load:\n",
"# for s in range(2*ste):\n",
"# Nll_list[s] = [np.min(Nll_list[s])]\n",
"\n",
"if bo5 and load: \n",
" CLs_values= []\n",
" for i in range(int(len(Nll_list)/2)):\n",
" CLs_values.append([])\n",
" for j in range(len(Nll_list[2*i])):\n",
" CLs_values[i].append(Nll_list[2*i][j]-Nll_list[2*i+1][j])\n",
"\n",
"\n",
"if not load:\n",
" \n",
" if not os.path.exists(dirName):\n",
" os.mkdir(dirName)\n",
" print(\"Directory \" , dirName , \" Created \")\n",
"\n",
" with open(\"{}/{}-{}_{}s--CLs_Nll_list.pkl\".format(dirName, mi,ma,ste), \"wb\") as f:\n",
" pkl.dump(Nll_list, f, pkl.HIGHEST_PROTOCOL)\n",
" \n",
" with open(\"{}/{}-{}_{}s--Ctt_list.pkl\".format(dirName, mi,ma,ste), \"wb\") as f:\n",
" pkl.dump(Ctt_list, f, pkl.HIGHEST_PROTOCOL)\n",
" \n",
" with open(\"{}/{}-{}_{}s--Ctt_error_list.pkl\".format(dirName, mi,ma,ste), \"wb\") as f:\n",
" pkl.dump(Ctt_error_list, f, pkl.HIGHEST_PROTOCOL)\n",
" \n",
"# CLs_values = []\n",
" \n",
"# for i in range(int(len(Nll_list)/2)):\n",
"# CLs_values.append([])\n",
"# for j in range(nr_of_toys):\n",
"# CLs_values[i].append(Nll_list[2*i][j]-Nll_list[2*i+1][j])\n",
"\n",
"# with open(\"{}/{}-{}_{}s--CLs_list.pkl\".format(dirName, mi,ma,ste), \"wb\") as f:\n",
"# pkl.dump(CLs_values, f, pkl.HIGHEST_PROTOCOL)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# # print(CLs_values)\n",
"# for Nll_l in Nll_list:\n",
"# if Nll_l:\n",
"# print(np.min(Nll_l))\n",
"# print(Nll_l)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Plot"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# l = []\n",
"\n",
"# if not os.path.exists('data/CLs/plots'):\n",
"# os.mkdir('data/CLs/plots')\n",
"# print(\"Directory \" , 'data/CLs/plots' , \" Created \")\n",
"\n",
"# for i in range(len(CLs_values)):\n",
"# plt.clf()\n",
"# plt.title('Ctt value: {:.2f}'.format(Ctt_steps[i]))\n",
"# plt.hist(CLs_values[0], bins = 100, range = (-25, 25), label = 'Ctt fixed to 0')\n",
"# plt.hist(CLs_values[i], bins = 100, range = (-25, 25), label = 'Ctt floating')\n",
"# plt.axvline(x=np.mean(CLs_values[0]),color='red', linewidth=1.0, linestyle = 'dotted')\n",
"# plt.legend()\n",
"# plt.savefig('data/CLs/plots/CLs-BR({:.1E}).png'.format(BR_steps[i]))\n",
" \n",
"# l.append(len(np.where(np.array(CLs_values[i]) < np.mean(CLs_values[0]))[0]))\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"# for s in range(len(l)):\n",
"# print('BR: {:.4f}'.format(BR_steps[s]))\n",
"# print(2*l[s]/len(CLs_values[s]))\n",
"# print()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# for i in range(len(Nll_list)):\n",
"# print(np.mean(np.array(Nll_list[i])))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
