{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Import"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\util\\execution.py:57: UserWarning: Not running on Linux. Determining available cpus for thread can failand be overestimated. Workaround (only if too many cpus are used):`zfit.run.set_n_cpu(your_cpu_number)`\n",
" warnings.warn(\"Not running on Linux. Determining available cpus for thread can fail\"\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"WARNING: The TensorFlow contrib module will not be included in TensorFlow 2.0.\n",
"For more information, please see:\n",
" * https://github.com/tensorflow/community/blob/master/rfcs/20180907-contrib-sunset.md\n",
" * https://github.com/tensorflow/addons\n",
"If you depend on functionality not listed there, please file an issue.\n",
"\n"
]
}
],
"source": [
"import os\n",
"\n",
"# os.environ[\"CUDA_VISIBLE_DEVICES\"] = \"-1\"\n",
"\n",
"import numpy as np\n",
"from pdg_const import pdg\n",
"import matplotlib\n",
"import matplotlib.pyplot as plt\n",
"import pickle as pkl\n",
"import sys\n",
"import time\n",
"from helperfunctions import display_time, prepare_plot\n",
"import cmath as c\n",
"import scipy.integrate as integrate\n",
"from scipy.optimize import fminbound\n",
"from array import array as arr\n",
"import collections\n",
"from itertools import compress\n",
"import tensorflow as tf\n",
"import zfit\n",
"from zfit import ztf\n",
"# from IPython.display import clear_output\n",
"import os\n",
"import tensorflow_probability as tfp\n",
"tfd = tfp.distributions"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# chunksize = 10000\n",
"# zfit.run.chunking.active = True\n",
"# zfit.run.chunking.max_n_points = chunksize"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Build model and graphs\n",
"## Create graphs"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"def formfactor(q2, subscript, b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2): #returns real value\n",
" #check if subscript is viable\n",
"\n",
" if subscript != \"0\" and subscript != \"+\" and subscript != \"T\":\n",
" raise ValueError('Wrong subscript entered, choose either 0, + or T')\n",
"\n",
" #get constants\n",
"\n",
" mK = ztf.constant(pdg['Ks_M'])\n",
" mbstar0 = ztf.constant(pdg[\"mbstar0\"])\n",
" mbstar = ztf.constant(pdg[\"mbstar\"])\n",
"\n",
"\n",
" mmu = ztf.constant(pdg['muon_M'])\n",
" mb = ztf.constant(pdg['bquark_M'])\n",
" ms = ztf.constant(pdg['squark_M'])\n",
" mB = ztf.constant(pdg['Bplus_M'])\n",
"\n",
" #N comes from derivation in paper\n",
"\n",
" N = 3\n",
"\n",
" #some helperfunctions\n",
"\n",
" tpos = (mB - mK)**2\n",
" tzero = (mB + mK)*(ztf.sqrt(mB)-ztf.sqrt(mK))**2\n",
"\n",
" z_oben = ztf.sqrt(tpos - q2) - ztf.sqrt(tpos - tzero)\n",
" z_unten = ztf.sqrt(tpos - q2) + ztf.sqrt(tpos - tzero)\n",
" z = tf.divide(z_oben, z_unten)\n",
"\n",
" #calculate f0\n",
"\n",
" if subscript == \"0\":\n",
" prefactor = 1/(1 - q2/(mbstar0**2))\n",
" _sum = 0\n",
" b0 = [b0_0, b0_1, b0_2]\n",
"\n",
" for i in range(N):\n",
" _sum += b0[i]*(tf.pow(z,i))\n",
"\n",
" return ztf.to_complex(prefactor * _sum)\n",
"\n",
" #calculate f+ or fT\n",
"\n",
" else:\n",
" prefactor = 1/(1 - q2/(mbstar**2))\n",
" _sum = 0\n",
"\n",
" if subscript == \"T\":\n",
" bT = [bT_0, bT_1, bT_2]\n",
" for i in range(N):\n",
" _sum += bT[i] * (tf.pow(z, i) - ((-1)**(i-N)) * (i/N) * tf.pow(z, N))\n",
" else:\n",
" bplus = [bplus_0, bplus_1, bplus_2]\n",
" for i in range(N):\n",
" _sum += bplus[i] * (tf.pow(z, i) - ((-1)**(i-N)) * (i/N) * tf.pow(z, N))\n",
"\n",
" return ztf.to_complex(prefactor * _sum)\n",
"\n",
"def resonance(q, _mass, width, phase, scale):\n",
"\n",
" q2 = tf.pow(q, 2)\n",
"\n",
" mmu = ztf.constant(pdg['muon_M'])\n",
"\n",
" p = 0.5 * ztf.sqrt(q2 - 4*(mmu**2))\n",
"\n",
" p0 = 0.5 * ztf.sqrt(_mass**2 - 4*mmu**2)\n",
"\n",
" gamma_j = tf.divide(p, q) * _mass * width / p0\n",
"\n",
" #Calculate the resonance\n",
"\n",
" _top = tf.complex(_mass * width, ztf.constant(0.0))\n",
"\n",
" _bottom = tf.complex(_mass**2 - q2, -_mass*gamma_j)\n",
"\n",
" com = _top/_bottom\n",
"\n",
" #Rotate by the phase\n",
"\n",
" r = ztf.to_complex(scale*tf.abs(com))\n",
"\n",
" _phase = tf.angle(com)\n",
"\n",
" _phase += phase\n",
"\n",
" com = r * tf.exp(tf.complex(ztf.constant(0.0), _phase))\n",
"\n",
" return com\n",
"\n",
"\n",
"def axiv_nonres(q, b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2):\n",
"\n",
" GF = ztf.constant(pdg['GF'])\n",
" alpha_ew = ztf.constant(pdg['alpha_ew'])\n",
" Vtb = ztf.constant(pdg['Vtb'])\n",
" Vts = ztf.constant(pdg['Vts'])\n",
" C10eff = ztf.constant(pdg['C10eff'])\n",
"\n",
" mmu = ztf.constant(pdg['muon_M'])\n",
" mb = ztf.constant(pdg['bquark_M'])\n",
" ms = ztf.constant(pdg['squark_M'])\n",
" mK = ztf.constant(pdg['Ks_M'])\n",
" mB = ztf.constant(pdg['Bplus_M'])\n",
"\n",
" q2 = tf.pow(q, 2)\n",
"\n",
" #Some helperfunctions\n",
"\n",
" beta = 1. - 4. * mmu**2. / q2\n",
"\n",
" kabs = ztf.sqrt(mB**2. + tf.pow(q2, 2)/mB**2. + mK**4./mB**2. - 2. * (mB**2. * mK**2. + mK**2. * q2 + mB**2. * q2) / mB**2.)\n",
"\n",
" #prefactor in front of whole bracket\n",
"\n",
" prefactor1 = GF**2. *alpha_ew**2. * (tf.abs(Vtb*Vts))**2. * kabs * beta / (128. * np.pi**5.)\n",
"\n",
" #left term in bracket\n",
"\n",
" bracket_left = 2./3. * tf.pow(kabs,2) * tf.pow(beta,2) * tf.pow(tf.abs(ztf.to_complex(C10eff)*formfactor(q2, \"+\", b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2)),2)\n",
"\n",
" #middle term in bracket\n",
"\n",
" _top = 4. * mmu**2. * (mB**2. - mK**2.) * (mB**2. - mK**2.)\n",
"\n",
" _under = q2 * mB**2.\n",
"\n",
" bracket_middle = _top/_under *tf.pow(tf.abs(ztf.to_complex(C10eff) * formfactor(q2, \"0\", b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2)), 2)\n",
" \n",
" #Note sqrt(q2) comes from derivation as we use q2 and plot q\n",
"\n",
" return prefactor1 * (bracket_left + bracket_middle) * 2 * q\n",
"\n",
"def vec(q, funcs, b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2):\n",
" \n",
" q2 = tf.pow(q, 2)\n",
"\n",
" GF = ztf.constant(pdg['GF'])\n",
" alpha_ew = ztf.constant(pdg['alpha_ew'])\n",
" Vtb = ztf.constant(pdg['Vtb'])\n",
" Vts = ztf.constant(pdg['Vts'])\n",
" C7eff = ztf.constant(pdg['C7eff'])\n",
"\n",
" mmu = ztf.constant(pdg['muon_M'])\n",
" mb = ztf.constant(pdg['bquark_M'])\n",
" ms = ztf.constant(pdg['squark_M'])\n",
" mK = ztf.constant(pdg['Ks_M'])\n",
" mB = ztf.constant(pdg['Bplus_M'])\n",
"\n",
" #Some helperfunctions\n",
"\n",
" beta = 1. - 4. * mmu**2. / q2\n",
"\n",
" kabs = ztf.sqrt(mB**2. + tf.pow(q2, 2)/mB**2. + mK**4./mB**2. - 2 * (mB**2 * mK**2 + mK**2 * q2 + mB**2 * q2) / mB**2)\n",
" \n",
" #prefactor in front of whole bracket\n",
"\n",
" prefactor1 = GF**2. *alpha_ew**2. * (tf.abs(Vtb*Vts))**2 * kabs * beta / (128. * np.pi**5.)\n",
"\n",
" #right term in bracket\n",
"\n",
" prefactor2 = tf.pow(kabs,2) * (1. - 1./3. * beta)\n",
"\n",
" abs_bracket = tf.pow(tf.abs(c9eff(q, funcs) * formfactor(q2, \"+\", b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2) + ztf.to_complex(2.0 * C7eff * (mb + ms)/(mB + mK)) * formfactor(q2, \"T\", b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2)),2)\n",
"\n",
" bracket_right = prefactor2 * abs_bracket\n",
"\n",
" #Note sqrt(q2) comes from derivation as we use q2 and plot q\n",
"\n",
" return prefactor1 * bracket_right * 2 * q\n",
"\n",
"def c9eff(q, funcs):\n",
"\n",
" C9eff_nr = ztf.to_complex(ztf.constant(pdg['C9eff']))\n",
"\n",
" c9 = C9eff_nr + funcs\n",
"\n",
" return c9"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def G(y):\n",
" \n",
" def inner_rect_bracket(q):\n",
" return tf.log(ztf.to_complex((1+tf.sqrt(q))/(1-tf.sqrt(q)))-tf.complex(ztf.constant(0), -1*ztf.constant(np.pi))) \n",
" \n",
" def inner_right(q):\n",
" return ztf.to_complex(2 * tf.atan(1/tf.sqrt(tf.math.real(-q))))\n",
" \n",
" big_bracket = tf.where(tf.math.real(y) > ztf.constant(0.0), inner_rect_bracket(y), inner_right(y))\n",
" \n",
" return ztf.to_complex(tf.sqrt(tf.abs(y))) * big_bracket\n",
"\n",
"def h_S(m, q):\n",
" \n",
" return ztf.to_complex(2) - G(ztf.to_complex(1) - ztf.to_complex(4*tf.pow(m, 2)) / ztf.to_complex(tf.pow(q, 2)))\n",
"\n",
"def h_P(m, q):\n",
" \n",
" return ztf.to_complex(2/3) + (ztf.to_complex(1) - ztf.to_complex(4*tf.pow(m, 2)) / ztf.to_complex(tf.pow(q, 2))) * h_S(m,q)\n",
"\n",
"def two_p_ccbar(mD, m_D_bar, m_D_star, q):\n",
" \n",
" \n",
" #Load constants\n",
" nu_D_bar = ztf.to_complex(pdg[\"nu_D_bar\"])\n",
" nu_D = ztf.to_complex(pdg[\"nu_D\"])\n",
" nu_D_star = ztf.to_complex(pdg[\"nu_D_star\"])\n",
" \n",
" phase_D_bar = ztf.to_complex(pdg[\"phase_D_bar\"])\n",
" phase_D = ztf.to_complex(pdg[\"phase_D\"])\n",
" phase_D_star = ztf.to_complex(pdg[\"phase_D_star\"])\n",
" \n",
" #Calculation\n",
" left_part = nu_D_bar * tf.exp(tf.complex(ztf.constant(0.0), phase_D_bar)) * h_S(m_D_bar, q) \n",
" \n",
" right_part_D = nu_D * tf.exp(tf.complex(ztf.constant(0.0), phase_D)) * h_P(m_D, q) \n",
" \n",
" right_part_D_star = nu_D_star * tf.exp(tf.complex(ztf.constant(0.0), phase_D_star)) * h_P(m_D_star, q) \n",
"\n",
" return left_part + right_part_D + right_part_D_star"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Build pdf"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"class total_pdf_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.5, upper_limit=1.5)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Load data"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"x_min = 2*pdg['muon_M']\n",
"x_max = (pdg[\"Bplus_M\"]-pdg[\"Ks_M\"]-0.1)\n",
"\n",
"# # Full spectrum\n",
"\n",
"obs_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/hook_numeric_integrate/mul_1:0' shape=() dtype=float64>"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"total_f_fit.normalization(obs_toy)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Test if graphs actually work and compute values"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"# def total_test_tf(xq):\n",
"\n",
"# def jpsi_res(q):\n",
"# return resonance(q, jpsi_m, jpsi_s, jpsi_p, jpsi_w)\n",
"\n",
"# def psi2s_res(q):\n",
"# return resonance(q, psi2s_m, psi2s_s, psi2s_p, psi2s_w)\n",
"\n",
"# def cusp(q):\n",
"# return bifur_gauss(q, cusp_m, sig_L, sig_R, cusp_s)\n",
"\n",
"# funcs = jpsi_res(xq) + psi2s_res(xq) + cusp(xq)\n",
"\n",
"# vec_f = vec(xq, funcs)\n",
"\n",
"# axiv_nr = axiv_nonres(xq)\n",
"\n",
"# tot = vec_f + axiv_nr\n",
" \n",
"# return tot\n",
"\n",
"# def jpsi_res(q):\n",
"# return resonance(q, jpsi_m, jpsi_s, jpsi_p, jpsi_w)\n",
"\n",
"# calcs = zfit.run(total_test_tf(x_part))\n",
"\n",
"\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",
"probs1 = total_f_fit.pdf(test_q, norm_range=False)\n",
"\n",
"calcs_test1 = 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:13: UserWarning: Creating legend with loc=\"best\" can be slow with large amounts of data.\n",
" del sys.path[0]\n",
"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\IPython\\core\\pylabtools.py:128: UserWarning: Creating legend with loc=\"best\" can be slow with large amounts of data.\n",
" fig.canvas.print_figure(bytes_io, **kw)\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAaUAAAD4CAYAAABMtfkzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOy9e3zcVZ3//3x/5pLJPWnSe9omJaX3UtpSUApquZWLXFaQqriorOyyqF9dXYHvrrLrF3YtP5VVkWURVBaVclmVikWUq8By6xXaQml6oU3vTdLmPpmZz/v3x+czySSZW9IkM0nO8/EozJzP+bzPO9P085r3Oe/zPqKqGAwGg8GQDViZdsBgMBgMhihGlAwGg8GQNRhRMhgMBkPWYETJYDAYDFmDESWDwWAwZA3eTDuQbZSXl2tlZWWm3TAYsp7QgS10ePLIHz+9W3vLiTryW/bSMeZU/IH8DHnnEOtja0szeSd20F44jUDhmIz6NRJZv379MVUde7J2jCj1oLKyknXr1mXaDYMh6zn8L9PZU3ImZ371kW7tb/zxV5z5+t+z9xP3M3X+ORnyzuHIv05nV9FSzvraat5e/woLfn8p7537HWYtvy6jfo1EROSDgbBjpu8MBkO/sLBBej9CxOMDwLYjQ+1SAsT5r+V+B88avwzxMKJkMBj6hYWCeHq3e5w2jYSG2qVeCEqXKDl+ZY9YGuJhRMlgMPQLJ1KSXu3iCpUdyfzDX1SjmoTlihKaeb8MiTFrSmkQCoWora2lvb09064YUhAIBKioqMDn82XalRGPYKNxIiXxOI+VbBAlAO0RKWmW+GWIjxGlNKitraWwsJDKykokzjdDQ3agqtTV1VFbW0tVVVWm3RnxONN38daUnMeK2tk1fdc5rWim77IaM32XBu3t7ZSVlRlBynJEhLKyMhPRDhGWxk90sDojEnuoXeqFoJ1TjJ2JDpp5vwyJMaKUJkaQhgfm72noSJx9507f2eGhdikOsZGS46uJlLIbI0oGg6FfJBIly41INJJ5URJA3S8qJtFheGBEaQTx4osvctlllwEQDAY5//zzWbhwIY8++mivvl/96lf5y1/+AjiJHLfeeiszZsxg3rx5LF26lKeffhqAf/u3f+u85/jx49x7770D4mswGOTaa6+lurqaM888kz179sTt98c//pGZM2dSXV3Nd7/73c72lStXsmPHjgHxxdA/LGzU6r0snU0RSeyaElGxzAK/DIkxojRC2bhxI6FQiE2bNnHttdd2u1ZfX8/rr7/OueeeC8C3vvUtDh48yJYtW9iyZQu///3vaWpqAgZPlB588EFKS0upqanha1/7GrfcckuvPpFIhJtvvpmnn36abdu28cgjj7Bt2zYAbrrpJu66664B8cXQP7zYcfcpRTfPZsPDP1aUPG6ig9k8m92kJUoiskJEtotIjYjcGud6jog86l5/Q0QqY67d5rZvF5GLUtkUkSrXxg7Xpj/ZGCJSJiIviEiziNyTwP81IrIlvY8k+9izZw+zZs3i+uuvZ8GCBVx99dW0trYCTiQxa9Ysli1bxm9+8xsAjhw5wnXXXcemTZtYuHAhO3fu7GbviSeeYMWKFQC0trby05/+lB//+Mfk5OQAMH78eD75yU9y66230tbWxsKFC/nMZz7Drbfeys6dO1m4cCH/+I//eFI/05NPPsn1118PwNVXX81zzz1Hz1OQ33zzTaqrq5k+fTp+v5+VK1fy5JNPAnDOOefw7LPPEg5nfopoVKKKJQpWnM2zWTR9B8TsUzIVHYYDKVPCxdkJ9xPgAqAWeEtE1qjqtphuNwANqlotIiuBVcC1IjIHWAnMBSYBz4rIqe49iWyuAu5W1dUicp9r+z8TjQG0A98C5rl/evr/V0Bznz6VJPzr77ey7UDjQJkDYM6kIm7/+NykfbZv386DDz7I2WefzRe+8AXuvfdevvSlL/HFL36R559/nurq6s6IaNy4cTzwwAN873vf46mnnupl69VXX+Xqq68GoKamhqlTp1JUVNSr33e/+13uueceNm3aBDjiuGXLls73PTnnnHM6I6xYvve973H++ed3a9u/fz9TpkwBwOv1UlxcTF1dHeXl5XH7AFRUVPDGG28AYFkW1dXVbN68mcWLFyf+4AyDgh0JYwEab02pM/U686IkMWnrYrnTiib7LqtJJ1JaCtSo6i5V7QBWA1f06HMF8JD7+gngPHHSoK4AVqtqUFV3AzWuvbg23XuWuzZwbV6ZbAxVbVHVV3DEqRsiUgD8A3BHGj9nVjNlyhTOPvtsAK677jpeeeUV3nvvPaqqqpgxYwYiwnXXpVdk8uDBg4wde9LFfHvx8ssvs2nTpl5/egoS0Csqgt6Zc6n6jBs3jgMHDgyA54a+EomWEIq7ppRd03fac/rOJDpkNelsnp0M7It5XwucmaiPqoZF5ARQ5ra/3uPeye7reDbLgOOqGo7TP9EYx5L4/v+A7wOtyX5AEbkRuBFg6tSpybqmjGgGi54P7Oj7/qRA5+bmdu7lqa6uZu/evTQ1NVFYWHhSPvYlUqqoqGDfvn1UVFQQDoc5ceIEY8aMidsnSm1tLZMmTep8397eTm5u7kn5bOgfkXAYH8RdU/JkWaQkPfYpqW0ipWwmnUgp3lOv51fYRH0Gqj1dP7ocElkIVKvqbxP16TSier+qLlHVJYMRQQwEe/fu5bXXXgPgkUceYdmyZcyaNYvdu3d3rhk98sgjyUx0Mnv2bGpqagDIy8vjhhtu4Ctf+QodHR2AE0n98pe/BMDn8xEKOd+KCwsL44pOlL5ESpdffjkPPeQEvk888QTLly/vJbBnnHEGO3bsYPfu3XR0dLB69Wouv/zyzuvvv/8+c+dm5kvCaKczCoqzphTdp5QNazdOpOQ85iyvSXQYDqQjSrXAlJj3FUDPOZPOPiLiBYqB+iT3Jmo/BpS4NnqOlWiMRHwIWCwie4BXgFNF5MWkP2kWM3v2bB566CEWLFhAfX09N910E4FAgPvvv59LL72UZcuWMW3atLRsXXrppbz44oud7++44w7Gjh3LnDlzmDdvHldeeWXn9N6NN97IggUL+MxnPkNZWRlnn3028+bNO+lEhxtuuIG6ujqqq6v5wQ9+0JnufeDAAS655BLAWWu65557uOiii5g9ezaf/OQnO0Xo8OHD5ObmMnHixJPyw9A/ItH6cVa8NaXsSXSwYiIlj6noMDxQ1aR/cKb4dgFVgB/YDMzt0edm4D739UrgMff1XLd/jnv/LsCTzCbwOLDSfX0f8PfJxojx4XPAPQl+hkpgS6qfVVVZvHix9mTbtm292oaS3bt369y5cwfU5tlnn60NDQ0DanMo+cEPfqAPPPBA3GuZ/vsaDZw4dkD19iJ99Vd39rq2f/8+1duLdMNj/5YBz7rTenu5vvKTm1RVtb45qHp7kW78xTcy7NXIBFinaTxjU/1JuaakzvrNl4BnXEH5mapuFZHvuE6sAR4EHhaRGpzoZaV771YReQzYBoSBm1WdVcZ4Nt0hbwFWi8gdwEbXNonGcG3tAYoAv4hcCVyo3bMDDT34/ve/z969eykpKcm0K/2ipKSEz372s5l2Y9RiR1Px46WEe7Mn0cHqdnSFEFbLJDpkOWlVCVfVtcDaHm3fjnndDlyT4N47gTvTsem278LJzuvZnmyMyhT+7yFOuvhwobKyki1bBnab1Zln9sxVGV58/vOfz7QLo5qo4EgcUeqcJsuSRIdoSrglEMEya0pZjqnoYDAY+kzEFZx4ZYYkiyondKvoYAk2lllTynKMKBkMhj6j7vSdxEkJ92bR9J2gnckYlgg2YqbvshwjSgaDoc9EOqfv4mTfebMnJdyKPbpCxJm+M5FSVmNEyWAw9JnOdO8403ceK3um7yzRzlJI0ek7yQK/DIkxojSCGIlHV1RWVjJ//nwWLlzIkiVLOtu/8Y1v8Pzzzw+IL4a+07V5Np4oCSH1ZH6aTLvvu+9MdMi0X4akGFEaoYyEoyuivPDCC2zatIl169Z1tn35y1/udr6SYWix3dp38bPv3ISCTGffRUUppiSXImb6LssxojQMGK1HVyRj2rRp1NXVcejQoZPyw9A/kqaEd67dZDpScsUnppK5WVPKftLap2SI4elb4dA7A2tzwny4OPm3/tF4dAU4324vvPBCRIS//du/5cYbb+y8tmjRIl599VU+8YlPJPzcDIODnUSULEsIZ8N+oE5R6qqpaGMhmRZLQ1KMKA0Teh5d8aMf/Yjzzz+/8+iKaPv999+f0tZgHl2RLvGiongVz1999VUmTZrEkSNHuOCCC5g1a1bntKM5uiJz2JHEa0qQLQ9/53cs9swnJ9HBRErZjBGlvpIiohksRuPRFUDnURXjxo3jqquu4s033+wUJXN0RebQ6JqSp3ekBLhrSpkVJbUjCCAxBwxE8JjpuyzHrCkNE0bj0RUtLS2d47W0tPCnP/2JefO6qkW9//773d4bho7ONaU4m2fBefhnOlKKRuOxkZKazbNZjxGlYcJoPLri8OHDLFu2jNNOO42lS5dy6aWXdiZohEIhampquqWJG4aOqChZCSKlbKgxFz32PPbLji0WYiKlrMZM3w0TLMvivvvu69W+YsUK3nvvvV7tH/3oR/noRz8a19Y555zDbbfdxvHjxykpKcHv93PXXXdx11139eq7atUqVq1a1fn+17/+df9/iBgCgQCPP/54r/ZJkyaxdq1Tp3f69Ols3rw57v1PPfUUV199NV6v+RXOBNE1pXi17wAi4sn4w9+2bRzJ7LGmZCKlrMZESqOU6NEVw5VwOMzXv/71TLsxenH3ICWKlJzCp5ndp9R57HmP7DswkVI2Y75mDgPM0RW9ueaauKeYGIYIOypKCbPvPBkv52PHXVMy03fZjomU0qQvGzsNmcP8PQ0NGomuKcUXpbB4sDQ0lC71JuGakpm+y2aMKKVBIBCgrq7OPPCyHFWlrq6OQCCQaVdGPBE3JTyhKOHDsjMrSrYdr6JD5rMCDckx03dpUFFRQW1tLUePHs20K4YUBAIBKioqMu3GiCcaKXkSJJrY4kEyXPtO40zTRcSDL9YvOxL3SHdD5jCilAY+n4+qqqpMu2EwZA12qpRw8eLNcKJDV0HWLh8jeLCikdL+9fDT5XDDszDljAw4aIiHmb4zGAx9JrpPyZMwJdyLleFIqbM+X8yebFu8WFGx3Po75/+7XxpizwzJSEuURGSFiGwXkRoRuTXO9RwRedS9/oaIVMZcu81t3y4iF6WyKSJVro0drk1/sjFEpExEXhCRZhG5J8ZOnoj8QUTeE5GtImLOOTAYBojoIX+We/R5TyLiw5PhRIeuig4xkVKsKJlpu6wkpSiJU0fkJ8DFwBzgUyIyp0e3G4AGVa0G7gZWuffOAVYCc4EVwL0i4klhcxVwt6rOABpc2wnHANqBbwHfiOP+91R1FnA6cLaIXJzq5zUYDKlRNwryJNqnJJ6uh3+GiO5Tii1e1c0vyxXUSIazBA3dSCdSWgrUqOouVe0AVgNX9OhzBfCQ+/oJ4Dxx8jCvAFaralBVdwM1rr24Nt17lrs2cG1emWwMVW1R1VdwxKkTVW1V1Rfc1x3ABsCsgBsMA0C0ooMnQfadLb7Mi1LnmlJMRQfxdq0peVxRynCWoKE76YjSZGBfzPtaty1uH1UNAyeAsiT3JmovA467NnqOlWiMlIhICfBx4LkE128UkXUiss5k2BkMqelcU0owfaeWF0+mC7JGN+9226cUZ/rOREpZRTqiFO9shJ4bdhL1Gaj2dP3ohYh4gUeAH6nqrnh9VPV+VV2iqksG45whg2GkoSnKDEUsH55siZSsmIoO4ukSy+g+pkwf227oRjqiVAtMiXlfAfQ8Wa2zjysCxUB9knsTtR8DSlwbPcdKNEYq7gd2qOp/pNHXYDCkg5vo4EsUKYk3C0Qpuk8pZvrOivErKkaRjqF1zJCUdETpLWCGmxXnx0lcWNOjzxrgevf11cDz6nxNWQOsdDPnqoAZwJuJbLr3vODawLX5ZIoxEiIid+CI11fT+DkNBkOadEZKvviiZFtePGRHokOv6TuikZI7bRdqx5A9pNw8q6phEfkS8AzgAX6mqltF5DvAOlVdAzwIPCwiNTjRy0r33q0i8hiwDQgDN6s6sXM8m+6QtwCrXUHZ6Nom0RiurT1AEeAXkSuBC4FG4J+A94ANbv2re1T1gb5/TAaDIRZx12G8Xn/c62r5smbzrMQWZI2N4KJrSWEjStlEWhUdVHUtsLZH27djXrcDccs2q+qdwJ3p2HTbd+Fk5/VsTzZGZQLX+35WuMFgSE00JdwXX5QQLx4yXSU8TqRkdUVK7cEOAkCwvYWcDPhniI+p6GAwGPqOHcZWwZsoJdzjw5s103cxjznLi9dNdNh5uAGAg8cahto1QxKMKBkMhr5jhwjhwUowF6GWN+OiFD3MT3qsKUXXumw3WcMTCQ69a4aEGFEyGAx9RuwwYTzdHvjdsHx4sbuKomYAO+KIknaLlHx4XLHyutN4lm1EKZswomQwGPqM2CHCyZaks6CEjxAn0cHyOKJk250JD14TKWUVRpQMCXnjdz+Bfymmva0l064Ysg03UkqEZkEJn85D/mLynTRa1dwOdSY8eEyklFUYUTIkpHLT9wE4fuxghj0xZBuWHSYiSapsR0UpgxtTo5tnrdiFr05RCiPumpJXzebZbMKIksFg6DOpp++ca3Y4c5FS/Oy7mGnF6PSdiZSyCiNKBoOhz4gmn74Tj7N/KRTK3ANfo6UxY5IxNJrCbkc6IyWfEaWswoiSwWDoM2JHsJNM34n78I+EMliB2+59HLrErCmJewihTzsymiVo6I4RJYPB0GfEDhFJNn3nRkqRTEZK8Y7OiJm+E/doCwvbHF+RRRhRMhgMfcbSMBFJLEriVg8PhTKYRBCtfRdz7Hk0K1DtMBKbGRhuG1LXDIkxomQwGPqMI0pJpu8sN1IKZ06UorXvYjf4WlZ0WrEDiY2kTKXwrMGIksFg6DOWHU46fWd5M59917mmFFuXObrWFe7ofly7qRSeNRhRMqTGLAIbepB6+s6pux3O4JqSHY2EYlLCOxMwImHENqKUjRhRMhgMfcbSMHZSUXKm7+xwBtOto/uUYo5Dj04r2iEnUoqoG0WFzJpStmBEyZCaREU3DaMWT8pIKRcAuyNzEYho782z4k4rht3pu2YcP8mkeBq6YUTJkBozfWfogWgEtRInOnj8AQAiHZmLQOxoyrcVI56eaATXgaWRGFEykVK2YETJkBA1B/caEuDRMJokUvL63Ugpg1lt6oqSxE7feZy1LjsUxNIwLeqKZ7B16B00xMWIkiEh0dL/BkNPLA1jW8lEKRopZYMoxUR03qgotWNphBY3UuoImkgpWzCiZDAY+oxHI6j4El735WRTpNQlnuJzRCkSCuLRCM1upBQ2kVLWkJYoicgKEdkuIjUicmuc6zki8qh7/Q0RqYy5dpvbvl1ELkplU0SqXBs7XJv+ZGOISJmIvCAizSJyTw+/FovIO+49P5KEx2Qa4mGm7wyJ8BDu3PMTD68rSpoVotT1mLN8jghpuL1booMRpewhpSiJiAf4CXAxMAf4lIjM6dHtBqBBVauBu4FV7r1zgJXAXGAFcK+IeFLYXAXcraozgAbXdsIxgHbgW8A34rj/n8CNwAz3z4pUP6+hCzN9Z0iEV8OolThS8kcjpQzu/4k3fedxRcnuaMfCphUTKWUb6URKS4EaVd2lqh3AauCKHn2uAB5yXz8BnOdGJVcAq1U1qKq7gRrXXlyb7j3LXRu4Nq9MNoaqtqjqKzji1ImITASKVPU1VVXgv2NsGQyGk8BDpOvAvDj4A25WWyYLsrr7lGJFqXOtKxTEo2Ga1PEzk2tfhu6kI0qTgX0x72vdtrh9VDUMnADKktybqL0MOO7a6DlWojGS+V2bwm8ARORGEVknIuuOHj2axOToIjp9Fz3B02CI4iVFpOT1EFQfmtFIyXmMWJ7eoqShdjxEaJPMp64bupOOKMVbWOg5r5Ooz0C1p+tHOj71blS9X1WXqOqSsWPHJjE5uohO34mYfBhDd3waRt09P/HI8XkI4sto+Z7o9B2xkVKOO30XDuIlgnj8BNWLPRQVHewIHNw8+OMMc9J52tQCU2LeVwAHEvURES9QDNQnuTdR+zGgxLXRc6xEYyTzuyKF3waDoa/YNj6JoO6en3j4PRZBvEgkg5US3AjfihEln7umRIezhuTx+gjiH5qEjDfug/86F/a8OvhjDWPSEaW3gBluVpwfJ3FhTY8+a4Dr3ddXA8+76zhrgJVu5lwVTrLBm4lsuve84NrAtflkijHioqoHgSYROctdq/rrGFsGg6GfdE7JJYmUfB6hA39Gy/doZ0WHLlHy+5xpRQm1AOD1+mjHPzSR0p5XnP/vfW3wxxrGJDk60kFVwyLyJeAZwAP8TFW3ish3gHWqugZ4EHhYRGpwopeV7r1bReQxYBsQBm5W9zjIeDbdIW8BVovIHcBG1zaJxnBt7QGKAL+IXAlcqKrbgJuAXwC5wNPuH4PBcBKEQ0F8AN7EoiQidODDymCk1Jl9F7Om1BnBhdxIyeejvdU3JAVZW+sPkAc0HdpJ4aCPNnxJKUoAqroWWNuj7dsxr9uBaxLceydwZzo23fZdONl5PduTjVGZoH0dMC/eNUP6mEQHQyyhYKsrSomn7wA6xJ9RUSLO5lm/1yKIr1OUvF4/Qfz4h2DtK3R8PwANB3cZUUqCWcE2GAx9Ihx0hEY9gaT9QuJH7MydPBv9MhWbfef3Wk4EF3am73w+H+1DlJCRG25yfGg9POhjDWeMKBlSYoqEG2IJu3t6LF/ySCkkPjxZECl1W1PyWnSoF080UvL5acePNdhrX+EgfnU+t0CoYXDHGuYYUTIYDH0i7GaqSZI1JYCw+LEyGSnFqejgrCn58UScNSTL46MDP1ZkkCOltuMANGgBhXZj1wGEhl4YUTKkxoRKhhii1Q8kxZpSxPLjzaAo4R6H7umR6NCBF2/YiZTE46FN8vC603mDRrsjSrt1Ah7szveG3hhRMhgMfSLiRkqWN/maUtjy47WzocxQV6KDZQkhfPgiUVHy0WrlkxNpHlxn2qKiNBGAjqYjgzveMMaIkiElSbaDGUYh4Q5HaCTFmlLECuDVzImSRNeUPN1PyA2JH687fYflpd2TT2DQRclZR9pvOZXOThw7NLjjDWOMKBnSwMx/G7qInpHkSZXoYAXIsTNYZkh7H10BTgKG325zr/lotwoI2K2diRGDQbjVKT7TVlwFQHO9EaVEGFEyGAx9IipKqbLvwt48AprB6tvRNaUe1cyDEsCrIQDE46XDV+BeaBo0V9qbHFHylFc7Q50waeGJMKJkSInaZvrO0EW6a0oRbx45dAxqBJKU6JpSj8MIg1Zu52vx+Ah53K2s7ScGzZUOV5QKJ80EINR0bNDGGu4YUTIkxBzyZ4iH7Z6RFD3yPBHqy3dedAxyZltCBxxR8vRYU4oVJcvjIeR3RSnYOGiuhFqP06wBpo4vp1kD2C1GlBJhRMmQEjXiZIghmhLuy0keKeHLc/4fysyprl0FWbs/5jqsvK43nhwIFDmv2wdPlCJtjTSTy7iiAMcpRFqTHXAwujGiZEhI9JA/TO07Qwy2W/3AnyJSwp/ZSEk0QlgtLKv70WohT4zfvhysQInzehAjJW0/QZPmUZzro8lTgi9YN2hjDXeMKBkSYqbvDPGIuNN3OYHkomTlOKIUahvkdOsEqNpEsLCkuyiFvV2RkngDePOjkdLgrSlJsIlmcinO9dHqLSbQYUoNJcKIkiElZp+SIZbogXipRcnJamtvHbyHfVLsCIrg6SFKGiNKeAN488sACDUPXvRidTTRpLkU5XoJ+kvJi2ToMxkGGFEyGAx9QsNBbBUCOclTwj0BR5Q6WjMTKYlGnEipx1NOo9OKOKWScgrLCKtFxyCmaXtCzbRaeeR4PUQCYyi2jSglwoiSISUmUDLEIqFWWskh4Et+HJs3Kkptg7f/Jyl2/Ok7/AWdLy1fgOK8HOopIjyIpX/84WaCHmdcO6+MAB3YwQxlJWY5RpQMqTGJDoYYJNxGOzm9Egh64st1Uq3D7ZmJlNAIETy9pu8kNlLyBSgKeDmmxWjz0UFzJSfSQtjriJKVPxaAJlPVIS5GlAwGQ5+wwm20k/zYCgB/VJQyFCmJHSaMRc9AiUDXua8ef4CiXB/HtAhpHSRRsiMEtI2IWznCX1QOQGOdEaV4GFEypMbM3xlisCLtBCX5ehJATp7z8I9kaJrK0hBhvEgPVYrkje98Lb5cigI+jlGMt22QNrS6qeaa43weuSXO+M0NptRQPIwoGQyGPuENt9EhKTbOAoE8dw0lmKFEBztMGE/vC262HYDPEopzfdRpEf7B2jsUramX46Se55c6ohQ8YY6viEdaoiQiK0Rku4jUiMitca7niMij7vU3RKQy5tptbvt2EbkolU0RqXJt7HBt+k9ijK+JyFYR2SIij4ik8S/J0AtT0cEQi8dup8NK/U8pL8dHk+YOaqHTZIi7ptST2KzBHJ+H0jwfdRTji7TDYAioWylCch1RKi5zzlQKNQ3eGtZwJqUoiYgH+AlwMTAH+JSIzOnR7QagQVWrgbuBVe69c4CVwFxgBXCviHhS2FwF3K2qM4AG13Z/xpgMfAVYoqrzAI/bz2AwnAS+SBshK/X0Xa7fQyN5WMHMpD+LHYorSnk+D5/puI2/Cv4LAZ+F12PR4neSD2g8MOB+hNucn9+bWwxA6ZhywmqZ+ncJSCdSWgrUqOouVe0AVgNX9OhzBfCQ+/oJ4DxxJnKvAFaralBVdwM1rr24Nt17lrs2cG1e2c8xALxAroh4gTxg4H/jRgFm86whFp/dTiiNSCnf76VR87A6MhMpWXaYsPROW8/ze3jVns8GPZUcryNarXmTnIsn9g64H62NTvUGX75Tzsjj8XBcirBaTamheKQjSpOBfTHva922uH1UNQycAMqS3JuovQw47troOVafxlDV/cD3gL3AQeCEqv4p3g8oIjeKyDoRWXf0qAmpe2FSwg0x+OwgYU8aa0o+i0by8XVkrqJDvEgp19/VFvA5j8BwYRl/MGsAACAASURBVIXTcHzgRamt2RGlQH5xZ1uTVYwvaIqyxiMdUYq3GaHnV+dEfQaqvc9jiEgpThRVBUwC8kXkujh9UdX7VXWJqi4ZO3ZsvC4Gg8HFr0EiVopirICI0GYV4OsYvEKnybA0hC29Rakw4Ot8HY2UvMWTnKSI4/t69T9ZOpqciCinqOvZYurfJSYdUaoFpsS8r6D3NFhnH3eqrBioT3JvovZjQIlro+dYfR3jfGC3qh5V1RDwG+DDafy8hh6Y2TtDLDnaTiTFAX9R2ryF+MMZmr7TMBF6T9+VF3TtsYpGSuVFeRzQMnQQIqWwW1Mvr7i8s609p4zCsImU4pGOKL0FzHCz4vw4yQJrevRZA1zvvr4aeF6dhYg1wEo3c64KmAG8mcime88Lrg1cm0/2c4y9wFkikueuPZ0HvJvex2IwGOKi6pwm60sdKQF0eIvIjWRu82wkzppSWUFXkkae37leXpBDrV1OpOGDAfdDW+to1DxKCroKwQZzJ1CmdeYbXxxSipK7fvMl4Bmch/pjqrpVRL4jIpe73R4EykSkBvgH4Fb33q3AY8A24I/AzaoaSWTTtXUL8A+urTLXdn/GeAMnIWID8I77s97fj89o1GMSHQydhIN4sLtOlU3V3V9ErrZCJJy68wBjaSTu9F1+zJqSxy2VNLEkwB4dD3U1Ay8UrQ00aAEleV3ThuGCiQQIETTHovcieUVFF1VdC6zt0fbtmNftwDUJ7r0TuDMdm277Lrqy52Lb+zPG7cDt8e4x9AGT6GBw0VArAlj+9CKliN89qyjYCHljBs+xOFgaJiJ5vdp7VngAmFKax++0Am/7C9ByFArGDZwf7fU0UMDk3C5RkiInf6vxyAeMLTLr2LGkJUoGg8EAEGxtJACQU5iqKwCa457q2n48I6Jkx5m+A3jqy8u6vZ8yJo/31c3AO/LugIqSr+M4TVKM19M1MeUf4yyDNx/5gLHVSwZsrJGAKTNkSImZvjNEaW9yMsasQFFa/cXdMErb8cFyKSGJpu8A5k0uZt7krhTt0jwf+71TnTdH3xtQP3JCJ2j1dv+8SiZWAtBydOATK4Y7RpQMCTHHoRt6Emx2xMWTpihZbnQUbh76tZNkkVJPRIRA6WRarAI4sm3gnFClMFxPq697lDhh0jTnYMH6gU9BH+4YUTKkgREng0PQPdrck1ecoqeDVeQUH20/PvTHNHj6IEoAU8vyeV+q4MDGgXOi/Tg5GqQ10H06sDg/wFEpHZSyRsMdI0qGhKi7L9lM3xmihFocUfKlKUo+V5Q6TmRIlKz0RWnG+AJe6zgFPbQFOlr7NtieV+HVH/bOMmw8CEBH3oRet9R7x5PXWtu3cUYBRpQMCTHTd4aeRNyK1zn56YlSfmExbeondGLozw6ysNEEa0rxmDWhiHWRakQjfYuWOlrgl38Ff/42vHFft0vauN/xpbhnZTY4kTuF8o796Y8zSjCiZEiNiZQMLhG34nWOW1w0FWUFAY5pMXbT0J8d5NUQtuVL3dFl1oRCNtnVzpvaN9MfaPfLEG53Xr/xX2B3baFoO+YkMgTGVPS6LVhUSbnWoxk62iNbMaJkMBjSRtsbsVXIK0gvUiovzOEYxUjL0IuST0OE0zhiI0pVeT7NnhKOBSph91/SH2j/ehALrrjXqTK++6XOS+2HthNUH/njpva6zSp3BLB+3/b0xxoFGFEypMSsKRk6aW+kmVzyctJbqynL93NUi/G2D3H2nSp+Ooj0QZS8HotZEwt5w3O6s0aU7rpS/U5acifz95ursHNKYOPDXW4c3c4uncj44t4VMIonzwLgyJ4BzPYbARhRMiTErCkZetHRTBO53SptJyPg83DCKiUwWEeNJyISwkKJWP7UfWNYPK2U3zTOgkgQ9ryS1j32sRo2tIxh7XsNvBhYDu/+HtwD/AL177FDJzO5pHcFjMmnzAWg9ZCJlGIxomRIAyNOBgcJNtFKLn5v+o+OFn8ZeeHjEO4YRM964K7x2J70IyWAMyrH8EpoJrYnAO8/nfoGVexjNdREJnBaRTGrjp4JkQ7Y9Gto+ID89kNslllMKOpdVb18zBgOMwap29EnH0c6RpQMCelMCbdN7TuDg6/jBM1Wehtno7TkTsJCoXEIM83CQaDvorSkspQgfnaXnQPbnoRIKPkNzUfwhlvYoxP4j5Wns5Op7Ck43UkPX/9zAPaVLMWyetfbExH2+6czpun9Pvk40jGiZEiImb4z9CQQOt6rZE4q2gvcdOhBOKsoIdFsuD6K0rjCANXjCngycja01sGuF5PfUL8TgKa8qVSV53PJ/Il8o2kl2n4CXrmbN2Q+BZPnJLy9sWQ2k8N70VB7n/wcyRhRMqREjTgZXPIiJ2j3ppcOHkWL3MyzIRUlJ1JSX3qHEcZy3uxx/PTgKWigBDavTt65rsYZp8zJpPubc6pYF5zC70+/n7ozvsHft/09SyoTF6K1Js7HR4Sju9/ps58jFSNKBoMhPVQptBvp8PdNlPLKpxJRIVS3Z3D8iocbKUmaJ+TGcv7s8bTZHvZMuhTeXQPNidPZtW4nITwUjp8GwIKKEs6sGsPtGwv4l8bLqJdiPjIz8dEU4049A4CD2/uwL2qEY0TJkBqTEm4ACDbhI0wkt7RPt00sK+QQY2g7unuQHOuNnoQoLZpaSlm+n4fCFzhJC+t/kbBvx5Ed7LXHMW1sl1DfedU8gmGb328+wOWnTYqbeRdl+qkLaNEcQrUDWG9vmGNEyZAaI0oGgLZ6ADTQt3ORJpfkUatj0YY9g+BUfELBNgDE17c1JXBOo71i4WR+tTOHUNVyeOvBzunAnthHa9ilE5le3rUPqXpcIU//n3P44cqFrPrEgqRj+X1edvpnUlq/qc9+jlSMKBkMhrQIuUd3S35Zn+6bXJpLjT2ZwPEdQ/YFJ9ThiJLVjzUlgGuWVBCKKM+XfhKaD8GG/+7dybbxN+5ht06kqrz75thpZflcsXAyAV/q2nsNZYup7Kgh7FZgH+0YUTKkxkRKBqClwVlb8RaU9+m+8YU51FBBTqgRmoemMGukvRlI/9j2nsyeWMSCimLuen8COuVD8PL3IdTWvVPjfjx2kL1MoKK0f+MA5Ez/MB5RPnj75X7bGEkYUTIYDGnR3OAcP5FX0rejwr0ei2O5pzhvjrw70G7FJdzmFDmVNI9tj8cNy6rYeayVddNvgqaD8NYD3TsccyoxtBRUdTvqvK/MWLyciAp17yWvt2fbSmtHOGmfkUBan6SIrBCR7SJSIyK3xrmeIyKPutffEJHKmGu3ue3bReSiVDZFpMq1scO16T+JMUpE5AkReU9E3hWRD/Xt4zGASQk3OHQ0OJtfC8f2rnidilDZTOfFAB81ngi7PSpKBf22cen8iUwZk8u/v1uOVp8PL66CpphzofY7yQkd4+adlK9lZeXs9lZRcPD1pP1W//J+HrvzevYdG9lVxVOKkoh4gJ8AFwNzgE+JSM/dYDcADapaDdwNrHLvnQOsBOYCK4B7RcSTwuYq4G5VnQE0uLb7PIZ7zw+BP6rqLOA0YGi+po0QonvQTUFWA0Ck8SCNmktZad/WlADGT6ygTovRg5sHwbPeRI+DkEB61czj4fVY/O25p7Bh73FeO/WbTj28Z/6pa4zat9itE5kwbvxJ+3u4/MNUt28h2HI87vVIxOaqnf/M5+QpNr/w2EmPl82kEyktBWpUdZeqdgCrgSt69LkCeMh9/QRwnoiI275aVYOquhuoce3Ftenes9y1gWvzyv6MISJFwLnAgwCq2qGq8f/GDQZDSjxNBzmsYygv7FuRU4Dq8YWst6uJ7B2a/Th2exNhtfD1c00pyrVnTGF6eT7ferkN++yvwZYn4O3HoaMV3f0XXo7Mo7K8dwXwvhKYfSF+ibDrzT/GvX5sfw254tQO9O958aTHy2bSEaXJwL6Y97VuW9w+qhoGTgBlSe5N1F4GHHdt9Byrr2NMB44CPxeRjSLygIjE/e0RkRtFZJ2IrDt69GjiT2K0YiIlA+BvO8IxKSXPn/4R41FmjCtgvX0q3oad0DIEFcM7mmkhgD+N7Ldk+DwWt1w8i51HW/hVzjUw9cOw5svw+Oewwm38PvIhTh3X/ynCKLOWXkCL5tC8Lb4o1e91Jnk68DK1eRPhyMitR5mOKPWuJNi7bHSiPgPV3p8xvMAi4D9V9XSgBei1Hgagqver6hJVXTJ2bOLd16MNI0WGWPKCRzjh69+/D0eUZjhv9r0xgF4lINhMM7nk9KGaeSIunDOes6vLWPWnXRxc8V8w8TTY8Qxbx32cdcxi3uT+TxFGyc/L4/28RUw+9mrcAshth51yRjvGraCaWnbsH7lfntP5G6sFpsS8rwAOJOojIl6gGKhPcm+i9mNAiWuj51j9GaNWVaP/Ap7AESlDHzFrSgZsm6JwHS3+vmXeRSkryOFg/mw6JAd2vTDAzvVGg820aID8NA8jTIaI8N2/WoCtyi1/PIT9uafhlj38u//LzBhXMCBjAISnn88kPcLu99b3umbX7aRdfYxZ+HG8YrPn3XUDMmY2ko4ovQXMcLPi/DhJBWt69FkDXO++vhp4Xp0n2RpgpZs5VwXMAN5MZNO95wXXBq7NJ/szhqoeAvaJiJv2w3mAOeLRYOgPzYfwEiFcMKnfJuZOHcdb1gLY/sdBnxL2tNfTQOGACcaUMXncdvEs/vL+Ue5+bgd1kTze2F3Hx2b2T6TjUbXsGmwVjrz+eK9rOU0fcMCawIRZZwLQsmfDgI2bbaT8G1PVsIh8CXgG8AA/U9WtIvIdYJ2qrsFJJnhYRGpwopeV7r1bReQxHDEIAzeragQgnk13yFuA1SJyB7DRtU1/xgC+DPzKFb5dwOf79SmNdkykNOqx63ZhAVpa2W8bp08t5ffbT+PsyANweCtMOLlU6mT4gvXU6xiq/Ce3phTLdWdNY8v+Rn78fA2Pr6slbCvXLOl7enwiyidOY5tvNuP3PwPc1e1aftshjvsmIKWVtFgFBI5tGbBxs420vkao6lpgbY+2b8e8bgeuSXDvncCd6dh023fhZOf1bO/PGJuAJfHuMaSPMnIXVQ3p0XRwB8VAYHx1v20smlrCg5HF2P5fYL29GibcMXAO9iAnWE+9Vg5YpATONN6dV81jQnGA5947zNcumEH1uP5vzo1HY9XFzNnxffbueJupM7rq5pWGj3CodA6IUF84k4qGHTS1h9I+ln44YSo6GAyGlLQeep+QeiiZOL3fNk6bUkKTt5T3ipc5x4UnKHJ60tgRAqHjHKOI3JPMvuuJ12PxtQtO5akvn8O1Z0wdUNsA08/9NAB7X+k6x6mjtYlSGrELnahMJp7GLNnL23uHIIsxAxhRMqTGTN+NeiJ1u9mv5VSU9T8yCPg8fOiUMn7e/lHnVNfNjwycg7G01iMozVZx3GPIs5lxU6qp8Z3KuH1/7EwwOrLfOd3WV+ac2TRmxpkEJMTeOAkRIwEjSoaEmOPQDVH8J3bzgY5nUpKzgdLhYzPH8fjxGbSPXwwv3dW7yOlAcMI54bbOO2HgbQ8BjdVXcqq9kx3vvAVA/YFdABRNqAQgr9I5GLDjg7cy4t9gY0TJkBoTKY1uImFKW3ay31/Vr42zsZw3exwgPDX2i9C4H17otRR88rjnNjXkDE9ROuW8zxNSD8de/TkAre4epfIKd5/XmOm0WoUU1r8zIrdrGFEyJETdfckj8Pfe0BfqavBpiKaSman7pqCiNI+zq8v44a7x6KLPwf/eA+8/c/I+xlLvnHDbmj8lRcfspLh8ElsLzuLUw3+gPRjEOvY+LRqgbEKV00GE46XzmRXZQW3DIESaGcaIkiEhZvrOAGAfctKPZfzApHBfs3gK++rbeHn6V2HiAnjsetg5gBtqD73DIWs8eflFA2dziMk547OUc4L1f/o1+ce3s983FSvmeAzvlMWcKvt4e/fBDHo5OBhRMqSBSQkfzTR9sJGQeiidNjCidMn8iUwqDvAffzmAfvpxGFMFv/wE/OX/g3BH34ypwsG3u8J5Vahdx1aqKc3re+HYbGHWsk+w35rAlA2rmBHcSt2Y07tdH3PqWXjF5sj7Q1PgdigxomQwGJIS+eB1tuo0qieOGRB7fq/FTR+rZsPe4zy7D/jCMzDncnj+DvjJGfDGf0HLsfSM7XwO/uscePHfnfcHNkBjLS+G51GaP3xFSbx+Gs78JlP1IH6JULqk+xZN7xRn+6XsH3nlhgZuZ5lh5GIWlUYv4SBFdW+zQc/n0xMHbjrs2iVT+NXrH/Ct323hzH84l6JrfgELPwMvrYKnvwlP3wIVS2DyYpgwH0qmQeFEyBsDlgfsCBzbAS//wDH40ioINsPe19CcIp48sZgb84b3xtJ5F93AFlXCKixcekH3i4UTqM+pYErjRjrCNv4BKDybLRhRMiQkuqZkJGkUc2AjXu3gUMnpBAZwI6rfa/HdTyzgE//5v3xt9Sbu/+sleGZcADMugIOb4b21TuHWDf8NodbEhjw5sGIVHNgIr98L3gB1F/yQxt8WMK4wMGD+Zop5K/4m4bXmCWexZM9a3t3fwGnT+n7wYrZiRMmQGhMpjVrsXS+BCta0Dw247YVTSrj943P49pNb+affvsOdV83HY4lzNMTE0+BjtzkRUcMeOFHrHEXeVg9qAwIlU2Hah53oCeCSu8DysmNfO/A6k0tPbk9VtlMw6yMUf/AEH7z7FqdNW5FpdwYMI0qGhERTwo0ojV7atvyB7XoKc2acMij2//pDlRxtCvLj52uoa+nge1efRnHstJvlgbJTnD+pcI8+33+8AYDJJ7nRN9spnf1ReAbad7wEK0aOKI2ciUjDgGNSwkc5jQfJP7aZZ+3FLKsuH7Rhvn7hTG7/+BxeeO8IF//wLzz9zsGT2hS6t74VS2BiyfCfvkuGlEylzjeRsrq3RtRJtEaUDAmJVg1TI06jk22/A+CD8o8yZpAz2T5/dhX/c9OHKQz4uOlXG7j6vtd4ZushInbff/e2HWhk+tgCcrwDW4w1G2mdeBaL9F3e3lefaVcGDCNKhiQYMRq1qBJe9xCb7elUzx2a019Om1LCH76yjDuvmsehE+387cPr+dj3XuRHz+2gtiFJskMMqsqW/SeYM4CZgtlM6fyLKJVm3t/4cqZdGTCMKBkS0jl9Z9aURh/71+M99i6PRj7G5Qv7f9psX/F6LD5z5jRe+sePcu9nFjGpJMAP/vw+y1a9wMr7X+PxdftobA8lvL/mSDOHGttZWjUwe6qynYI5F2AjUPNcpl0ZMEyigyEhXWtKI2e+2pAe+vL3aZYCasav4JSxBUM+vtdjccn8iVwyfyL76lv57cb9/M+GWv7xibf5p99u4SMzx3LZgomcP3t8t4P8Hl9fiyVw/uzxQ+5zRsgv50j+TKqb3qSxPUTRCDj0z4iSISHD6yQaw4Bx6B1k+1oeCH2Cy5eefBHWk2XKmDy+ct4Mvry8mo37jvPU5oP84Z0D/HnbYQI+i/NmjWf5rHHUt3Tws1d2c+XCyUwoHtlJDrFEpi9n4dv38dJ7ezhv4YxMu3PSGFEyJETcCMnM3o0iVGHtN2m2Cvmt/zKeWVSRaY86EREWTS1l0dRS/vnS2by1p56n3j7I2ncO8od3nMKkS6vGcPvlczPs6dAy7vRL8b5zL4c2PgNGlAwjGVFMuDTa2PQr2Pu//L/QF7l6+Xxy/dmZwWZZwpnTyzhzehm3f3wOe+pa8FoW08ryEBldv7S+aWfSJnkU1L6E6s3D/udPK9FBRFaIyHYRqRGRW+NczxGRR93rb4hIZcy129z27SJyUSqbIlLl2tjh2vT3dwz3mkdENorIU+l/LAYwiQ6jjqPvo2u/yRbvPF7Ku5C/Oacq0x6lhddjUT2ukMry/GH/QO4XHh9Hx57F4vAGdh9tzrQ3J01KURIRD/AT4GJgDvApEZnTo9sNQIOqVgN3A6vce+cAK4G5wArgXlckktlcBdytqjOABtd2n8eI8e3/AO+m93EYYjGiNIporYdHr6MdHzc0/x23XTr3pE+ZNQwdBXMupEKOsXnT8D8iPZ1IaSlQo6q7VLUDWA1c0aPPFcBD7usngPPE+cpyBbBaVYOquhuoce3Ftenes9y1gWvzyn6OgYhUAJcCD6T3cRgMo5BgM/z6k9j1u/li65c4Y8FcLj9t6NLADSfPmNMuBqDtvT9n2JOTJx1Rmgzsi3lf67bF7aOqYeAEUJbk3kTtZcBx10bPsfo6BsB/AN8kRU6ziNwoIutEZN3Ro0eTdR1VWJ2JDiYlfMTScgwe+ji6fz1f169woPQM7rxy/uicBhvOlFZyNGcqFXX/SzAcybQ3J0U6ohTvt7PnfE6iPgPV3ucxROQy4Iiqro9zvXtn1ftVdYmqLhk7dmyq7qMG81ga4RzeBg9egH14K9/0/CMvWWfx88+f0b0gqmHY0Db1YyxlKxtqhvcR6emIUi0wJeZ9BXAgUR8R8QLFQH2SexO1HwNKXBs9x+rrGGcDl4vIHpzpweUi8ss0fl6DiynIOkJRhfUPwU8/Rqj1BF+wv8WLspRHvngW08ryM+2doZ+MO/1SAhJiz4bhPYWXjii9Bcxws+L8OEkFa3r0WQNc776+GnhenTK/a4CVbuZcFTADeDORTfeeF1wbuDaf7M8YqnqbqlaoaqVr/3lVvS7Nz8UAmNp3I5Dje+GRlfD7r1BbuIBljXewN38+T/zdh5g5oTDT3hlOgsCMc+nAT84Hz2falZMiZXqNqoZF5EvAM4AH+JmqbhWR7wDrVHUN8CDwsIjU4EQvK917t4rIY8A2IAzcrKoRgHg23SFvAVaLyB3ARtc2/RnDcHJY0ZNnTfbd8CcchDfugxe/i63wcOEX+deDH2HF/Ems+sQCCkdAeZpRjy+Xw2MWs+DYOg6daB+2VS3EPHC6s2TJEl23bl2m3cgKIreX4BFl44d+xOkXXZ/6BkP2YUfgncfhhTvh+F5qxpzLDYc/SYNvHP906Ww+uWSKSWoYQRx65gdMeO1feWr5M1x27llDOraIrFfVky4pb6qEGxJiHlXDGNuGd5+C+86B3/4tdZE8brb+mfMP/B0L5s3j2a9/hGvPmGoEaYQxftGlADRt+WOGPek/ZnecISGWmOm7YUe4w4mMXv0hHNtOY+4Uvuf7Og8fPZ0zp5fzu4tns3BKSaa9NAwSUn4q9b4JjD/yChFb8VjD70uHESWDYSTQ3ggbH4bXfgKN+zmSP4Pvy1d5omExiyrL+fnV1Xzk1LEmMhrpiNA4+SOcsXsN73xwhIVVw+8IDyNKhtSYSCl7OfIevPVTdPNqpKOZmryFfDfyWZ6tm88Fcybw2Eems3ja6DjwzuBQvvASCvY8yo4Nz7Ow6lOZdqfPGFEypMRoUpYRCcP2P8CbP4U9LxOx/LzgXcaPgh9jp30qVy2ezLMfrqR6nEnxHo0UzFpOGA/WzucAI0oGg2GwOFELm36Nrvs50nSABt8Efm5/mofbz2XixAquO28aly+cREGO+Wc9qgkUcbDoNGYdf5MTraFhV6HD/PYa0sDUvssY4SBsXwsbHkZ3Po+grPcs5L6OT/G/kcVcPL+Cn501lYVTSsx6kaETmXE+c9ffxbPvbOP8M0/LtDt9woiSwZCNHHoHNv4S3fwo0t7AMc9YfhW+iscj5zKpchZXL6rg7vkTzKZXQ1wmLLoM1t9F3dtPgxElw4jDLCoNDW0N8M4T6MZfIgc3ERYff7KX8OvQR6ktXsJVH5rGI4smM2VMXqY9NWQ53kkLOOEZQ+nBl1G9ZVhF0UaUDClRUwNv8Ai1w45n0LcfQ9//E5bdwXYq+XXoep7zfoRzFszg/yyuYMm00mH1YDFkGBHqJy7jjH3P8/7BE8ycNHz2phlRMhiGGtuGD15F334Me+vv8HQ0UkcJvwufx1Ocy+TZZ/Hx0ybyf2eOI+DzpLZnMMShZMEllNau4aX1LzFzUs9zWbMXI0qG1Jjpu4Hh8FZ4+zFCmx/D13yANgI8HTmD3+syfNUf5bKFU/jl7PEme84wIJTOuwh7rRDZ8Wd6HxaevZjffoNhMDm+D7b+huCGR8ipe5cIFn+JnMYa+69onnYBF51+Cj+cO2HYpe0ahgF5YziYP5uq46/THooMm6jbiJIhNSZS6huNB9Ctv6Vt0/+Qd9g5+HirXc3vItezf9IKPrJoDv88byJjC3My7KhhpBOeei7zt/2U9TtrOWvWtEy7kxZGlAyGgaDpELr1d7RsfIKCw28hwB57GmvtT7J34goWLVzETfMmMLE4N9OeGkYRYxecj+/d+zjw9gsw63OZdictjCgZ4mJHIl3nmphIKT7NR7C3PknzhscpPPwmgrLPnsLT9jUcrFjB6acv5XNzx1NeYCIiQ2bIO+VsQnjxfPAy8LlMu5MWRpQMcQmHQ/g73xlR6qTlGJGtT9K04XGKDr2Bhc1hezK/0E9wdOrFnLboLL4wexwlef7UtgyGwcafx8HCeUw/sWHYrCsZUTLEJRLqyLQL2cPxfQS3rKF5028pPbYeDzZ19kR+zVXUTbuEBYvO4vOzx5vqCoasJDL1HOZsuYd1NR9w5uzpmXYnJUaURjjrf3M3xdMXU73w3D7dFwqHiK5+jMpD/o5up3nTbwm+8yRljdvIAXbbU3jCuoqmqkuYt+jDfH7meHL92f/N0zC6GX/a+Xi2/piDm58DI0qGTLP47X+Bt4GFJ/p0XyjYOij+ZC2q6IGN1K/7H2T7U4xp3UMBsMOu5rGcv8aeeRmLFy3hhmmleD1WSnMGQ7aQN/1DBPHj2/cq8MVMu5OStERJRFYAPwQ8wAOq+t0e13OA/wYWA3XAtaq6x712G3ADEAG+oqrPJLMpIlXAamAMsAH4rKp29HUMEZni9p+AU+b6flX9YV8/oOGM2jb9LUwTbG2JMTRCIyU7QnjP/3LszSfI3fU0xR2HKVaL1+3ZbC26Cf+8j/Ph0xfwd+MLK7GwnAAAEYlJREFUTIkfw/DFm8OBogWccmL9sFhXSilKIuIBfgJcANQCb4nIGlXdFtPtBqBBVatFZCWwCrhWROYAK4G5wCTgWRE51b0nkc1VwN2qulpE7nNt/2c/xggDX1fVDSJSCKwXkT/38HtEE2xvJdDPezvammLejSBRCjbT9t6fqdvwJCW1z1MQOUGJ+nhFF7Cr/K8pWfhxzj1tJsuK+/vJGQzZhz3tHGa9czdvvr+TpXNPTX1DBkknUloK1KjqLgARWY1TsyL24X4F8C/u6yeAe8T5ankFsFpVg8BuEalx7RHPpoi8CywHPu32eci1+599HUNVXwMOAqhqk2t7cg+/RzRtLU39F6X2ltSdhgl6fC/HNqyhfcsfmFD/FrmEKNQ8XpJFHJ18AROXXMbZcyo535T3MYxQJiy8EN65m8NvPwcjQJQmA/ti3tcCZybqo6phETkBlLntr/e4d7L7Op7NMuC4qobj9O/PGACISCVwOvBGvB9QRG4EbgSYOnVqvC7DkraWRkr7eW84RpSGXaKDbdOxbx0H3/wdObv+xIS2HYwFdtkTeDJwCcHpF1K95HwurBqHz6wPGUYB+ZVn0Cq5+Pe9AtyUaXeSko4oxZtM7/mUStQnUXu8J0Gy/v0Zw7lJpAD4H+CrqtoYpy+qej9wP8CSJUuG2RM4McHWuD9uWoTbmwfQkyGgo5WGLX+ibuMayg+8QEmkngoVNjCTl0q/SGDeZSxatJSrzVlEhtGIx8fBwvlUnnibYDhCjjd715XSEaVaYErM+wrgQII+tSLiBYqB+hT3xms/BpSIiNeNlmL793kMEfHhCNKvVPU3afysI4pgS98y7mIJtcUKWnbqdOT4fmrf+C3hd9dScfwtSunAq7m86VlE/dTljF98GWfMnsEZJm3bYMCuOIsZJ37M27v2svDUqky7k5B0ROktYIabFbcfJ6ng0z36rAGuB14DrgaeV1UVkTXAr0XkBzhJCDOAN3Gim1423XtecG2sdm0+2Z8x3PWmB4F3VfUHff1gRgIdriiF1ENft3VGmo52vpZsmb6zIzTvepP9bz1J3gfPMaX9faYBe3Ucz+VfTLh6BTPPvIjlk8aYbDmDoQdj5y3H2vYjDm15AYazKLnrN18CnsFJ3/6Zqm4Vke8A61R1Dc7D/2E3yaAeR2Rw+z2Gk1wQBm5W1QhAPJvukLcAq0XkDmCja5u+jiEiy4DPAu+IyCbXxv9V1bX9+6iGHyF3+i6Iv8+iZLfUDbxD/UBbG9i/7ilatj7NxCOvUKQnqFZhs8xk89gvkr/g4yxa9CEuyTdlfQyGZJTMOIsQXvjgNeALmXYnIWmlG7kP8rU92r4d87oduCbBvXcCd6Zj023fRVeGXmx7n8ZQ1VeIv940aoi0OZFSUPwU9PFeq61LlIY00UGV5n1vs//NJ/Htepapre9QgU29FrA+Z8n/397dB1V13gkc//64XBSMCgIaBBUV5EURwTd8iU1jqiaNtmNMaptOskmnzWyyM012u3lZd7e17UybpLPd7Uy6mbTNNNlto01bE2u0hiTGrOIb4gsggkB8BUQRUEBA4Nk/zmMhegkvuXgvl99n5s4997nPec5zfnruj3Puc59D46SlxM1dSXpCPJlBQ/qfV6m+cYdyLjSZ8fWHMMb47dUEHQMbwDqanaTUSt/PIoY3nqOO2whn4Ac8mJYGzh3aTv3RrYyr2klUxwWSgGNMJjviG4SkriBt3l18cfSIAe+LUoHs6vh5pJS+TnnFRabGRvu6Ox5pUgpkDc73Qq1Bfb91QkTzWaqC4whvO85ADHRorCrl9N5NBJVlE38ljziuEW6GcyQkg9z4bzNu9irSkpNI1SHbSnlNePKduMte49TRnUyNXePr7nikSSmABTdWARBkOvq03tWGemI6qsgbNRPqj3ulL6athXNHdnDp8F+IrNxJbNsZUoBPTAy7wlfhSl5O8rwVLIoc7ZXtKaVudvuMJXS8K1wr34UzXsz/aFIKYCOuOqPpXbT1UPPTyvM+YLp0YCYuhPy/9nvuu2uNtZTveZvWwi1Mrs0hjiaiTTAF7jSOTVpDVMZKZqRlMFnPhpS6JSQ0grMhU4iqOejrrnRLk1KA6mhvZ0JrOQgEm74lpea8jTSa4cTM/CLk922gw+XzpziZ8xbuE38loSmPJNqpMaM4NHIJHdPuYVrWfcweG9XX3VFKecnlsXNIOfM21XVXGBs+0tfduYkmpQBVvH87KXKVyyYMl7T3er1zJw6TVpvNwciVJIwY1at1Gi9VUrrjfwkr2URiSyEzgVPEsDvqQcLSVpI2/0vcMVyHbCvlD8ISlxB2diOHj+Qw9gvLfd2dm2hSCkCFu95h5IfPc4mRlEbcQVLdzl6tV3/hHG1vfpOrDCfhgR91/mjWw3dSpqOdkt2baNv7K5Ia9pMuHZTJBD6KfZzoeQ+QPGM2k/SynFJ+JzZ9KeyAxpKPQZOSGkjlR3fTtPVfmdGcRyXRnL3rZToKtxJsej5TKs7ZTHj204zrqKfkS79lZsxE6mvOA2DaOm+NbtrbOLrtVaLzfkFSRyXVRJAz7huMWfAQqekLmKq/HVLKr4WEx1DpGs+oCwd83RWPNCkFgIryY1RtWkfmlQ+pZSQ5if9E5urvERMaxp6i9wim+6RUdvBDWt5bT2rLYc4xjvJVm5g5+w4A3MPsjS/amp3tnDhM88bHSG8rozgogZyMl8hY9jBLQvXeQ0oNJhfHzGZa9Q6aWloJG+Zfl9Y1KQ1i50+XcGrzT8i48A7hBJMT9yjTH/g3FoZHdlYKchPsYfTdyYI9XN76fWY27aOG0exO+B4Zq58mNqxz7oeQYaGAc6Z0Ijeb8Vu+yXDc7Ml4kXn3fZskvTyn1KDknrKIiAt/4VB+LhlzFvq6O5+iSWkQOl2cR/W2F0mvfY8I4GDkfUy5/0csjJ10c2WXG5cY2tvbcblcVJYXULVpHRlXPqLejGBX/JOk3/8Mi0aF37RqsDuENhOE+8pporf8HTUyhuBHt7Bg4tSB30ml1ICJnbUU9kFt0U7QpKT6o7npCoXZrxNW8HtSrhUSbULIHbua+JXPkjUxsdv1JMiZirXlagP5G39A5unXGUUwu+MeY/rqf2Fx5GdPNdKKmzm1W2k2bmrXvslETUhKDXojb0/kokQSVpHj667cRJOSH2tuaqA4ZzPXCt5hWt3HzKaJMzKePVO+S+Lyx1kwLrbnRlzOP3HJL7/G/KY97Bu9nMlrf8ai8b27w26YtABwKHoVC5Jn9XtflFJ+RIRzEXNIqNlL67V2Qtz+c88xTUp+puFyLcW7/oQc20zylb2kSwuXGUHx6DsInf8IqfPvYUJfvstxOWdKs5r2kDPpCRY++pN+9St2+VP9Wk8p5Z+CJi8h6tJ2CgoOMCMjy9fd+RtNSn6g7kIlJ3a9RUjJu6Q05TJb2qhhNPmRKwhN/yrJWfcyd1j/RrgF3+bMntBs3Mx/+Md9Xn9/2g9or69iQeLMfm1fKeWfJmYuh4PrqCl4HzQpqQsVJyn/vw2MKN9GcvNR5koHlUSTN24NozJXkzRnKfODP/8/T9ysuzlzaDx1S39Kmqvvp+jz7n/6c/dBKeV/Rscmcl7GElaxx9dd+RRNSrdQ5aliTu3aSPjJbUxrLSJaDKeDYjkQ9zBRc9eQMHMRMUHeHWZ9+8RE+H4RE7zaqlIqEFREzCGx5mNar7UR4vaPdOAfvQhgVxuvUJD9BmHH3mR6az4xQJlrMvviH2d81oNMSplN74YcKKWUd7mmfoHwS1s5lr+X1MzFvu4OoElpwNRdrOL42y+SenYDc2nkrMSwN/4J4hY/xNSEGejAaqWUr02cvQIOPMulgmzQpBSY2tvayP3jS6QW/RdZcpVDIxYTsuhJUrNWEOflS3NKKfV5hN8ez0nXJCLOvA+s93V3AOjVp6SIrBCRYhEpFZHnPLw/TEQ22vf3iUh8l/eet+XFIrK8pzZFZLJt44RtM8Tb2xgo9TXnOfbS3cw//lM+CU3l5Nc+IOOf32X6wnsRTUhKKT90IW4ZKa2FVFec9nVXgF4kJRFxAS8D9wCpwNdFJPWGat8Cao0xCcDPgRfsuqnAWmA6sAL4pYi4emjzBeDnxphEoNa27e1teF3dxSpqX15KUnM++9PWk/bM+8SnzBmozSmllFeMm/8gQWIo++gNX3cF6N2Z0jyg1BhTboxpBTYAX7mhzleA1+3yH4GlIiK2fIMxpsUY8wlQatvz2KZd5y7bBrbNr3pzG70LS9+0tjRT/cpKYtqrKFn2W+bd/5SeGSmlBoWJKXMpck8npeS/qTpT6uvu9Oo7pVjgTJfXZ4H53dUxxrSJSD0Qacv33rDu9blxPLUZCdQZ87f7d3et761t3EREvgN8x75sEJEa4KKnuj364ap+reanouhvHAKPxqKTxsIReHFY3/08mj2IAjzMCN13vUlKnu7aZnpZp7tyT6cRn1Xfm9u4udCYV4FXr78WkVxjzJC/9qZx6KSx6KSxcGgcOtlYxHujrd5cYzoLn/rtZRxQ0V0dEQkGRgOXPmPd7sovAuG2jRu35a1tKKWU8lO9SUoHgEQ7Ki4EZ1DB5hvqbAYesctrgA+NMcaWr7Uj5yYDicD+7tq06+ywbWDbfMeb2+hdWJRSSvlCj5fv7Pc3/wBsB1zAa8aYQhH5IZBrjNkM/Ab4HxEpxTl7WWvXLRSRPwDHgDbgSWNMO4CnNu0mnwU2iMiPgUO2bby8jZ682nOVIUHj0Elj0Ulj4dA4dPJaLMQ52VBKKaV8T8ctK6WU8hualJRSSvkNTUpd3OppiXxBRF4TkWoRKehSNkZEsu3UTtkiEmHLRUR+YeNxVEQyu6zziK1/QkQe8bQtfyYiE0Rkh4gUiUihiHzXlg/FWAwXkf0icsTGYr0t99qUX4OJnRHmkIhssa+HahxOiki+iBwWkVxbNvDHhzFGH873ai6gDJgChABHgFRf92sA9nMJkAkUdCl7EXjOLj8HvGCX7wW24fwWLAvYZ8vHAOX2OcIuR/h63/oYhxgg0y6PBEpwpqMairEQ4Da77Ab22X38A7DWlr8C/L1dfgJ4xS6vBTba5VR73AwDJtvjyeXr/etHPP4R+D2wxb4eqnE4CUTdUDbgx4eeKXW6ZdMS+ZIx5mOc0YtddZ3C6capnd4wjr04vyGLAZYD2caYS8aYWiAbZ97BQcMYU2mMybPLV4AinJlAhmIsjDGmwb5024fBe1N+DRoiEgd8Gfi1fe3Nqc8CwYAfH5qUOnmaTim2m7qBZpwxphKcD2tgrC3vLiYBFSt72SUD5wxhSMbCXrI6DFTjfHCU0cspv4CuU34N9lj8J/AM0GFf93rqMwIrDuD8YfKeiBwUZyo2uAXHh95PqVNvplMaavo6tdOgIyK3AX8CnjLGXHb+0PVc1UNZwMTCOL/tmyUi4cAmIMVTNfsckLEQkfuAamPMQRG583qxh6oBHYcuFhljKkRkLJAtIsc/o67XYqFnSp2G8rRE5+2pNva52pYH9BROIuLGSUi/M8b82RYPyVhcZ4ypAz7C+V7AW1N+DRaLgFUichLn8v1dOGdOQy0OABhjKuxzNc4fKvO4BceHJqVOQ3laoq5TON04tdPDdmRNFlBvT9m3A8tEJMKOvllmywYNe+3/N0CRMeY/urw1FGMRbc+QEJFQ4G6c79i8NeXXoGCMed4YE2eciUXX4uzXQwyxOACIyAgRGXl9Gef/dQG34vjw9QgPf3rgjCApwbmevs7X/RmgfXwTqASu4fwV8y2c6+AfACfs8xhbV3BulFgG5ANzurTzGM4XuKXAo77er37EYTHOZYSjwGH7uHeIxmImzpReR+0Hz7/b8ik4H6alwFvAMFs+3L4ute9P6dLWOhujYuAeX+/b54jJnXSOvhtycbD7fMQ+Cq9/Ht6K40OnGVJKKeU39PKdUkopv6FJSSmllN/QpKSUUspvaFJSSinlNzQpKaWU8hualJRSSvkNTUpKKaX8xv8DM7Nwv4XCekcAAAAASUVORK5CYII=\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, 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": null,
"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": null,
"metadata": {},
"outputs": [],
"source": [
"# 0.213/(0.00133+0.213+0.015)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Adjust scaling of different parts"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"total_f.update_integration_options(draws_per_dim=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": null,
"metadata": {},
"outputs": [],
"source": [
"# # print(\"jpsi:\", inte_fl)\n",
"# # print(\"Increase am by factor:\", np.sqrt(pdg[\"jpsi_BR\"]/pdg[\"NR_BR\"]*pdg[\"NR_auc\"]/inte_fl))\n",
"# # print(\"New amp:\", pdg[\"jpsi\"][3]*np.sqrt(pdg[\"jpsi_BR\"]/pdg[\"NR_BR\"]*pdg[\"NR_auc\"]/inte_fl))\n",
"\n",
"# # print(\"psi2s:\", inte_fl)\n",
"# # print(\"Increase am by factor:\", np.sqrt(pdg[\"psi2s_BR\"]/pdg[\"NR_BR\"]*pdg[\"NR_auc\"]/inte_fl))\n",
"# # print(\"New amp:\", pdg[\"psi2s\"][3]*np.sqrt(pdg[\"psi2s_BR\"]/pdg[\"NR_BR\"]*pdg[\"NR_auc\"]/inte_fl))\n",
"\n",
"# name = \"phi\"\n",
"\n",
"# print(name+\":\", inte_fl)\n",
"# print(\"Increase am by factor:\", np.sqrt(pdg[name+\"_BR\"]/pdg[\"NR_BR\"]*pdg[\"NR_auc\"]/inte_fl))\n",
"# print(\"New amp:\", pdg[name][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": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
" \n",
"# print(list_of_borders[:9])\n",
"# print(list_of_borders[-9:])\n",
"\n",
"\n",
"class UniformSampleAndWeights(zfit.util.execution.SessionHolderMixin):\n",
" def __call__(self, limits, dtype, n_to_produce):\n",
" # n_to_produce = tf.cast(n_to_produce, dtype=tf.int32)\n",
" low, high = limits.limit1d\n",
" low = tf.cast(low, dtype=dtype)\n",
" high = tf.cast(high, dtype=dtype)\n",
"# uniform = tfd.Uniform(low=low, high=high)\n",
"# uniformjpsi = tfd.Uniform(low=tf.constant(3080, dtype=dtype), high=tf.constant(3112, dtype=dtype))\n",
"# uniformpsi2s = tfd.Uniform(low=tf.constant(3670, dtype=dtype), high=tf.constant(3702, dtype=dtype))\n",
"\n",
"# list_of_borders = []\n",
"# _p = []\n",
"# splits = 10\n",
"\n",
"# _ = np.linspace(x_min, x_max, splits)\n",
"\n",
"# for i in range(splits):\n",
"# list_of_borders.append(tf.constant(_[i], dtype=dtype))\n",
"# _p.append(tf.constant(1/splits, dtype=dtype))\n",
" \n",
"# mixture = tfd.MixtureSameFamily(mixture_distribution=tfd.Categorical(probs=_p[:(splits-1)]),\n",
"# components_distribution=tfd.Uniform(low=list_of_borders[:(splits-1)], \n",
"# high=list_of_borders[-(splits-1):]))\n",
" mixture = tfd.MixtureSameFamily(mixture_distribution=tfd.Categorical(probs=[tf.constant(0.05, dtype=dtype),\n",
" tf.constant(0.93, dtype=dtype),\n",
" tf.constant(0.05, dtype=dtype),\n",
" tf.constant(0.065, dtype=dtype),\n",
" tf.constant(0.04, dtype=dtype),\n",
" tf.constant(0.05, dtype=dtype)]),\n",
" components_distribution=tfd.Uniform(low=[tf.constant(x_min, dtype=dtype), \n",
" tf.constant(3090, dtype=dtype),\n",
" tf.constant(3681, dtype=dtype), \n",
" tf.constant(3070, dtype=dtype),\n",
" tf.constant(1000, dtype=dtype),\n",
" tf.constant(3660, dtype=dtype)], \n",
" high=[tf.constant(x_max, dtype=dtype),\n",
" tf.constant(3102, dtype=dtype), \n",
" tf.constant(3691, dtype=dtype),\n",
" tf.constant(3110, dtype=dtype),\n",
" tf.constant(1040, dtype=dtype),\n",
" tf.constant(3710, dtype=dtype)]))\n",
"# dtype = tf.float64\n",
"# mixture = tfd.MixtureSameFamily(mixture_distribution=tfd.Categorical(probs=[tf.constant(0.04, dtype=dtype),\n",
"# tf.constant(0.90, dtype=dtype),\n",
"# tf.constant(0.02, dtype=dtype),\n",
"# tf.constant(0.07, dtype=dtype),\n",
"# tf.constant(0.02, dtype=dtype)]),\n",
"# components_distribution=tfd.Uniform(low=[tf.constant(x_min, dtype=dtype), \n",
"# tf.constant(3089, dtype=dtype),\n",
"# tf.constant(3103, dtype=dtype), \n",
"# tf.constant(3681, dtype=dtype),\n",
"# tf.constant(3691, dtype=dtype)], \n",
"# high=[tf.constant(3089, dtype=dtype),\n",
"# tf.constant(3103, dtype=dtype), \n",
"# tf.constant(3681, dtype=dtype),\n",
"# tf.constant(3691, dtype=dtype), \n",
"# tf.constant(x_max, dtype=dtype)]))\n",
"# mixture = tfd.Uniform(tf.constant(x_min, dtype=dtype), tf.constant(x_max, dtype=dtype))\n",
"# sample = tf.random.uniform((n_to_produce, 1), dtype=dtype)\n",
" sample = mixture.sample((n_to_produce, 1))\n",
"# sample = tf.random.uniform((n_to_produce, 1), dtype=dtype)\n",
" weights = mixture.prob(sample)[:,0]\n",
"# weights = tf.broadcast_to(tf.constant(1., dtype=dtype), shape=(n_to_produce,))\n",
" # sample = tf.expand_dims(sample, axis=-1)\n",
"# print(sample, weights)\n",
" \n",
"# weights = tf.ones(shape=(n_to_produce,), dtype=dtype)\n",
" weights_max = None\n",
" thresholds = tf.random_uniform(shape=(n_to_produce,), dtype=dtype)\n",
" return sample, thresholds, weights, weights_max, n_to_produce"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# total_f._sample_and_weights = UniformSampleAndWeights"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# 0.00133/(0.00133+0.213+0.015)*(x_max-3750)/(x_max-x_min)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# zfit.settings.set_verbosity(10)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"# # zfit.run.numeric_checks = False \n",
"\n",
"# nr_of_toys = 1\n",
"# nevents = int(pdg[\"number_of_decays\"])\n",
"# nevents = pdg[\"number_of_decays\"]\n",
"# event_stack = 1000000\n",
"# # zfit.settings.set_verbosity(10)\n",
"# calls = int(nevents/event_stack + 1)\n",
"\n",
"# total_samp = []\n",
"\n",
"# start = time.time()\n",
"\n",
"# sampler = total_f.create_sampler(n=event_stack)\n",
"\n",
"# for toy in range(nr_of_toys):\n",
" \n",
"# dirName = 'data/zfit_toys/toy_{0}'.format(toy)\n",
" \n",
"# if not os.path.exists(dirName):\n",
"# os.mkdir(dirName)\n",
"# print(\"Directory \" , dirName , \" Created \")\n",
"\n",
"# for call in range(calls):\n",
"\n",
"# sampler.resample(n=event_stack)\n",
"# s = sampler.unstack_x()\n",
"# sam = zfit.run(s)\n",
"# # clear_output(wait=True)\n",
"\n",
"# c = call + 1\n",
" \n",
"# print(\"{0}/{1} of Toy {2}/{3}\".format(c, calls, toy+1, nr_of_toys))\n",
"# print(\"Time taken: {}\".format(display_time(int(time.time() - start))))\n",
"# print(\"Projected time left: {}\".format(display_time(int((time.time() - start)/(c+calls*(toy))*((nr_of_toys-toy)*calls-c)))))\n",
"\n",
"# with open(\"data/zfit_toys/toy_{0}/{1}.pkl\".format(toy, call), \"wb\") as f:\n",
"# pkl.dump(sam, f, pkl.HIGHEST_PROTOCOL)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# with open(r\"data/zfit_toys/toy_0/0.pkl\", \"rb\") as input_file:\n",
"# sam = pkl.load(input_file)\n",
"# print(sam[:10])\n",
"\n",
"# with open(r\"data/zfit_toys/toy_0/1.pkl\", \"rb\") as input_file:\n",
"# sam2 = pkl.load(input_file)\n",
"# print(sam2[:10])\n",
"\n",
"# print(np.sum(sam-sam2))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# print(\"Time to generate full toy: {} s\".format(int(time.time()-start)))\n",
"\n",
"# total_samp = []\n",
"\n",
"# for call in range(calls):\n",
"# with open(r\"data/zfit_toys/toy_0/{}.pkl\".format(call), \"rb\") as input_file:\n",
"# sam = pkl.load(input_file)\n",
"# total_samp = np.append(total_samp, sam)\n",
"\n",
"# total_samp = total_samp.astype('float64')\n",
"\n",
"# data2 = zfit.data.Data.from_numpy(array=total_samp[:int(nevents)], obs=obs)\n",
"\n",
"# data3 = zfit.data.Data.from_numpy(array=total_samp, obs=obs)\n",
"\n",
"# print(total_samp[:nevents].shape)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# plt.clf()\n",
"\n",
"# bins = int((x_max-x_min)/7)\n",
"\n",
"# # calcs = zfit.run(total_test_tf(samp))\n",
"# print(total_samp[:nevents].shape)\n",
"\n",
"# plt.hist(total_samp[:nevents], bins = bins, range = (x_min,x_max), label = 'data')\n",
"# # plt.plot(test_q, calcs_test*nevents , label = 'pdf')\n",
"\n",
"# # plt.plot(sam, calcs, '.')\n",
"# # plt.plot(test_q, calcs_test)\n",
"# # plt.yscale('log')\n",
"# plt.ylim(0, 200)\n",
"# # plt.xlim(3080, 3110)\n",
"\n",
"# plt.legend()\n",
"\n",
"# plt.savefig('test2.png')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# sampler = total_f.create_sampler(n=nevents)\n",
"# nll = zfit.loss.UnbinnedNLL(model=total_f, data=sampler, fit_range = (x_min, x_max))\n",
"\n",
"# # for param in pdf.get_dependents():\n",
"# # param.set_value(initial_value)\n",
"\n",
"# sampler.resample(n=nevents)\n",
"\n",
"# # Randomise initial values\n",
"# # for param in pdf.get_dependents():\n",
"# # param.set_value(random value here)\n",
"\n",
"# # Minimise the NLL\n",
"# minimizer = zfit.minimize.MinuitMinimizer(verbosity = 10)\n",
"# minimum = minimizer.minimize(nll)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# jpsi_width"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# plt.hist(sample, weights=1 / prob(sample))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Fitting"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# start = time.time()\n",
"\n",
"# for param in total_f.get_dependents():\n",
"# param.randomize()\n",
" \n",
"# # for param in total_f.get_dependents():\n",
"# # print(zfit.run(param))\n",
" \n",
"# nll = zfit.loss.UnbinnedNLL(model=total_f, data=data2, fit_range = (x_min, x_max))\n",
"\n",
"# minimizer = zfit.minimize.MinuitMinimizer(verbosity = 5)\n",
"# # minimizer._use_tfgrad = False\n",
"# result = minimizer.minimize(nll)\n",
"\n",
"# # param_errors = result.error()\n",
"\n",
"# # for var, errors in param_errors.items():\n",
"# # print('{}: ^{{+{}}}_{{{}}}'.format(var.name, errors['upper'], errors['lower']))\n",
"\n",
"# print(\"Function minimum:\", result.fmin)\n",
"# # print(\"Results:\", result.params)\n",
"# print(\"Hesse errors:\", result.hesse())"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# print(\"Time taken for fitting: {}\".format(display_time(int(time.time()-start))))\n",
"\n",
"# # probs = total_f.pdf(test_q)\n",
"\n",
"# calcs_test = zfit.run(probs)\n",
"# res_y = zfit.run(jpsi_res(test_q))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# plt.clf()\n",
"# # plt.plot(x_part, calcs, '.')\n",
"# plt.plot(test_q, calcs_test, label = 'pdf')\n",
"# # plt.plot(test_q, res_y, label = 'res')\n",
"# plt.legend()\n",
"# plt.ylim(0.0, 10e-6)\n",
"# # plt.yscale('log')\n",
"# # plt.xlim(3080, 3110)\n",
"# plt.savefig('test3.png')\n",
"# # print(jpsi_width)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# _tot = 4.37e-7+6.02e-5+4.97e-6\n",
"# _probs = []\n",
"# _probs.append(6.02e-5/_tot)\n",
"# _probs.append(4.97e-6/_tot)\n",
"# _probs.append(4.37e-7/_tot)\n",
"# print(_probs)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# dtype = 'float64'\n",
"# # mixture = tfd.Uniform(tf.constant(x_min, dtype=dtype), tf.constant(x_max, dtype=dtype))\n",
"# mixture = tfd.MixtureSameFamily(mixture_distribution=tfd.Categorical(probs=[tf.constant(0.007, dtype=dtype),\n",
"# tf.constant(0.917, dtype=dtype),\n",
"# tf.constant(0.076, dtype=dtype)]),\n",
"# components_distribution=tfd.Uniform(low=[tf.constant(x_min, dtype=dtype), \n",
"# tf.constant(3080, dtype=dtype),\n",
"# tf.constant(3670, dtype=dtype)], \n",
"# high=[tf.constant(x_max, dtype=dtype),\n",
"# tf.constant(3112, dtype=dtype), \n",
"# tf.constant(3702, dtype=dtype)]))\n",
"# # for i in range(10):\n",
"# # print(zfit.run(mixture.prob(mixture.sample((10, 1)))))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# print((zfit.run(jpsi_p)%(2*np.pi))/np.pi)\n",
"# print((zfit.run(psi2s_p)%(2*np.pi))/np.pi)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# def jpsi_res(q):\n",
"# return resonance(q, _mass = jpsi_mass, scale = jpsi_scale,\n",
"# phase = jpsi_phase, width = jpsi_width)\n",
"\n",
"# def psi2s_res(q):\n",
"# return resonance(q, _mass = psi2s_mass, scale = psi2s_scale,\n",
"# phase = psi2s_phase, width = psi2s_width)\n",
" \n",
"# def p3770_res(q):\n",
"# return resonance(q, _mass = p3770_mass, scale = p3770_scale,\n",
"# phase = p3770_phase, width = p3770_width)\n",
" \n",
"# def p4040_res(q):\n",
"# return resonance(q, _mass = p4040_mass, scale = p4040_scale,\n",
"# phase = p4040_phase, width = p4040_width)\n",
" \n",
"# def p4160_res(q):\n",
"# return resonance(q, _mass = p4160_mass, scale = p4160_scale,\n",
"# phase = p4160_phase, width = p4160_width)\n",
" \n",
"# def p4415_res(q):\n",
"# return resonance(q, _mass = p4415_mass, scale = p4415_scale,\n",
"# phase = p4415_phase, width = p4415_width)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# 0.15**2*4.2/1000\n",
"# result.hesse()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Constraints"
]
},
{
"cell_type": "code",
"execution_count": null,
"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": null,
"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": null,
"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": null,
"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": null,
"metadata": {},
"outputs": [],
"source": [
"# print('{0}/{1} fits converged'.format(len(Ctt_list), nr_of_toys))\n",
"# print('Mean Ctt value = {}'.format(np.mean(Ctt_list)))\n",
"# print('Mean Ctt error = {}'.format(np.mean(Ctt_error_list)))\n",
"# print('95 Sensitivy = {}'.format(((2*np.mean(Ctt_error_list))**2)*4.2/1000))"
]
},
{
"cell_type": "code",
"execution_count": null,
"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": null,
"metadata": {},
"outputs": [],
"source": [
"# Ctt.floating = False"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# zfit.run(nll.value())"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# result.fmin"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# BR_steps = np.linspace(0.0, 1e-3, 11)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# CLS Code"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"# zfit.run.numeric_checks = False \n",
"\n",
"load = True\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 = 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",
"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",
" \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))\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:\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",
" if result.converged:\n",
" Nll_list[-1].append(result.fmin)\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",
"# 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",
"# print(Nll_list)"
]
},
{
"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": {},
"outputs": [],
"source": [
"for s in range(len(l)):\n",
" print('BR: {:.4f}'.format(BR_steps[s]))\n",
" print(2*l[s]/len(CLs_values[0]))\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
}
