{
"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",
"import random\n",
"import numpy as np\n",
"from pdg_const1 import pdg\n",
"import matplotlib\n",
"import matplotlib.pyplot as plt\n",
"import pickle as pkl\n",
"import sys\n",
"import time\n",
"from helperfunctions import display_time, prepare_plot\n",
"import cmath as c\n",
"import scipy.integrate as integrate\n",
"from scipy.optimize import fminbound\n",
"from array import array as arr\n",
"import collections\n",
"from itertools import compress\n",
"import tensorflow as tf\n",
"import zfit\n",
"from zfit import ztf\n",
"# from IPython.display import clear_output\n",
"import os\n",
"import tensorflow_probability as tfp\n",
"tfd = tfp.distributions\n",
"\n",
"from matplotlib.pyplot import figure"
]
},
{
"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=-2.5, upper_limit=2.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": [],
"source": [
"# total_f_fit.normalization(obs_fit)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Test if graphs actually work and compute values"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"# def total_test_tf(xq):\n",
"\n",
"# def jpsi_res(q):\n",
"# return resonance(q, jpsi_m, jpsi_s, jpsi_p, jpsi_w)\n",
"\n",
"# def psi2s_res(q):\n",
"# return resonance(q, psi2s_m, psi2s_s, psi2s_p, psi2s_w)\n",
"\n",
"# def cusp(q):\n",
"# return bifur_gauss(q, cusp_m, sig_L, sig_R, cusp_s)\n",
"\n",
"# funcs = jpsi_res(xq) + psi2s_res(xq) + cusp(xq)\n",
"\n",
"# vec_f = vec(xq, funcs)\n",
"\n",
"# axiv_nr = axiv_nonres(xq)\n",
"\n",
"# tot = vec_f + axiv_nr\n",
" \n",
"# return tot\n",
"\n",
"# def jpsi_res(q):\n",
"# return resonance(q, jpsi_m, jpsi_s, jpsi_p, jpsi_w)\n",
"\n",
"# calcs = zfit.run(total_test_tf(x_part))\n",
"\n",
"\n",
"\n",
"test_q = np.linspace(x_min, x_max, int(2e6))\n",
"\n",
"probs = total_f_fit.pdf(test_q, norm_range=False)\n",
"\n",
"calcs_test = zfit.run(probs)\n",
"\n",
"Ctt.set_value(0.5)\n",
"\n",
"calcs_test1 = zfit.run(probs)\n",
"\n",
"Ctt.set_value(0.0)\n",
"\n",
"Dbar_s.set_value(0.3)\n",
"\n",
"DDstar_s.set_value(0.3)\n",
"\n",
"calcs_test2 = zfit.run(probs)\n",
"# res_y = zfit.run(jpsi_res(test_q))\n",
"# b0 = [b0_0, b0_1, b0_2]\n",
"# bplus = [bplus_0, bplus_1, bplus_2]\n",
"# bT = [bT_0, bT_1, bT_2]\n",
"# f0_y = zfit.run(tf.math.real(formfactor(test_q,\"0\", b0, bplus, bT)))\n",
"# fplus_y = zfit.run(tf.math.real(formfactor(test_q,\"+\", b0, bplus, bT)))\n",
"# fT_y = zfit.run(tf.math.real(formfactor(test_q,\"T\", b0, bplus, bT)))"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAaUAAAEMCAYAAABkwamIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOydeXzU1dX/32cmmZCNACEBwr6E3Q0i7lahKrYo1qVVu9jWyvO02k1al650sU/porY/tX14tIutFqnaipRiVcBdIKAoWyCEfUsgkITsM3N/f8x3kslklu+EzJJw3q9XX8zc77n3nvnWzGfOved7rhhjUBRFUZRUwJFsBxRFURTFj4qSoiiKkjKoKCmKoigpg4qSoiiKkjKoKCmKoigpg4qSoiiKkjKkJduBVGPgwIFm1KhRyXZDUVKerYdqye2TzrD+mR3aTzS0su94A+MH5ZKRltzfvR8eqKEwN4NBfftQXd/CgRONTBzcl3SnJNWv3sj69euPGmMKTnUcFaUgRo0aRWlpabLdUJTUZvt2Spbs5YrJg/if68/ocGnZBwe56+n3eO6blzK+5hCMH580H0f/cQd3XT6O+VdOYPHavdz3/IesuH8mQ47sS55fvRQR2dMd4+jynaIosTNnDsYYnCG+QZzii0I8XgNz5iTYsQDmzMEYEMsfR6r4pURERUlRlNjZvh2PMW1f9IE4HAFf/tu3J9qzNkxZGQB+D/1+GUNS/VIio6KkKErsLFiA1xtalNKsL3+vMbBgQYIdC8Ca29EWKfmak+6XEhEVJUVRuoTXEDFScnuTW1fTX9bT72KH5TslZVFRUhQldhYswBtlT8nrTW5E4vnhD4HOy3deg0ZKKYyKkqIosTN+PJ4wy3fOwD2lJGa4pU2aCLSLUYflO828S1lUlBRFiZ1lyzCm/Qs/kLZlMmNg2bJEe9ZGyz9e6PC+LYJLsl9KZFSUFEWJnbo6K/uu86U0p3/5zmeXNE765vaLpEiK+KVEREVJUZTYmTcPj9e0RR+B+EXA7fXCvHmJ9qyN9C//NxCY6OD712tMUv1SImNLlERktoiUiUi5iNwX4nqGiDxjXV8jIqMCrt1vtZeJyFXRxhSR0dYYO6wxXZHmEJF8EVklIidF5JEw/i8VkU32bomiKNEwVtWTUMt3zsCU8PXrE+pXII3vrAXaEx1SxS8lMlFFSUScwKPA1cBk4BYRmRxkdjtw3BgzDngIWGj1nQzcDEwBZgOPiYgzypgLgYeMMcXAcWvssHMATcD3gW+F8f964GS0z6koin3M3fOB0Cnh7RUdgPnzE+lWB9Lv/TYQ+JxSQPZdEv1SImMnUpoBlBtjKowxLcBiYG6QzVzgz9brZ4FZ4lvAnQssNsY0G2N2AeXWeCHHtPrMtMbAGvO6SHMYY+qNMW/iE6cOiEgOcDfwUxufU1EUm3iLioD26COQDtl3ll0y8A4eArQv3/n/TbZfSmTsiNJQYF/A+/1WW0gbY4wbqAHyI/QN154PnLDGCJ4r3ByR+Anwa6AhkpGIzBORUhEpraqqijKkoijub3wTaP+iD6SDKCUxImn+ut9H6eCXMcn1S4mMHVEKVeM9+JHocDbd1W7Xj3aHRM4Gxhlj/hHOpm0QYxYZY0qMMSUFBadceV1Rej2ukcMBQiY6+B+o9ZjkRiS5Y0YCAQ/PSmpEcEpk7IjSfmB4wPthwMFwNiKSBuQB1RH6hms/CvSzxgieK9wc4bgAmC4iu4E3gfEisjriJ1UUxRb1b74DhF6+cwRWdEjiMTC1b7wNdC4z5DUk1S8lMnZEaR1QbGXFufAlLiwNslkK3Ga9vhFYaYwxVvvNVubcaKAYWBtuTKvPKmsMrDFfiDJHSIwxvzPGFBljRgEXA9uNMZfZ+LyKokSjzFdlO+qeUhKrcTt27PD9G1SQ1Zjk+qVEJqooWfs3dwEvAVuBJcaYzSLyYxG51jJ7AsgXkXJ8iQX3WX03A0uALcAK4E5jjCfcmNZY9wJ3W2PlW2OHnQPAioYeBD4vIvtDZAcqitKN9PnZT4AoopTkatyZlo9tkVKK+KVExtbJs8aY5cDyoLYfBLxuAm4K0/cB4AE7Y1rtFfiy84LbI80xKor/u4GpkWwURbFP9b/+Az97NXqktHp1gj1r58Tyl+GBVzod8uc1JNUvJTJa0UFRlJjJ+dpXgDCJDoEJBUmsnJDzVZ+P7YkOvn+1okNqo6KkKErMNJ81DQiT6BBYOaGkJKF+BdJ6js/HTokO3uT6pURGRUlRlJip/ewXgNCilOZIjUip4fO+YjCOoOeUvAaNlFIYFSVFUWJmxKjBQORIyeM1kJubUL8CKRzme+bQ72GHig5J9EuJjIqSoigxs3NjGRAm+y5wT+lg8CONiePwlp1A50jJmOT6pURGRUlRlJjp8+YbQPtSXSAdUsKTmOWW/sbrvhehHp7V7LuURUVJUZSYyfvrH4HQVcI7JBQsWpRQvwLJefIPQOfsO49Jrl9KZFSUFEWJmV1/WAy0nzIbSHuiA/Dii4l0qwNVT/kOGwg+usIYk1S/lMioKCmKEjPD7/oSECZSahMlL9x6a0L9CmTAvM8DnVPCPV6TVL+UyKgoKYoSM8dnXglAmiP0V4jTIb5lsjlzEulWBxquvBoIc8hfEv1SIqOipChKzFTOuQGAMJqEU8S3fJfEiKTuel9Vsvbad75/vUYjpVRGRUlRlJi5YNxAIHyk5HBYX/6hTgFMEGMLOz6L1CEBI4l+KZFRUVIUJWZeK6sE2g/0CybN4cDtMRD+dJm4U36kFghT0SGJfimRUVFSFCVmBi71ZbY5w0VKYkVKTz+dSLc6kP3cEqA9KPL/m2y/lMioKCmKEjMDVr0MhK4SDlaig9fAsmWJdKsD2S+vAEIlOiTXLyUyKkqKosTMxoWPAaHLDPnaHbi9yY1IDjzqOx/U76EzcE9JI6WURUVJUZSYOefOzwLhRSndKbg9XrjmmkS61YFhn/sUEOI5JUNS/VIio6KkKErM7L7h00B4UUpzii9SSuIREdWfvg2g7eRZsb7tjB7yl9KoKCmKEjOV084HIkRKDgetHi9cdlkCvepI/QUXAyGW74xJql9KZFSUFEWJmdlX+U5uDVUlHKxIyWOgqCiRbnVgwvRJQOdEB4+XpPqlRMaWKInIbBEpE5FyEbkvxPUMEXnGur5GREYFXLvfai8TkauijSkio60xdlhjuiLNISL5IrJKRE6KyCMB42SJyL9EZJuIbBaRn8d+exRFCcXzq7cA7XXugklzOHB7vVBXl0i3OrBp6z4gTEWHJPqlRCaqKImIE3gUuBqYDNwiIpODzG4HjhtjxgEPAQutvpOBm4EpwGzgMRFxRhlzIfCQMaYYOG6NHXYOoAn4PvCtEO7/yhgzETgHuEhEro72eRVFic6o533Za+EipXSn0OpJ7hERA576E9AuSv7qE54kH6mhRMZOpDQDKDfGVBhjWoDFwNwgm7nAn63XzwKzxLe7OBdYbIxpNsbsAsqt8UKOafWZaY2BNeZ1keYwxtQbY97EJ05tGGMajDGrrNctwAZgmI3PqyhKFAZs2QiErhIOkOa0IqXS0kS61YHMD94D2hMd/Prp9pqk+qVExo4oDQX2Bbzfb7WFtDHGuIEaID9C33Dt+cAJa4zgucLNERUR6QdcA7wa5vo8ESkVkdKqqio7QyrKac1b9/wMiLCn5BBa3cmNSPb8z8NAe6KDiJDmEN+RGhoppSx2RCnUf3XBhaPC2XRXu10/OiEiacDfgN8aYypC2RhjFhljSowxJQUFBdGGVJTTnqvuugUIv6fkSnPQ6k1u9t34T/meRQqM5pwOKwFDs+9SFjuitB8YHvB+GHAwnI0lAnlAdYS+4dqPAv2sMYLnCjdHNBYBO4wxD9uwVRTFBqWf/yoQOVJyewwsWJBArzpy4Bv3AB0Lgqc5rOenkuiXEhk7orQOKLay4lz4EheWBtksBW6zXt8IrDTGGKv9ZitzbjRQDKwNN6bVZ5U1BtaYL0SZIywi8lN84vUNG59TURSbHBs6Coj08Kz1nNL48Qn0qiMNI8cAHevzpTkdvkSHJPqlRCaqKFn7N3cBLwFbgSXGmM0i8mMRudYyewLIF5Fy4G7gPqvvZmAJsAVYAdxpjPGEG9Ma617gbmusfGvssHMAiMhu4EHg8yKyX0Qmi8gw4Lv4svs2iMj7IvKl2G+RoijB3PTf1wNRygx5DZSUJNKtDky97qNAe6ID+CMlb1L9UiKTFt0EjDHLgeVBbT8IeN0E3BSm7wPAA3bGtNor8GXnBbdHmmNUGNf1JC9FiQP/t/gNeGVH2CrhvvOUvHAweKU/cax7axM8voZA3WyrXp5Ev5TIaEUHRVFi5pwlvgWMsA/P+p9T+vWvE+lWB4Y8/ijQMZpr2+tKol9KZFSUFEWJmaxjlWGX7sBX+87tTW6klH7kMNBx+c7p1Egp1VFRUhQlZl65/Z6wmXcQUPsuiRFJ+b0/AoIjJeucJ42UUhYVJUVRYua2r91IujP810e6P/tu+vQEetWR6TdcARB6TymJfimRUVFSFCVmXvjyD0h3RoiU/M8DJbFywuYf/gLo+PBsW/adVnRIWVSUFEWJmYaMTNIiREppTodv+S43N4FedaQ1KxsIU9EhiX4pkVFRUsJS3+xm/Z7jyXZDSUE+97Ovkh4p0cEpvjJDc+Yk0KuOTL/L96x9p+w7r0mqX0pkVJSUsHzjmfe54XdvU13fkmxXlBTjgQf/GTlScjgwBjzbyhLoVUdWvvAG0HFPqa2iw/btSfJKiYaKkhKWTQdqAGhq9STZEyXVuOLvvyct0p6Sdc388IeJcqkTxb/7FdDxWSqnf09Ja9+lLCpKiqLEjNdrSHdEyr7zCYE3ah3/OGLNHZzo4EmqU0o0VJQURYmZf8y9g/S0SNl3vq+W+vu/lyiXOrF53t1A55RwrRKe2qgoKYoSMw/cd0Ob8ITCHynlnjUlUS514srrLwXCRErjx/Pe3uOMuu9fVFSdTJaLSghUlBRFiZlf3PmriM8p+R+sPbb4uUS51InXf/0HIHhPyUpVX7aMv7yzB4DS3ZphmkqoKCmKEjOuhpMRIyV/Zp63pi5RLnXCWV/v+zdUpFRXR0a6z8cWjzcp/imhUVFSFCVm5j21MGL2nT+KGnD3XYlyqRPn/o/vyLUOe0pOK/tu3rw2UfVGPitUSTAqSkpU9E9WCebr3348Yu07/xf+7pdeT5RLnVj+xxeBjst3bQ/Prl/f9lBtq0f/C08lVJQURYmZLzz/SMQq4a4031fLgB/cnyiXOnHOb34KhCkzNH9+myjpc3iphYqSEhU9vlcJpjInn/S08F8fGda1psJBiXKpEw0DfXOH3FMqKmprU1FKLVSUlKjo4oYSzN8uuSli7Tt/pLTv819OlEud+OBTtwMgAd9yaU7rPKX5831HawANLSpKqYSKkqIoMfPcT2+KWPvOHylNv3BqolzqxC2fuAAIFSl5oaiIFrdPlBo1UkopbImSiMwWkTIRKReR+0JczxCRZ6zra0RkVMC1+632MhG5KtqYIjLaGmOHNaYr0hwiki8iq0TkpIg8EuTXdBH50OrzWwk8F1mxjd40JZjPffmRiM8pZaQ5AXh78YpEudSJ5/73H0CIPSWvgdLStlRwXb5LLaKKkog4gUeBq4HJwC0iMjnI7HbguDFmHPAQsNDqOxm4GZgCzAYeExFnlDEXAg8ZY4qB49bYYecAmoDvA98K4f7vgHlAsfW/2dE+r9IZXb5TghlauS/ic0r+5bv0ih2JcqkTuXt3ARD4U7RtT2n79rasOxWl1MJOpDQDKDfGVBhjWoDFwNwgm7nAn63XzwKzrKhkLrDYGNNsjNkFlFvjhRzT6jPTGgNrzOsizWGMqTfGvIlPnNoQkSFAX2PMO8YYAzwZMJaiKKfAf636S8TnlPzLd+N/92CiXOrEjD//Fgg6T8l/TPuCBbS4fWLUqHtKKYUdURoK7At4v99qC2ljjHEDNUB+hL7h2vOBE9YYwXOFmyOS3/uj+A2AiMwTkVIRKa2qqoow5OmJLt8pwXz2swsjPqfkF6X/PLo4US514tlf/RXouHyX7nTQ6jGYVasCIiWt6JBK2BGlUN9JwSs64Wy6q92uH3Z86txozCJjTIkxpqSgoCDCkKcnunynBLNg2W9s7SmV/KzTFnTCmPmgr0J5YJKgXyy9d8zTRIcUxY4o7QeGB7wfBhwMZyMiaUAeUB2hb7j2o0A/a4zgucLNEcnvYVH8VhQlRowxbBxUHHFPyV9X7nBx8rLvDo2fgggE5je5rOiu9ZxpmuiQotgRpXVAsZUV58KXuLA0yGYpcJv1+kZgpbWPsxS42cqcG40v2WBtuDGtPqusMbDGfCHKHCExxhwC6kTkfGuv6nMBYymK0kVaPYa/nT07YqTk//J//+qbEuVWJzZe/ckO6eDQXpOv4fO3tz2npKKUWkQVJWv/5i7gJWArsMQYs1lEfiwi11pmTwD5IlIO3A3cZ/XdDCwBtgArgDuNMZ5wY1pj3QvcbY2Vb40ddg4AEdkNPAh8XkT2B2TyfRl4HF+CxU7g37HcHEVROuP2etn0UOTnlBwOId0pfPnacxLoWUe+9onpHfaTAFzWsmK/wfltoqTLd6lFWnQTMMYsB5YHtf0g4HUTEPInkTHmAeABO2Na7RX4svOC2yPNMSpMeymQvPUDRemFtLi9XPSVP/OtCGWGwLev9PCfVnNvgvwK5pEnX8OxsWPikj9V/cCH5bj/WQZookOqoRUdFEWJiWa3l/P3fdj2BR8OV5qDoo3vJsirzgz7YE2ISMnns/ON13wP0aKRUqqhoqQoSky0uL3c8v6Ktgy7cGSkOZi+4tmINvHknBV/7yxK/mPan/wjbmv5rsXt9T1Qq6QEKkpKVCLkkyinIc1uD1+68YdRI6WMNAf/++2HE+RVZ57+3qME14z1+7zzD3/rcI5Sszv+0dK63dVc+otVlFcm7zTenoCKkhIV1SQlkKZWL79Z+su2Z37C4Upz8OnfJO88pet/9e0O1RwAXE5fdDf8ri/5TqC1SERVhxfeP8De6gb+s+VI3OfqyagoKYoSEy0eL6+OO9dGpOTkvTMvSpBXndky/dJOGYJ+n49+5Ao8XtOWup6IfaU9xxoAKK88Gfe5ejIqSoqixERzq5elky+LGillpDl4/dwrEuRVZzZcdHWn03H9zykd/NgnaPUYcvv4EpATkYF3oqEVgEMnmqJYnt6oKClh0WU7JRQtHi+7F86xtXz31B0XJMirzvziprM7L99ZPl8+aRBuj5ecNlGKf6RU02iJUk1j3OfqyagoKVFRcVICaW71MOreZbay7+b89vUEedWZr/9tQ6dIyS+kL75/gFavIScjcaJU2+QXpSZNHoqAipISFn82rdGSrEoALR4v125ZHXVPKdPl5IJ3X0qQV505+43lnSIlf2XzQcuew+NtX76L956S12uobWzF5XTQ7PbqEewRUFFSwqI/5pRQNLd6mVW+LuryXWZ6GuduejtBXnVmyoY3OhWN9Qtp4euv4PEacjLSgfhn39W3uPEaGD0wG4Dq+pa4zteTUVFSoqLipATS4vHy9Wu/HTVSynI5ueeG5B1dsei/f9LpIEJ/tt2r3/c9P9WW6OCOb6KDfz/JL0pHTzbHdb6ejIqSoigx0dzq4fFnfxR1Tykrw8mDf/l+grzqzF0P3d1pT8kvpDO/fTtA+55SnCOl2kbfuaWjCzRSioaKkhIVDZSUQFo8Xv529uyoy3dZ6Wk8ddZVSSvhs/KSuWH3lNZf5avtnKg9pbZIKd8nSsdUlMKioqRERTOFlECaW728O/wMW8t37w4/g4YWd4I868iHxed02lPyC+m2idMAEpYS7s+880dKx06qKIVDRUlRlJho8XhZ89htnZbGgsl0OVnz2G1JyzR77J5rOkVKIkKfdAff+sJMAHIzEhspDe7bhz7pDqrrdU8pHCpKSlQ0TlICaXZ7KbnnuQ7HjIciy+Vk6jf/njRRuvXBlzslOgBkudK483erAd8ek8vpiHtFh1pLlPpmppOfnaHLdxFQUVLCos8nKaFocXu59f3ozx9luZzc8v6KpC3fffSNF0JGc5npTma87DtSI83hoE+6I/7Ld42tiPgiswHZLl2+i4CKkhIV3VJSAmlq9XDG4R1R7bJcaZxxeEdCKnCHYsyebTgdnb/islxOhldsASDNKWS6nHEXztomN7kZaTgcQn6Oi2O6fBcWFSUlLIL/V6aqktJOY6uHhz/5rah2WS4n35n91aQt3z1807dCRkpZLic/vfYbgC9Syu2TTl1TfEWpprGVvCzfg7oDczI4WqeRUjhUlJSw6PKdEoqGFg//73/vjmqX6XKy+On7krZ8t/CRr+EMsaeU6XLym//9JgBOh9C3T1pbdly8qG1spW+fdlE6Vt+sWa1hsCVKIjJbRMpEpFxEOj2iLSIZIvKMdX2NiIwKuHa/1V4mIldFG1NERltj7LDGdJ3CHN8Ukc0isklE/iYifWK7PQro8p3SkcYWD0vmfDGqXZYrjYcvvjVpkdIfrwidIZiZ7uTXF94CQEa6g7zM9LaHW+NFTWMreZl+UXLR6jFxn7OnElWURMQJPApcDUwGbhGRyUFmtwPHjTHjgIeAhVbfycDNwBRgNvCYiDijjLkQeMgYUwwct8buyhxDga8BJcaYqYDTslMU5RRobPVwfNioqHZZLicV/YcmTZR2Dxja6Tkl8IllRf+hgO+5pb6Z6W0p2/Gitqk9UirIzQCgSksNhcROpDQDKDfGVBhjWoDFwNwgm7nAn63XzwKzxJcvOhdYbIxpNsbsAsqt8UKOafWZaY2BNeZ1XZwDIA3IFJE0IAs4aOPzKkFooKQE0tDi4YEH7ERKTl588ptJW757/JGvhI6ULL/AdzpuXmZ63Jfvahpb6ZvpeyYqP9snSlr/LjR2RGkosC/g/X6rLaSNMcYN1AD5EfqGa88HTlhjBM8V0xzGmAPAr4C9wCGgxhjzn1AfUETmiUipiJRWVVWFvRGnK7p8pwTS2OLmB4+F/FPqQLYrjfPvepKTcU4iCMfV9ywOuaeU5XJy3p1PAlak1Ced2sZWvHEsh1Tb6G5fvst1ASpK4bAjSqGekAv+fy+cTXe1xzyHiPTHF0WNBoqAbBH5TAhbjDGLjDElxpiSgoKCUCaKolg0tHi4csVfo9o5HMJX1r9AbZJE6dY3loSNlL609nkA+lh7Sl7jO14iHrS4vTS2ejokOgAcrVNRCoUdUdoPDA94P4zOy2BtNtZSWR5QHaFvuPajQD9rjOC5Yp3jo8AuY0yVMaYVeB640MbnVYLQLDwlkMZWDwNqjtmyHd54vK2aQaLJrz3WqcwQ+ArFDjpZDfiW7/zLavHaV/KP608J75/lwiFwVB+gDYkdUVoHFFtZcS58yQJLg2yWArdZr28EVhpfvuNS4GYrc240UAysDTem1WeVNQbWmC90cY69wPkikmXtPc0Cttq7LYqihKOxxcOb/23vnKQnb/xa3PdrwvHArC+FLBqbl5nGAzO/BPiW7/zLanEXJWsep0MYkJ2hy3dhiCpK1v7NXcBL+L7UlxhjNovIj0XkWsvsCSBfRMqBu4H7rL6bgSXAFmAFcKcxxhNuTGuse4G7rbHyrbG7MscafAkRG4APrc+6qAv36LRH95QUPy1uL26v4fNfuzG6MfDIg3ckJfXZGMNzj3+VDGfnr7j+2S5e/NPXAV/tuwFW4kG8zjiqCah752dgjksjpTCkRTcBY8xyYHlQ2w8CXjcBN4Xp+wDwgJ0xrfYK2rPnAtu7MscPgR+G6qPYR0VJ8eOvpv3WPQ/wCRv2T3/xu0mJlFo9hvtnf5XZISKlAdku7p/9VQCyM9LaUrQra+MTudQGRUrgSwvXSCk0WtFBURTb+OvYpeXl2bJP65eXlD2lFo+XeldmyOW7/lku6l2ZgO/Qv0K/KMUp8SB4+Q58olRZ2xSX+Xo6KkpKVDTRQfHjf+bo8m9Ff04J4I5ffC0p2Xctbi9PPPujtpNmAxmQ7eKJZ3/U9j47I41sl5OqOItSvwBRGtovk8O1Tbg98T0yoyeioqSERZftlGD81Rne/Nfbtuz/8uf/cLLZnfAv3xa3l5nzFoWMlApyM5g5r+P2cmHfPlTWxSdyCbWnNLRfJl4DhzVa6oSKkhIVFSfFj7+a9pRFD9qy/8jTjwJwsjmx0VKL28s33nwKV4hIKd3p4IWql1n02eltbQW5GXFbvjvR0Eq2y9khaivq51s+PHhCRSkYFSUlLFEOFlVOQ+qspIVQEUgoMtKcQPzSrcPR4vFFdOH8PGt4P66cMrjtfWFuBkfiFLUEFmP1M7S/T5QOnGiIy5w9GRUlJSwaISnB+COlxu9835b94bvvBeKXbh2OZreXhy/+NBnhxHPBgg5vRwzI4sDxxrgsM55oaCEvy9WhrShPI6VwqCgpUVFxUvz4l+GGn3+2LfuLPnYRQMKP/25xe1m5aF74iG78+A5vR+Zn4fYaDtV0v0hUnWxuy/Dzk+lykp/tYv/xxm6fr6ejoqQoim38y3eepS/asj+xxFdjLtHP5LR6DLff+ENcTmdog2XLOrwdmZ8NwO5j9d3uS1Vdc9uzUIEM7Z/J/uO6fBeMipISFU0JV/zUNbnJSHPgarT35d3P44s8jiV4+a7F7SW7pTF8pFRX1+HtyPwsAPYc616R8HoNR0+GFqXRA7OpqOp+EezpqCgpUdHlO8VPbZOb3D5pMG+eLfuMr3yZ3Iy0uD0DFI4Wj4f/WfH/wotSkP+DcvvQJ93R7SJR09hKq8d0Wr4DGFeQw4ETjW0PJCs+VJQURbFNXVMruX3SYf16ex3Wr2dgEkrqtLi9XPP534RMCff7FYjDIUwY3Jcth2q6NN+f3trFz/+9DU/QmUz+NHP/cRWBjC3MAWBn1ckuzdlbUVFSoqKBkuKnzh8pzZ9vr8P8+eRnuxKe6NDU6uW7Kx8nIz3MV1wI/6cW9WXzwVpMjEsD5ZUnWfDiFn7/2k6eWbevw7V91b7lwGFWCnggYwtUlEKhoqSERcVICcYXKaVBUZG9DkVFDMxJfKTU0OLhSM4AslxhEh1C+D91aB51TW72VRjGMHAAACAASURBVMeWEffGDt9p1QW5GTz+ZkWHE2z3WKI0ykqkCGRkfhYOgZ26r9QBFSUlKrH+clR6L7VNbt8JqjFESvGslhCOhhY3j8+4nixXmIMQQkZKviKz7+8/EdNc5ZUnyctM5zsfm0hFVT1v7Tzadm3PsXpy+6TRLyu9U78+6U5G5Wez7VBtTPP1dlSUlLD4CzqoJCl+qutbGJDtiilSKuqXSU1ja1s6eSJobPGw5tHPxRQpTRqSS25GGu/stHeqrp8dlScZV5jDx84YQn62iyff2dN2bfuROsYU5CBhyqNMHZrHpgNd28fqragoKWFRMVIC8XgNJxpayM92QWmpvU6lpW37KQdOJO5B0YZWD9d/4eGQVcL9fgWT5nRw3pgBvB0Q6dhhZ+VJigtzyEhzcsuMEby69Qj7qhvweA2bDtRy1rDwx3ycOSyPgzVNCc9OTGVUlJSo6OqdAr70Zq/xndzK9u32Om3f3l7nLYHVCxpbPEyoPRTeIIz/F44dyJ5jDW0JCtGorm/hWH0L46xMulvPGwHAU2v28sH+E5xsdjNtRP+w/c8Y6hOsDw/EtmTYm1FRUhTFFv76dQOyXZ1qx4VlwQKG9UtCpNTi5q7XnwpvEMb/j04aBMDyDyMIWgDllb7MOX96d1G/TK6cPJin3t3DwhXbcDkdXD6xMGz/qUPzEIGN+3QJz4+KkmIDDZWUIFFavdpep9WrGZiTgcvpSGik1NDi4Vtf+U1Ev0IxIj+Ls4blseyD2ESp2BIlgO9+fBKuNAfvVlQz79IxnSqEB5KdkcaEQbms210dda5dR+upOA3Sx1WUlKjo8p0CQaJks6ID8+bhcAhD+2eyL4F13hpbPNzzfIQznyL4f81ZRXx4oIayw3VhbfzsqKwjM93ZVvUbYPiALF6d/xGW3nUR868cH6G3jwvHDmT9nuM0tYav7NDY4uHjv32DuY+81etPq7UlSiIyW0TKRKRcRO4LcT1DRJ6xrq8RkVEB1+632stE5KpoY4rIaGuMHdaYrlOYo5+IPCsi20Rkq4hcENvtURTFTwdRKimx18myGzMwuy2qSAQNLR72jJ4U3iCC/zdMG0ZGmoM/vb0r6jzllScZW5iNw9Exu65floszh/ULm3UXyIVj82l2e3lvb/h9pff3naChxUNds5v39/Xu/aeooiQiTuBR4GpgMnCLiEwOMrsdOG6MGQc8BCy0+k4GbgamALOBx0TEGWXMhcBDxphi4Lg1dsxzWH1+A6wwxkwEzgK22r0xSjsaKCkA1fW+DLH+WbFFSgDFg3LZdbSe1gT9yq9rbmXNFTdG9SsU/bNdXD9tGM9vOBD1od/yypOMK8iJaBONGWMG4BB4J0LWX2Dixcb9vXv/yU6kNAMoN8ZUGGNagMXA3CCbucCfrdfPArPE9xNhLrDYGNNsjNkFlFvjhRzT6jPTGgNrzOu6MoeI9AUuBZ4AMMa0GGN6908MRYkjh2ub6J+VTp90J+Tm2utk2RUX5tDqMd1ehTscJxpaeezLl0X1Kxx3XDIat9fw21d3hLWpbWrlUE0T4wfbvBdh6NsnnbOG92P19qqwNnuq63E6hIE5Lj6M8eHenoYdURoKBBZ02m+1hbQxxriBGiA/Qt9w7fnACWuM4LlinWMMUAX8UUTeE5HHRaRzrQ9AROaJSKmIlFZVhf8P43RF95QUgMM1TQzq28f35uBBe50su/GDfF/cO45E36fpDmoaW3nwj6vCG0Txf0xBDrfMGM7Ta/aGrU3n/yzjC09NlACumDyID/bXcKgmdDLI3upGhvbL5Mxh/djSyytA2BGlUIuiwV9T4Wy6q70rc6QB04DfGWPOAeqBTvthAMaYRcaYEmNMSUFBQSiT0xotM6SAL1IakmeJUgzZdwDjCnMQgW02kgdOFY/XUNfkZsr2DVH9isTXZ40n0+Xkvuc+6FDPzk/ZYZ9YTTjFSAngysmDAXhly5GQ1/ceq2dkfhaThuSys6o+YlJET8eOKO0Hhge8HwYE/8xosxGRNCAPqI7QN1z7UaCfNUbwXF2ZY78xZo3V/iw+kVIUpQscrmlisF+UFi2y18myy3Q5mTAol/cSsElf2+grZ3TOS3+P6lckCnIz+OE1U1i3+zh/eKtz0sOHB2rIzUhjaL/OFcBjZVxhDmMKsvlPOFGqbmD4gCwmDemLx2vYcaT3pobbEaV1QLGVFefCl1SwNMhmKXCb9fpGYKXx/bxeCtxsZc6NBoqBteHGtPqsssbAGvOFrsxhjDkM7BORCVafWcAWG59XCULjJKXZ7eHoyRYG97W+gF+0dxx6oN05I/rz3t7jIaOO7qTGEqXS//fn8EY2/b9h2lCumDyIn/97G+9WdKyJt3bXMUpG9e+UeddVZk8ZzNs7j3UqOVTb1MrxhlZGDMhi8pC+AGztxUt4UUXJ2r+5C3gJX/baEmPMZhH5sYhca5k9AeSLSDlwN9YymTFmM7AEnxisAO40xnjCjWmNdS9wtzVWvjV2zHNYfb4KPCUiHwBnAz+L9QadzuiqneKnstb3Rdm2fHfrrfY6BthNG9GPuiZ33M8POt7gS12/8Adfs+VXJESEX3/yLEbkZ/Hlv66nvNK3/LjraD07q+q5cOzAU/bXzyfOGYrHa3jh/QMd2vdaySEjBmQxMj+bzHRnr95XClPXvSPGmOXA8qC2HwS8bgJuCtP3AeABO2Na7RX4svOC27syx/uAzQcqlHCoOCn+B1+L/EtVc+bY6xhgVzJqAADvVhyjeNCp78OEw39MRstVH7PlVzT69knn8c+V8KlF7/LJ/32XX9xwJss3HcLpEK45y2a1dBsUD8rlrGF5PLfhAF+6ZExb+14rHXxkfhZOhzBhcO7pHSkpitEFvNOe3Uetw+oGZvkauhApjcrPYmR+Fiu3VXa3ex04UtsEQJ/Pf8aWX3YYU5DD3//rAvplpvOlJ0t5fsMB7rhkTPseWzdxw/RhbD1Uy5aD7aLjT6MfaR0UOGlIX7Yeiv2E3J6CipKiKFHZc6weV5qjvZyOjUoFwXYiwsyJhby18xgNLe4InU6NI7VNOB1CYd8ICQh2/Q9g1MBs/v2NS3jk1nN44rYS7p09IXqnGLn2rCIy0hz8dU37mUx7q+vJz3aRk+Fb2Jo8JJfaJjcHa5q6ff5UQEVJiU7v/EGmxMCuo/WMHJDVvqlv91d6kN0VkwbR4vby6tb4RUuHa5opzM2I7GMXo4yMNCdzzixi1qRBtkoIxUq/LBfXTxvKc+v3t5V1qqjypYP7meRPdjjYO5fwVJQURYnK7mP1bctHADz9tL2OQXbnj8lnaL9MlpTuC9Ph1Nl3vMGXph3JR7v+J4EvXjSaZreXp97dg9dr2HKwlilF7QcFTuzlGXgqSkpUNFA6vWn1eNl9tIGxBQGitGyZvc5Bdg6H8MmS4byx42jcjmGoqDrJ2IKcyD7a9T8JFA/KZdbEQha9UcGb5Uepa3ZzZsDptTkZaYwYkMXWwypKymmHypECO6tO0uLxMrmob3tjFyMl8J3OmpnujFhXrqucaGjh6MkWxhZm99hICeD+j02iqdXD5/6wFqdDOh0UOHlIX7YeSkzJpkSjoqREpZcm+Sg28WeC+R/cBOCaa+x1DmFXkJvB5y4cyQsbD3b7MQz+L+riwtzIPtr1P0mMK8zhlzeexYgBWdw3eyIDczI6XJ80pC+7j9XHNWEkWagoKRHwbeRqSvjpzeaDtfRJdzAm8IiGGI+uCObOy8cxuG8fvvX3jTS2dF8dt/V7fCe4ThvRP7KPdv1PItedM5TX77mcOy4d0+napCG5GEOvjJZUlJQIqBgpvhpvEwb3xRlYTueyy+x1DmPXt086C284k51VJ/na4ve67TTVt8qPMWFQLnlZ6ZF9tOt/inLmsH4AvLf3eJI96X5UlJSo6PLd6UtTq4f3953g3JH9O14oslnJIILdpeML+OGcyby85Qj/9Zf1nGw+taWoyrom1uw6xpVTBkX30a7/KcrgvD6MzM9iza7qZLvS7agoKYoSlo37TtDi9nLemPyOF+psLhtFsfv8RaP5yXVTWVVWyVUPvc7LW450uVLBE2/uwuBb9oo6t13/U5gZowawbnd13AvcJhoVJSUqves/eSUW1uyqRgTOHRUUKcV4dEUkPnv+SJb81wX0SXdwx5OlXP2bN3jizV0cPBH6wLtQvLPzGE+8sYvrzh7qSwePNrdd/1OY88bkc6KhlbIEHZyYKGwVZFVOb3prjS0lOiu3VTK1KI9+Wa6OF0pL7SUL2LQrGTWAf3/9UpZuPMgf3tzFT5Zt4SfLtjCsfybnjOjPuIIcRuRnMiQvk5yMNLJcThpbPRyuaWJVWSXPrNvHiPwsFlwzxd7cdv1PYc4b7Stwu3ZXdVuVh96AipKiKCE5UtvE+/tO8K0rx3e+2I2Rkh9XmoMbpw/jxunDqKg6ycptlby39wQb9hxn2QcHw+5tupwOrj9nGPd/bKIvwcHO3L0gUhrWP5OivD68s/MYt104KtnudBsqSkpY/EvVGiednvxn82EArpwyuPPFyy6zdyS6XbsgxhTkdEhBb2r1cOBEI4dONFHf4qaxxUNGmoOBuRlMLcoj0+WMbe4u+pVKiAiXFBewfNMh3B4vac7esRujoqSERZftTm+WlO5nwqBcigtzOl9csMDeIHbtotAn3cnYgpz2/aJTnbub/Eo2H5lQwDOl+3h/34m286p6Or1DWpW4YDq9UE4XPtxfw4cHarj1vBGhq2GPD7GkFwq7dvEg0tzJ9KsbuWjsQBwCr22vSrYr3YaKkhIW07Z8p6p0uvHEmxVkpjvb06uDKbF5oLNdu3gQae5k+tWN5GWlc86I/ipKyumBLt+dnpRXnmTpxoN87oKR5GWmhzY6eNDeYHbt4kGkuZPpVzfzkfEFfHighmMnm5PtSregoqSEpS1SUm06rfj5v7eRkeYMWXOtjV//2t5gdu3iQaS5k+lXN/OR8QUYA2+WH022K92CLVESkdkiUiYi5SJyX4jrGSLyjHV9jYiMCrh2v9VeJiJXRRtTREZbY+ywxnR1dQ7rmlNE3hOR1D1AJUVRLTr9eGnzYV7ZeoRvfLS4U2XqDmiklDKcMTSPAdkuVpf1jiW8qKIkIk7gUeBqYDJwi4hMDjK7HThujBkHPAQstPpOBm4GpgCzgccskYg05kLgIWNMMXDcGjvmOQJ8+zqw1d7tUALxL99ppHR6cLimie/+40MmDs7lixePjmyskVLK4HAIlxYP5LXtVb2i5JCdSGkGUG6MqTDGtACLgblBNnOBP1uvnwVmiS9lZy6w2BjTbIzZBZRb44Uc0+oz0xoDa8zrujgHIjIM+DjwuL3boQRigv5Vei9NrR7ufHoDDS0eHrn1HNKjPfMyfbq9ge3axYNIcyfTrzhw+cRCqutb+OBATbJdOWXsiNJQYF/A+/1WW0gbY4wbqAHyI/QN154PnLDGCJ4r1jkAHgbuASLWxReReSJSKiKlVVW9IwTuDjRCOj1o9Xi56+kNrN9znF/ceCbjCnOjd4pDRYdup5dXdAjkkuICRGDVtspku3LK2BGlEA8pdPrxHM6mu9pjnkNE5gCVxpj1Ia53NDZmkTGmxBhTUlBQEM38tMHbtnyn6tRbaWzxcOdTG3hlayU/mTuFOWfaPNIh14ZwxWIXDyLNnUy/4sCAbBdnD+/H6l6QGm5HlPYDwwPeDwOCdwnbbEQkDcgDqiP0Ddd+FOhnjRE8V6xzXARcKyK78S0PzhSRv9r4vIqFSlHvprKuiZv/711e3nqEH14zmc9eMMp+5zlzutcuHkSaO5l+xYnLJxTywf4TPT413I4orQOKraw4F76kgqVBNkuB26zXNwIrje/n9VLgZitzbjRQDKwNN6bVZ5U1BtaYL3RlDmPM/caYYcaYUdb4K40xn7F5XxRoUyUVp97Hym1HuPrhNyg7XMvvPzOdL1wUJbEhmO3bu9cuHkSaO5l+xYnLJvhSw1/f0bOjpaiiZO3f3AW8hC+LbYkxZrOI/FhErrXMngDyRaQcuBu4z+q7GVgCbAFWAHcaYzzhxrTGuhe42xor3xo75jm6ekOUdvyVHHT1rvdwoqGF+5//kC/+qZSC3AxevOtirgpVcDUaCa591yVOg9p3gUwtymNgjotV23q2KNkqyGqMWQ4sD2r7QcDrJuCmMH0fAB6wM6bVXoGVPRfUHvMcAddXA6vDXVdCo2LUe/B6DX9fv4+f/3sbNY2tfOni0Xzrqgn0SQ9RXVvpkTgcwqXjC1i5rRKP1+B0hNpuT320SrgSFq/RpPCejjGGl7cc4cGXt7PtcB0lI/vz47lTmVx0iofCaaSUklw+oZDnNxzg/X0nmD6yf/QOKYiWGVLColLUczHGsKqskrmPvsW8v6ynqdXDb24+m7//9wWnLkigVcJTlEuLC3AIrC7ruanhGikpYdHadz2PZreHpe8f5Ik3d7HtcB1D+2XyixvP5PpzhnbvIXDLbFbtsmsXDyLNnUy/4kheVjrTRvRnVVkl86+ckGx3uoSKkhIV1aTU50RDC0+t2cuf3t5NVV0zEwfn8ssbz2Tu2UNxpcVhQaSurnvt4kGkuZPpV5y5fGIhv3ypjMq6Jgpz+yTbnZjR5TtF6aEYYyjdXc3dz7zPeT97lV++VMakIX35y+0z+PfXL+GmkuHxESSAefO61y4eRJo7mX7Fmcsm+AoAvNZDC7RqpKRERZfvUouahlaef28/f1u7l+1HTpKbkcYnS4bzmfNHMmFwgioVrI9aKCU2u3gQae5k+hVnJg/pS2FuBqu3V3FTyfDoHVIMjZQUpQfg9RrW7qpm/pKNzPjZK/zoxS1kpjtZeMMZrPnuLH5y3dTECRLA/PndaxcPIs2dTL/ijIhw2YQCXt9ehdsTsexnSqKRkhIVPQ49eew+Ws/z7x3gH+/tZ191I9kuJzdMH8atM0YwdWhe8hwrslkjz65dPIg0dzL9SgCXTShkSel+Nuw9wYzRA5LtTkyoKClR0eW7xFLT2Mq/PjjE8xv2U7rnOCJw0diBfPOj45k9dTBZrhT4s9VIKaW5uHggToewuqxSRUlRlNhp9Xh5Y0cVz204wMtbjtDi9jKuMId7Z0/kunOKGJKXmWwXO1JUZO/0Vrt28SDS3Mn0KwH07ZNOycj+rCqr4p7ZE5PtTkyoKClR0UApPni9hnW7q3nxg4Ms//Aw1fUt9M9K59YZI7h+2lDOGJqH7xzLFKS0tHvt4kGkuZPpV4K4bEIhC1ds43BNE4Pzek5quIqSoiQQYwwb99fw4saD/OuDQxyubaJPuoOPThrEtWcVcdmEwvilcXcn27fb25exaxcPIs2dTL8SxOUTC1i4Yhuvba/kU+eOSLY7tlFRUqKih/ydGsYYth2u48WNB3nxg4Psq27E5XTwkQkFfOesScyaWEh2Rg/7U1ywAFav7j67eBBp7mT6lSAmDMplSF4fVm2rUlFSFAUqqk7y4sZDvPjBQcorT+J0CBeNG8jXZhZz5ZTB5GWmJ9vFrmP3Cz2ZX/yR5u7lggTtqeEvbjxEq8dLeneWmYojPcNLRekBGGPYfqSOh1/ZzuyHX2fmr1/j4Ve3MyDbxU+vm8ra78ziyS/O4KaS4T1bkEArOvQQLi0u4GSzm437TiTbFdtopKSEJHDJTlfvwmOMYdOBWv696RArNh+moqoeEZg+oj/f+/gkPn7mkNTLnOsOSkq61y4eRJo7mX4lkPPH5CMC7+w8RsmonpEarqKkhMSrQhQWr9fw3r4TrNh0iH9vOsz+4404HcJ5owfwhQtHcdWUwRT27TnZTl1CI6UeQf9sF5MG9+Xtncf46qziZLtjCxUlJSStAeVJtKIDeKwyPyusiOhIbTPpTt8e0VdnjuOKyYMZkO1KtpuJIzfXXqVtu3bxINLcyfQrwVwwNp+/vLuHplZPjzhpWEVJCYnbq8t39c1u3iw/ysqtlbyy9QjH6lvISHPwkfEFXH3GYGZOHNTz94a6it0HT5P5gGqkuXvxg7PBXDg2nyfe3MWGvce5cOzAZLsTFRWlXs5f3t3DGUPzOHt4v5j69cRCjt3B3mMNrNx2hFe3VbKmopoWj5fcjDQ+MqGAq6cO4bIJBT0vfTserF4N11zTfXbxINLcyfQrwcwYPQCnQ3hn57EeIUq2su9EZLaIlIlIuYjcF+J6hog8Y11fIyKjAq7db7WXichV0cYUkdHWGDusMV1dmUNEhovIKhHZKiKbReTrsd+ens/3/7mJ6x59K+Z+ja2ette9OVJye7ysqTjG/yzfykcffI1Lf7mKBS9u4cCJRm67cCR/u+N8NvzgCh65dRofP3OICpKfRYu61y4eRJo7mX4lmNw+6Uwdmsc7O48l2xVbRP0LExEn8ChwBbAfWCciS40xWwLMbgeOG2PGicjNwELgUyIyGbgZmAIUAa+IyHirT7gxFwIPGWMWi8jvrbF/14U53MB8Y8wGEckF1ovIy0F+92q8p5Ct0NjiiW7UQzle38Jr26t4dVslr5VVUtvkJt0pnDc6n1tnjGDmxEJGDcxOtpupzYsvdq9dPIg0dzL9SgIXjs3n/16voL7ZnfI/rOxESjOAcmNMhTGmBVgMzA2ymQv82Xr9LDBLfEW75gKLjTHNxphdQLk1XsgxrT4zrTGwxryuK3MYYw4ZYzYAGGPqgK3AUHu3pXfQ0Np1YekQKXWHM0nEGEPZ4ToeW13OTb9/m+k/fZlvPPM+7+w8ypVTBvO7T09jw/ev4K9fOo8vXjxaBckOt97avXbxINLcyfQrCVwwJh+3VWsx1bEjmUOBfQHv9wPnhbMxxrhFpAbIt9rfDerrF4ZQY+YDJ4wx7hD2XZkDAGup7xxgTagPKCLzgHkAI0b0nHIc0TjZ5I5uFIbASKknlhlqavXwTsUxVm2r5NWtlRw40QjAlKK+3HX5OGZOGsSZQ/NwOFK04GmqM2dO99rFg0hzJ9OvJFAyqj/pTuGdimNcNqEw2e5ExI4ohfqrDf6WCmcTrj1UhBbJvitz+DqJ5ADPAd8wxtSGsMUYswhYBFBSUtLzvoHDUN9yCqJ0ClFWsjhc08SqMp8IvVV+lMZWD5npTi4aN5C7Zo7j8gmFPapackqjkVKPIsuVxtnD+/WIfSU7orQfCDzofRgQnE/pt9kvImlAHlAdpW+o9qNAPxFJs6KlQPuY5xCRdHyC9JQx5nkbn7VX0dDcdWFpaEn95Tuv1/DBgRpWbvVly20+6PvNMbRfJjeVDGPmxELOH5PfI57N6HGI2MuAsWsXDyLNnUy/ksR5o/P53Ws7U35fyY5n64BiERkNHMCXVBD8M2MpcBvwDnAjsNIYY0RkKfC0iDyILwmhGFiLL7rpNKbVZ5U1xmJrzBe6Moe13/QEsNUY82CsN6Y3cLK565HSiYaWbvSk+6hrauXNHUd5dVslq8sqOXqyBYfA9JH9uXf2RGZOLGT8oJzUPYeot2D3Cz2ZX/yR5j7NBAng3NEDeGRVORv2HueS4oJkuxOWqKJk7d/cBbwEOIE/GGM2i8iPgVJjzFJ8X/5/EZFyfNHLzVbfzSKyBNiCLxvuTmOMByDUmNaU9wKLReSnwHvW2MQ6h4hcDHwW+FBE3rfG+I4xZnnXblXPo8FavnN1oTpwdX1r+5sk//3uPlrPq9sqWbntCGt3VdPqMfTtk8ZHJhQya2IhHxlfQP/TqZpCKvD00/aWwOzaxYNIcyfTryQxbUQ/HALrdlX3bFECsL7Ilwe1/SDgdRNwU5i+DwAP2BnTaq/Al50X3B7THMaYNwm933Ta4I+UunJo3PGASCnRZYZaPV7W7a5m5dZKVpZVUlFVD8C4why+eNFoZk4sZPrI/qT1kFL8vZJly+x9qdu1iweR5k6mX0kit086U4ryWLMrtTPwUndhUTll/KKU0QVROlLbhNMheBJUmfXYyWZWl1WxsqyS18uqqGt243I6OG/MAD53/khmThzEiPyshPii2ODpp7vXLh5EmjuZfiWRc0cN4Kk1e2h2e8hIS829VhWlXsyxk75op28X6rPtPlrPqPwsdlbVx2X53X8a68ptlby69Qjv7TuBMVCQm8HHzhjCzEmFXDxuYEpvyJ7WXHONvQdQ7drFg0hzJ9OvJDJjdH/+8NYuNh2oYfrI1DzKQv/iezFVdc1d6uf1Giqq6jlvTD47raWz7qCp1cPbO4/y6tZKVm6r5FBNEwBnDsvj67OKmTmxkKlF+uxQj0CPruiR+M9UWrvruIqSknj8X/qtMRZXLTtSR12zm3NG9OOVrUdOaUepsq7JqrJdyZvlVTS1eslyObl43EC+8dFiLp9Q2PvPHuqNXHZZ99rFg0hzJ9OvJDIwJ4OxBdms213NlxmbbHdCoqLUi9l+xHdejNsTm6y8vOUI4KuXBbFlz/qOBD/JK1uP8PKWI7xvHcM8tF8mnywZzkcnDeK8MQNSdj1bsUlRkb3ziOzaxYNIcyfTryQzY/QAln1wCI/X4EzBVQkVpV7K4Zom9lY3AOD22o+U6pvd/PXdPVw8biBF/ewd4+3fH1r2wUGWfXCIPcd88541LI/5V4zno5MHMXFwrj471Juw+4WezC/+SHOfpoIEvmSHv63dR9nhOiYX9U22O51QUeql/N8bFQBcUjyQjVa0Eg2v1/C9f26isq6Z331mels+vTdMqFRV18yS0n38470DlFeexOkQLhybz7xLx/DRSYMYpMtyvZdFi+zty9i1iweR5k6mX0nmXGtfad3u6pQUJX3QoxdhjOHdimPc+n/v8sSbu/j0eSOYMCi3wymy4ahtauXOpzfwj/cOMP+K8Uwf2Z906zmg4AP/th+p42t/e48Lf/4qv3ypjAFZLn5y3VTWfGcWf7n9PD593kgVpN5OaWn32sWDSHMn068kM6x/JkV5fVibohXDNVLqBRhjWF1WxaOryindc5yBOS6+9/FJfOGi0fzqP2UR95QaWzw8+c5ufv/actzsTwAADURJREFUTmqb3Hz3Y5P40iWjgfaHbpvdPlGqa2rlgX9t5ZnSfWS70vjM+SP59HkjGVeYE/fPqKQYeshfj0VEOHf0AN7eeQxjTMotq6so9WA8XsOKTYd5bHU5mw/WUpTXhx9dO4VPnTu8rQhpukNo9Xo7/cfX7PaweO0+HllVTlVdM5eOL+DbV07gjGF5bTYZAaK0r7qBzzyxhn3VDdx+0WjuvHyclvY5nbnsMt+R4t1lFw8izZ1Mv1KAc0cN4IX3D7LnWEPKnR+motQDaXF7+ed7B/j9azupOFrP6IHZ/OLGM7nu7KGdSgqlOR0Y4xOwNKevQsNzG/bzm1d2cOBEIzNGD+DRW6cxY3TnZxbSnA6cDuF4Qwuf+8Najte38Mx/XdC2Jq2cxixY0L128SDS3Mn0KwXw/72v3V2toqR0nYqqkzyzbh/Prt/PsfoWphT15dFbpzF76uCwqZ1t+0Jew4a91Xzvnx+y/chJzhqWx89vOIOLxw2MGL67nA7++NZuRGDxHeerICk+xo/vXrt4EGnuZPqVAowryKF/VjprKqr5ZMnw6B0SiIpSilNZ28RLW46wbONB1uyqxukQZk0s5DPnj+SS4siCApDu9F3/67t7+NnyrQzrn8XvPu0TMjtryU1u37lKc88q4rwx+af+gZTeQUkJHAw+Vu0U7OJBpLmT6VcK4HAI54/J592K1NtXUlFKQfYea2DF5kOs2HSYDXt96dxjBmbz7asmcNP0YTFVQEizIqif/msrl44v4LFPTyMnhnpy/mzwOy4dY/8DKL0fu1/oyfzijzT3aSxIfi4Ym8+/Nx1mX3VjShU7VlFKAfxVEFZsOsyKzYfZesh3guqUor7Mv2I8s6cOZlxh1w6uSw/YY/p/t5wTkyABfH/OZA4cb2RKUV50Y+X04de/hvnzu88uHkSaO5l+pQgXWCsf71QcZUT+iCR7046KUpIwxrBxfw0rNh3mpc2H2XW0HhEoGdmf7318EldNGczwAaf+62X8oFwAvnLZWPK6UC389otHn7IPSi9EI6Uez7jCHAbmZPDOzmN86tzUESUxp+GxwJEoKSkxpXF6sM7jNazfc5zlHx7ipc2HOVTTRJpDuGBsPrOnDuaKyYMozO3+h063Ha5lXEGOHoqnKEoH7np6A+t2V/Pu/bNOeV9JRNYbY0pO1SeNlBLAjiN1PLNuHy9sPEhVXTOuNAeXFhfw7asmMGviIPKyYo9gYmHi4NQrJaL0cKZPh/Xru88uHkSaO5l+pRAXjM1n2QeH2HW0njEFqfEQvIpSnPB6Da9sPcLvX9vJhr0nSHMIsyYVMufMIi6fWBjz3o6ipBRa0aFX0L6vdExFqTdTurua77+wma2Hahk+IJPvfmwSn5g2lIE5Gcl2TVG6h9zc7rWLB5HmTqZfKcTogdkM7ZfJqm2VfPq8kcl2B7BZkFVEZotImYiUi8h9Ia5niMgz1vU1IjIq4Nr9VnuZiFwVbUwRGW2NscMa09Xdc8QLt8fLz5Zv5cbfv0NNQwsPfvIsVs2/jDsuHaOCpPQu5szpXrt4EGnuZPqVQogIV00ZzOs7jlLX1JpsdwAboiQiTuBR4GpgMnCLiEwOMrsdOG6MGQc8BCy0+k4GbgamALOBx0TEGWXMhcBDxphi4Lg1dnfP0e20erx8+akNLHq9gk+fN4KX7/4I108bpskFSu9k+/butYsHkeZOpl8pxtVnDKbF7eXVrZXJdgWwFynNAMqNMRXGmBZgMTA3yGYu8Gfr9bPALPGlcswFFhtjmo0xu4Bya7yQY1p9ZlpjYI15XXfOYe+2xIbXa7jn2Q94ecsRFlwzmQc+cQbZumek9Ga09l2vYfqI/owemM1vXt3B0ZPNyXbH1p7SUGBfwPv9wHnhbIwxbhGpAfKt9neD+g61XocaMx84YYxxh7Dvrjk6ISLzAP+JXydF5BhwNJRtNL6wEL7QlY6pyUC6eB96IXov2vHdix/9yJ61Xbt4EGnuU/er1/03UfDtLncdCHTLppQdUQqVvB78cFM4m3DtoSK0SPbdOUfnRmMWAW3pOCJS2h359j0dvQ/t6L1oR++FD70P7Vj3YlR3jGVn+W4/EFhGdhgQ/Dh0m42IpAF5QHWEvuHajwL9rDGC5+quORRFUZQUxY4orQOKraw4F76kgqVBNkuB26zXNwIrja9UxFLgZitzbjRQDKwNN6bVZ5U1BtaYL3TnHPZui6IoipIMoi7fWfs3dwEvAU7gD8aYzSLyY6DUGLMUeAL4i4iU44tebrb6bhaRJcAWwA3caYzxAIQa05ryXmCxiPwUeM8am26eIxr6ZJ0PvQ/t6L1oR++FD70P7XTbvdDad4qiKErKoA/RKIqiKCmDipKiKIqSMqgoBZDoskTJQET+ICKVIrIpoG2AiLxslXZ6WUT6W+0iIr+17scHIjItoM9tlv0OEbkt1FypjIgMF5FVIrJVRDaLyNet9tPxXvQRkbUistG6Fz+y2rut5FdPwqoI856ILLPen673YbeIfCgi74tIqdUW/78PY4z+z7ev5gR2AmMAF7ARmJxsv+LwOS8FpgGbAtp+Adxnvb4PWGi9/hjwb3zPgp0PrLHaBwAV1r/9rdf9k/3ZYrwPQ4Bp1utcYDu+clSn470QIMd6nQ6ssT7jEuBmq/33wJet118Bfm+9vhl4xno92fq7yQBGW39PzmR/vi7cj7uBp4Fl1vvT9T7sBv5/e/cWKlUVx3H8+6PEThesMCNSKKnU4ngBUcEgCVG6EBJBgSBldIHsQoERRvhaQYS9lfZgRA9hN4IoCcygetAwLU/a0aQHD56yMgOpyH8Pa+1mczhzmpNzZvbM/D4wzOx19uw9e3H2/GetvfZ/TR1RNuHnh1tKNS1LS9ROEbGTNHqxrJzCaWRqp62RfEG6h+wyYCWwPSJ+johfgO2kvIMdIyKGIuLL/PokMEDKBNKLdRER8XtenJQfQfNSfnUMSdOBW4DNebmZqc+6wYSfHw5KNaOlU7q8zrrd5tKIGIL0ZQ1My+X16qSr6ip3uywgtRB6si5yl9UeYJj0xXGIBlN+AeWUX51eFy8C64HTebnh1Gd0Vz1A+mHykaTdSqnYoAXnh7OG1jSSTqnXjDe1U8eRdD6wDXgsIn5T/Smhu7ouIt3bN1/ShcDbwJzRVsvPXVkXkm4FhiNit6RlRfEoq3Z1PZQsjYijkqYB2yV9O8a6TasLt5Rqejkt0bHc1CY/FznsuzqFk6RJpID0ekS8lYt7si4KEfErsIN0XaBZKb86xVLgNklHSN33N5JaTr1WDwBExNH8PEz6obKIFpwfDko1vZyWqJzCaWRqpzV5ZM0S4ERusn8IrJB0UR59syKXdYzc978FGIiIF0p/6sW6uCS3kJDUBywnXWNrVsqvjhART0XE9EiJRe8iHddqeqweACSdJ+mC4jXp//prWnF+tHuER5UepBEkB0n96Rva/Xkm6BjfAIaAv0i/Yu4l9YN/DHyXny/O64o0UeIhYB+wsLSdtaQLuIPAPe0+rv9RD9eTuhH2Anvy4+YerYu5pJRee/MXzzO5fCbpy3QQeBOYnMvPycuD+e8zS9vakOvoAHBTu4/tDOpkGbXRdz1XD/mYv8qPb4rvw1acH04zZGZmleHuOzMzqwwHJTMzqwwHJTMzqwwHJTMzqwwHJTMzqwwHJTMzqwwHJbMWkLRK0iuS3pW0olf2bTZeDkpmLRAR70TEfcDdwJ2QEsFKOpUToZLLHpAUkm4ola3LZcvrbV/SVZL2jSibLOl74OAo++7L8+T8KWlqEw/V7Iw4KJm11tOkO98LhyJifml5LimzwhwASeeSsm78SLpTvp7DwAxJ5XP6fuCTiNg/ct8RcSrvt+Nysll3c1AyazJJV0vaIWmXpOfybJyS9CzwQeR5nOroJ6WCmp2XHyGlsjkdEcfy9q/MXXG7lGaMnRURp4EfgCvyOn3AE8DGcezbrO0clMyaSNJZwFbg8YhYCPSRcoc9TEp0eoekB8fYxBzSTKezJU0hdbd9RspJV2Q231za/kbSDKCQkqgWwewh4L2IODKOfZu1nedTMmuuVcD+UotkgDRJ3CZg01hvlDQDOB4Rh/McNuuBl4BrSF16xfavA7bluZ/OBj4t7WuWpJ2koLQEoJF9m1WFg5JZcy0gZRwvzCPN5NqIudSuG50kTRu9iDSnTxHk5pEyNm8Z5f0DpDmAHiXNEXVsfB/drP3cfWfWXMfJXWiSFgNrqLVy/ks/taD0PLAu0oyw/aVtDAEriwENkvpVmy53gBTE1ub3m3UcByWz5noNWJiHZ99OClKDDb63n3ztKCLej4jPc/m1QDGC7lXSeTuQh5I/GbX5Zw7kbbwcESfO+EjM2sDdd2ZNFBE/AYvh32tEy/LIuEbeu7pO+bTS61PUZkEdud4f+Jy2DueWktnEmcfYXXd/A1PKN8+2SnHzLDAJaChomrWCZ541M7PKcEvJzMwqw0HJzMwqw0HJzMwqw0HJzMwqw0HJzMwqw0HJzMwqw0HJzMwqw0HJzMwqw0HJzMwq4x9J2syy0xG2/AAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.clf()\n",
"# plt.plot(x_part, calcs, '.')\n",
"plt.plot(test_q, calcs_test)#, label = 'pdf (Ctt = 0.0)')\n",
"# plt.plot(test_q, calcs_test1, label = 'pdf (Ctt = 0.5)')\n",
"# plt.plot(test_q, calcs_test2, label = 'pdf (D-contribs = 0.3)')\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.axvline(x=jpsi_mass -70,color='red', linewidth=1.0, linestyle = 'dotted')\n",
"plt.axvline(x=jpsi_mass +70,color='red', linewidth=1.0, linestyle = 'dotted')\n",
"plt.axvline(x=psi2s_mass -50,color='red', linewidth=1.0, linestyle = 'dotted')\n",
"plt.axvline(x=psi2s_mass +50,color='red', linewidth=1.0, linestyle = 'dotted')\n",
"# plt.legend()\n",
"plt.ylim(0.0, 1.5e-6)\n",
"plt.xlabel(r'$q^2 [MeV^2]$')\n",
"# plt.yscale('log')\n",
"# plt.xlim(770, 785)\n",
"plt.savefig('test.png')\n",
"# print(jpsi_width)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"\n",
"\n",
"# probs = mixture.prob(test_q)\n",
"# probs_np = zfit.run(probs)\n",
"# probs_np *= np.max(calcs_test) / np.max(probs_np)\n",
"# plt.figure()\n",
"# plt.semilogy(test_q, probs_np,label=\"importance sampling\")\n",
"# plt.semilogy(test_q, calcs_test, label = 'pdf')\n"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"# 0.213/(0.00133+0.213+0.015)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Adjust scaling of different parts"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"total_f.update_integration_options(draws_per_dim=2000000, mc_sampler=None)\n",
"# total_f_fit.update_integration_options(draws_per_dim=2000000, mc_sampler=None)\n",
"# inte = total_f.integrate(limits = (950., 1050.), norm_range=False)\n",
"# inte_fl = zfit.run(inte)\n",
"# print(inte_fl/4500)\n",
"# print(pdg[\"jpsi_BR\"]/pdg[\"NR_BR\"], inte_fl*pdg[\"psi2s_auc\"]/pdg[\"NR_auc\"])"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"# # print(\"jpsi:\", inte_fl)\n",
"# # print(\"Increase am by factor:\", np.sqrt(pdg[\"jpsi_BR\"]/pdg[\"NR_BR\"]*pdg[\"NR_auc\"]/inte_fl))\n",
"# # print(\"New amp:\", pdg[\"jpsi\"][3]*np.sqrt(pdg[\"jpsi_BR\"]/pdg[\"NR_BR\"]*pdg[\"NR_auc\"]/inte_fl))\n",
"\n",
"# # print(\"psi2s:\", inte_fl)\n",
"# # print(\"Increase am by factor:\", np.sqrt(pdg[\"psi2s_BR\"]/pdg[\"NR_BR\"]*pdg[\"NR_auc\"]/inte_fl))\n",
"# # print(\"New amp:\", pdg[\"psi2s\"][3]*np.sqrt(pdg[\"psi2s_BR\"]/pdg[\"NR_BR\"]*pdg[\"NR_auc\"]/inte_fl))\n",
"\n",
"# name = \"phi\"\n",
"\n",
"# print(name+\":\", inte_fl)\n",
"# print(\"Increase am by factor:\", np.sqrt(pdg[name+\"_BR\"]/pdg[\"NR_BR\"]*pdg[\"NR_auc\"]/inte_fl))\n",
"# print(\"New amp:\", pdg[name][0]*np.sqrt(pdg[name+\"_BR\"]/pdg[\"NR_BR\"]*pdg[\"NR_auc\"]/inte_fl))\n",
"\n",
"\n",
"# print(x_min)\n",
"# print(x_max)\n",
"# # total_f.update_integration_options(draws_per_dim=2000000, mc_sampler=None)\n",
"# total_f.update_integration_options(mc_sampler=lambda dim, num_results,\n",
"# dtype: tf.random_uniform(maxval=1., shape=(num_results, dim), dtype=dtype),\n",
"# draws_per_dim=1000000)\n",
"# # _ = []\n",
"\n",
"# # for i in range(10):\n",
"\n",
"# # inte = total_f.integrate(limits = (x_min, x_max))\n",
"# # inte_fl = zfit.run(inte)\n",
"# # print(inte_fl)\n",
"# # _.append(inte_fl)\n",
"\n",
"# # print(\"mean:\", np.mean(_))\n",
"\n",
"# _ = time.time()\n",
"\n",
"# inte = total_f.integrate(limits = (x_min, x_max))\n",
"# inte_fl = zfit.run(inte)\n",
"# print(inte_fl)\n",
"# print(\"Time taken: {}\".format(display_time(int(time.time() - _))))\n",
"\n",
"# print(pdg['NR_BR']/pdg['NR_auc']*inte_fl)\n",
"# print(0.25**2*4.2/1000)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Sampling\n",
"## Mixture distribution for sampling"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"\n",
" \n",
"# print(list_of_borders[:9])\n",
"# print(list_of_borders[-9:])\n",
"\n",
"\n",
"class UniformSampleAndWeights(zfit.util.execution.SessionHolderMixin):\n",
" def __call__(self, limits, dtype, n_to_produce):\n",
" # n_to_produce = tf.cast(n_to_produce, dtype=tf.int32)\n",
" low, high = limits.limit1d\n",
" low = tf.cast(low, dtype=dtype)\n",
" high = tf.cast(high, dtype=dtype)\n",
"# uniform = tfd.Uniform(low=low, high=high)\n",
"# uniformjpsi = tfd.Uniform(low=tf.constant(3080, dtype=dtype), high=tf.constant(3112, dtype=dtype))\n",
"# uniformpsi2s = tfd.Uniform(low=tf.constant(3670, dtype=dtype), high=tf.constant(3702, dtype=dtype))\n",
"\n",
"# list_of_borders = []\n",
"# _p = []\n",
"# splits = 10\n",
"\n",
"# _ = np.linspace(x_min, x_max, splits)\n",
"\n",
"# for i in range(splits):\n",
"# list_of_borders.append(tf.constant(_[i], dtype=dtype))\n",
"# _p.append(tf.constant(1/splits, dtype=dtype))\n",
" \n",
"# mixture = tfd.MixtureSameFamily(mixture_distribution=tfd.Categorical(probs=_p[:(splits-1)]),\n",
"# components_distribution=tfd.Uniform(low=list_of_borders[:(splits-1)], \n",
"# high=list_of_borders[-(splits-1):]))\n",
" mixture = tfd.MixtureSameFamily(mixture_distribution=tfd.Categorical(probs=[tf.constant(0.05, dtype=dtype),\n",
" tf.constant(0.93, dtype=dtype),\n",
" tf.constant(0.05, dtype=dtype),\n",
" tf.constant(0.065, dtype=dtype),\n",
" tf.constant(0.04, dtype=dtype),\n",
" tf.constant(0.05, dtype=dtype)]),\n",
" components_distribution=tfd.Uniform(low=[tf.constant(x_min, dtype=dtype), \n",
" tf.constant(3090, dtype=dtype),\n",
" tf.constant(3681, dtype=dtype), \n",
" tf.constant(3070, dtype=dtype),\n",
" tf.constant(1000, dtype=dtype),\n",
" tf.constant(3660, dtype=dtype)], \n",
" high=[tf.constant(x_max, dtype=dtype),\n",
" tf.constant(3102, dtype=dtype), \n",
" tf.constant(3691, dtype=dtype),\n",
" tf.constant(3110, dtype=dtype),\n",
" tf.constant(1040, dtype=dtype),\n",
" tf.constant(3710, dtype=dtype)]))\n",
"# dtype = tf.float64\n",
"# mixture = tfd.MixtureSameFamily(mixture_distribution=tfd.Categorical(probs=[tf.constant(0.04, dtype=dtype),\n",
"# tf.constant(0.90, dtype=dtype),\n",
"# tf.constant(0.02, dtype=dtype),\n",
"# tf.constant(0.07, dtype=dtype),\n",
"# tf.constant(0.02, dtype=dtype)]),\n",
"# components_distribution=tfd.Uniform(low=[tf.constant(x_min, dtype=dtype), \n",
"# tf.constant(3089, dtype=dtype),\n",
"# tf.constant(3103, dtype=dtype), \n",
"# tf.constant(3681, dtype=dtype),\n",
"# tf.constant(3691, dtype=dtype)], \n",
"# high=[tf.constant(3089, dtype=dtype),\n",
"# tf.constant(3103, dtype=dtype), \n",
"# tf.constant(3681, dtype=dtype),\n",
"# tf.constant(3691, dtype=dtype), \n",
"# tf.constant(x_max, dtype=dtype)]))\n",
"# mixture = tfd.Uniform(tf.constant(x_min, dtype=dtype), tf.constant(x_max, dtype=dtype))\n",
"# sample = tf.random.uniform((n_to_produce, 1), dtype=dtype)\n",
" sample = mixture.sample((n_to_produce, 1))\n",
"# sample = tf.random.uniform((n_to_produce, 1), dtype=dtype)\n",
" weights = mixture.prob(sample)[:,0]\n",
"# weights = tf.broadcast_to(tf.constant(1., dtype=dtype), shape=(n_to_produce,))\n",
" # sample = tf.expand_dims(sample, axis=-1)\n",
"# print(sample, weights)\n",
" \n",
"# weights = tf.ones(shape=(n_to_produce,), dtype=dtype)\n",
" weights_max = None\n",
" thresholds = tf.random_uniform(shape=(n_to_produce,), dtype=dtype)\n",
" return sample, thresholds, weights, weights_max, n_to_produce"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"# total_f._sample_and_weights = UniformSampleAndWeights"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Constraints"
]
},
{
"cell_type": "code",
"execution_count": 20,
"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/3.)\n",
" sigma2 = ztf.constant(0.128/3.)\n",
" sigma3 = ztf.constant(0.548/3.)\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": 21,
"metadata": {},
"outputs": [],
"source": [
"param_values_dic = {\n",
" 'jpsi_m': jpsi_mass,\n",
" 'jpsi_s': jpsi_scale,\n",
" 'jpsi_p': jpsi_phase,\n",
" 'jpsi_w': jpsi_width,\n",
" 'psi2s_m': psi2s_mass,\n",
" 'psi2s_s': psi2s_scale,\n",
" 'psi2s_p': psi2s_phase,\n",
" 'psi2s_w': psi2s_width,\n",
" 'p3770_m': p3770_mass,\n",
" 'p3770_s': p3770_scale,\n",
" 'p3770_p': p3770_phase,\n",
" 'p3770_w': p3770_width,\n",
" 'p4040_m': p4040_mass,\n",
" 'p4040_s': p4040_scale,\n",
" 'p4040_p': p4040_phase,\n",
" 'p4040_w': p4040_width,\n",
" 'p4160_m': p4160_mass,\n",
" 'p4160_s': p4160_scale,\n",
" 'p4160_p': p4160_phase,\n",
" 'p4160_w': p4160_width,\n",
" 'p4415_m': p4415_mass,\n",
" 'p4415_s': p4415_scale,\n",
" 'p4415_p': p4415_phase,\n",
" 'p4415_w': p4415_width,\n",
" 'rho_m': rho_mass,\n",
" 'rho_s': rho_scale,\n",
" 'rho_p': rho_phase,\n",
" 'rho_w': rho_width,\n",
" 'omega_m': omega_mass,\n",
" 'omega_s': omega_scale,\n",
" 'omega_p': omega_phase,\n",
" 'omega_w': omega_width,\n",
" 'phi_m': phi_mass,\n",
" 'phi_s': phi_scale,\n",
" 'phi_p': phi_phase,\n",
" 'phi_w': phi_width,\n",
" 'Dstar_m': Dstar_mass,\n",
" 'DDstar_s': 0.0,\n",
" 'DDstar_p': 0.0,\n",
" 'D_m': D_mass,\n",
" 'Dbar_m': Dbar_mass,\n",
" 'Dbar_s': 0.0,\n",
" 'Dbar_p': 0.0,\n",
" 'tau_m': pdg['tau_M'],\n",
" 'Ctt': 0.0,\n",
" 'b0_0': 0.292,\n",
" 'b0_1': 0.281,\n",
" 'b0_2': 0.150,\n",
" 'bplus_0': 0.466,\n",
" 'bplus_1': -0.885,\n",
" 'bplus_2': -0.213,\n",
" 'bT_0': 0.460,\n",
" 'bT_1': -1.089,\n",
" 'bT_2': -1.114}\n",
"\n",
"\n",
"def reset_param_values(variation = 0.05):\n",
" for param in total_f_fit.get_dependents():\n",
" if param.floating:\n",
" param.set_value(param_values_dic[param.name] + random.uniform(-1, 1) * param_values_dic[param.name]* variation)\n",
"# print(param.name)\n",
"# for param in totalf.get_dependents():\n",
"# param.set_value()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Pull dictionary"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"# BR_steps = np.linspace(0.0, 1e-3, 11)\n",
"pull_dic = {}\n",
"\n",
"mi = 2e-4\n",
"ma = 6e-4\n",
"ste = 5\n",
"\n",
"# mi = 1e-4\n",
"# ma = 3e-3\n",
"# ste = 20\n",
"\n",
"for param in total_f_fit.get_dependents():\n",
" if param.floating:\n",
" pull_dic[param.name] = []\n",
" for step in range(2*ste):\n",
" pull_dic[param.name].append([])\n",
" \n",
"\n",
"\n",
"def save_pulls(step):\n",
" for param in total_f_fit.get_dependents():\n",
" if param.floating:\n",
" pull_dic[param.name][step].append((params[param]['value'] - param_values_dic[param.name])/params[param]['minuit_hesse']['error'])\n",
"\n",
"\n",
"\n",
"# for key in pull_dic.keys():\n",
"# print(np.shape(pull_dic[key]))\n",
"# save_pulls(New_step=True)\n",
"# params[Ctt]['value']"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"36668\n",
"5404696\n"
]
}
],
"source": [
"# for param in total_f_fit.get_dependents():\n",
"# if param.floating:\n",
"# print(param.name)\n",
"\n",
"# print(params[Ctt])\n",
"\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",
"nevents = int(pdg[\"number_of_decays\"]*cut_BR)\n",
"\n",
"print(nevents)\n",
"print(int(pdg[\"number_of_decays\"]))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# CLS Code"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[0. 0.0003 0.0004 0.0005 0.0006]\n",
"[0. 0.26726124 0.3086067 0.34503278 0.37796447]\n",
"WARNING:tensorflow:From C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\core\\sample.py:163: to_int32 (from tensorflow.python.ops.math_ops) is deprecated and will be removed in a future version.\n",
"Instructions for updating:\n",
"Use tf.cast instead.\n",
"Step: 0/5\n",
"Current Ctt: 0.0\n",
"Ctt floating: True\n",
"Toy 0/1 - Fit 0/1\n",
"------------------------------------------------------------------\n",
"| FCN = 2.977E+05 | Ncalls=1078 (1078 total) |\n",
"| EDM = 0.013 (Goal: 5E-06) | up = 0.5 |\n",
"------------------------------------------------------------------\n",
"| Valid Min. | Valid Param. | Above EDM | Reached call limit |\n",
"------------------------------------------------------------------\n",
"| False | True | True | False |\n",
"------------------------------------------------------------------\n",
"| Hesse failed | Has cov. | Accurate | Pos. def. | Forced |\n",
"------------------------------------------------------------------\n",
"| False | True | False | False | True |\n",
"------------------------------------------------------------------\n",
"Function minimum: 297652.2218218978\n",
"----------------------------------------------------------------------------------------------\n",
"| | Name | Value | Hesse Err | Minos Err- | Minos Err+ | Limit- | Limit+ | Fixed |\n",
"----------------------------------------------------------------------------------------------\n",
"| 0 | phi_s | 20.2 | 0.7 | | | 14.8182 | 23.5818 | |\n",
"| 1 | rho_p | -0.39 | 0.19 | | |-6.28319 | 6.28319 | |\n",
"| 2 | rho_s | 1.59 | 0.24 | | |0.0253049| 2.0747 | |\n",
"| 3 | p4415_p | -2.37 | 0.18 | | |-6.28319 | 6.28319 | |\n",
"| 4 | psi2s_p | 2.08 | 0.04 | | |-6.28319 | 6.28319 | |\n",
"| 5 | Dbar_p | 0.13 | 0.85 | | |-6.28319 | 6.28319 | |\n",
"| 6 | omega_p | 0.16 | 0.21 | | |-6.28319 | 6.28319 | |\n",
"| 7 | DDstar_p | -2.6 | 7.0 | | |-6.28319 | 6.28319 | |\n",
"| 8 | jpsi_p | -1.506 | 0.022 | | |-6.28319 | 6.28319 | |\n",
"| 9 | p3770_s | 2.43 | 0.21 | | |0.918861 | 4.08114 | |\n",
"| 10| p4415_s | 1.50 | 0.19 | | |0.126447 | 2.35355 | |\n",
"| 11| bplus_1 | -0.889 | 0.026 | | | -2 | 2 | |\n",
"| 12| bplus_2 | -0.24 | 0.04 | | | -2 | 2 | |\n",
"| 13| p3770_p | -2.79 | 0.12 | | |-6.28319 | 6.28319 | |\n",
"| 14| bplus_0 | 0.478 | 0.013 | | | -2 | 2 | |\n",
"| 15| p4160_p | -2.19 | 0.07 | | |-6.28319 | 6.28319 | |\n",
"| 16| p4160_s | 2.26 | 0.20 | | | 0.71676 | 3.68324 | |\n",
"| 17| Dbar_s | 0.16 | 0.15 | | | -0.3 | 0.3 | |\n",
"| 18| omega_s | 6.1 | 0.9 | | | 4.19232 | 9.40768 | |\n",
"| 19| phi_p | 0.75 | 0.14 | | |-6.28319 | 6.28319 | |\n",
"| 20| Ctt | 0.09 | 0.16 | | | -2.5 | 2.5 | |\n",
"| 21| p4040_p | -2.88 | 0.14 | | |-6.28319 | 6.28319 | |\n",
"| 22| DDstar_s | 0.04 | 0.18 | | | -0.3 | 0.3 | |\n",
"| 23| p4040_s | 1.28 | 0.20 | | |0.00501244| 2.01499 | |\n",
"----------------------------------------------------------------------------------------------\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"| | phi_s rho_p rho_s p4415_p psi2s_p Dbar_p omega_p DDstar_p jpsi_p p3770_s p4415_s bplus_1 bplus_2 p3770_p bplus_0 p4160_p p4160_s Dbar_s omega_s phi_p Ctt p4040_p DDstar_s p4040_s |\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"| phi_s | 1.000 -0.067 -0.002 -0.063 -0.056 0.057 0.005 0.072 0.003 -0.037 -0.043 0.059 -0.078 0.042 0.057 -0.023 -0.062 -0.067 -0.011 0.427 0.060 -0.055 -0.071 -0.058 |\n",
"| rho_p | -0.067 1.000 0.152 0.105 0.089 -0.085 -0.018 -0.117 -0.031 0.060 0.067 -0.089 0.171 -0.069 -0.088 0.038 0.102 0.109 0.199 -0.169 -0.091 0.088 0.115 0.094 |\n",
"| rho_s | -0.002 0.152 1.000 -0.187 -0.172 0.195 0.041 0.239 -0.052 -0.087 -0.144 0.225 -0.015 0.143 0.245 -0.050 -0.186 -0.218 -0.287 0.046 0.180 -0.150 -0.237 -0.182 |\n",
"| p4415_p | -0.063 0.105 -0.187 1.000 0.553 -0.599 0.050 -0.775 0.245 0.284 0.431 -0.722 0.040 -0.428 -0.712 0.267 0.562 0.700 0.075 -0.036 -0.562 0.511 0.748 0.544 |\n",
"| psi2s_p | -0.056 0.089 -0.172 0.553 1.000 -0.480 0.046 -0.719 0.254 0.148 0.444 -0.668 -0.062 -0.421 -0.661 0.153 0.537 0.638 0.070 -0.033 -0.479 0.415 0.695 0.561 |\n",
"| Dbar_p | 0.057 -0.085 0.195 -0.599 -0.480 1.000 -0.049 0.834 -0.024 -0.172 -0.526 0.780 0.002 0.649 0.772 -0.010 -0.649 -0.790 -0.077 0.040 0.646 -0.397 -0.784 -0.644 |\n",
"| omega_p | 0.005 -0.018 0.041 0.050 0.046 -0.049 1.000 -0.063 0.008 0.024 0.038 -0.057 0.014 -0.038 -0.059 0.014 0.049 0.057 0.390 -0.014 -0.049 0.040 0.062 0.048 |\n",
"| DDstar_p | 0.072 -0.117 0.239 -0.775 -0.719 0.834 -0.063 1.000 -0.286 -0.343 -0.624 0.935 0.016 0.601 0.924 -0.189 -0.781 -0.933 -0.097 0.044 0.746 -0.611 -0.960 -0.769 |\n",
"| jpsi_p | 0.003 -0.031 -0.052 0.245 0.254 -0.024 0.008 -0.286 1.000 0.172 0.152 -0.273 -0.026 -0.071 -0.276 0.161 0.200 0.281 0.015 -0.007 -0.112 0.225 0.239 0.207 |\n",
"| p3770_s | -0.037 0.060 -0.087 0.284 0.148 -0.172 0.024 -0.343 0.172 1.000 0.218 -0.313 -0.020 -0.316 -0.307 0.097 0.281 0.264 0.034 -0.017 -0.290 0.229 0.326 0.322 |\n",
"| p4415_s | -0.043 0.067 -0.144 0.431 0.444 -0.526 0.038 -0.624 0.152 0.218 1.000 -0.552 -0.085 -0.405 -0.546 0.136 0.585 0.591 0.058 -0.026 -0.378 0.411 0.619 0.509 |\n",
"| bplus_1 | 0.059 -0.089 0.225 -0.722 -0.668 0.780 -0.057 0.935 -0.273 -0.313 -0.552 1.000 -0.039 0.571 0.823 -0.184 -0.714 -0.849 -0.089 0.037 0.687 -0.567 -0.926 -0.704 |\n",
"| bplus_2 | -0.078 0.171 -0.015 0.040 -0.062 0.002 0.014 0.016 -0.026 -0.020 -0.085 -0.039 1.000 0.005 -0.040 0.050 0.010 0.034 0.010 -0.025 0.276 0.008 -0.017 -0.020 |\n",
"| p3770_p | 0.042 -0.069 0.143 -0.428 -0.421 0.649 -0.038 0.601 -0.071 -0.316 -0.405 0.571 0.005 1.000 0.566 0.026 -0.479 -0.559 -0.059 0.029 0.395 -0.269 -0.600 -0.480 |\n",
"| bplus_0 | 0.057 -0.088 0.245 -0.712 -0.661 0.772 -0.059 0.924 -0.276 -0.307 -0.546 0.823 -0.040 0.566 1.000 -0.179 -0.703 -0.839 -0.094 0.038 0.681 -0.558 -0.915 -0.695 |\n",
"| p4160_p | -0.023 0.038 -0.050 0.267 0.153 -0.010 0.014 -0.189 0.161 0.097 0.136 -0.184 0.050 0.026 -0.179 1.000 0.077 0.129 0.019 -0.009 -0.242 0.131 0.148 -0.050 |\n",
"| p4160_s | -0.062 0.102 -0.186 0.562 0.537 -0.649 0.049 -0.781 0.200 0.281 0.585 -0.714 0.010 -0.479 -0.703 0.077 1.000 0.720 0.075 -0.036 -0.503 0.584 0.767 0.552 |\n",
"| Dbar_s | -0.067 0.109 -0.218 0.700 0.638 -0.790 0.057 -0.933 0.281 0.264 0.591 -0.849 0.034 -0.559 -0.839 0.129 0.720 1.000 0.088 -0.040 -0.743 0.521 0.916 0.723 |\n",
"| omega_s | -0.011 0.199 -0.287 0.075 0.070 -0.077 0.390 -0.097 0.015 0.034 0.058 -0.089 0.010 -0.059 -0.094 0.019 0.075 0.088 1.000 -0.072 -0.074 0.060 0.096 0.074 |\n",
"| phi_p | 0.427 -0.169 0.046 -0.036 -0.033 0.040 -0.014 0.044 -0.007 -0.017 -0.026 0.037 -0.025 0.029 0.038 -0.009 -0.036 -0.040 -0.072 1.000 0.038 -0.029 -0.045 -0.034 |\n",
"| Ctt | 0.060 -0.091 0.180 -0.562 -0.479 0.646 -0.049 0.746 -0.112 -0.290 -0.378 0.687 0.276 0.395 0.681 -0.242 -0.503 -0.743 -0.074 0.038 1.000 -0.483 -0.732 -0.482 |\n",
"| p4040_p | -0.055 0.088 -0.150 0.511 0.415 -0.397 0.040 -0.611 0.225 0.229 0.411 -0.567 0.008 -0.269 -0.558 0.131 0.584 0.521 0.060 -0.029 -0.483 1.000 0.581 0.401 |\n",
"| DDstar_s | -0.071 0.115 -0.237 0.748 0.695 -0.784 0.062 -0.960 0.239 0.326 0.619 -0.926 -0.017 -0.600 -0.915 0.148 0.767 0.916 0.096 -0.045 -0.732 0.581 1.000 0.753 |\n",
"| p4040_s | -0.058 0.094 -0.182 0.544 0.561 -0.644 0.048 -0.769 0.207 0.322 0.509 -0.704 -0.020 -0.480 -0.695 -0.050 0.552 0.723 0.074 -0.034 -0.482 0.401 0.753 1.000 |\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"Hesse errors: OrderedDict([(<zfit.Parameter 'phi_s' floating=True>, {'error': 0.7471080388438764}), (<zfit.Parameter 'rho_p' floating=True>, {'error': 0.18514176261667714}), (<zfit.Parameter 'rho_s' floating=True>, {'error': 0.24325558408812809}), (<zfit.Parameter 'p4415_p' floating=True>, {'error': 0.17951439955418258}), (<zfit.Parameter 'psi2s_p' floating=True>, {'error': 0.036277261386909565}), (<zfit.Parameter 'Dbar_p' floating=True>, {'error': 0.8521601760503357}), (<zfit.Parameter 'omega_p' floating=True>, {'error': 0.20697201816491884}), (<zfit.Parameter 'DDstar_p' floating=True>, {'error': 6.981778412516394}), (<zfit.Parameter 'jpsi_p' floating=True>, {'error': 0.022270745277642945}), (<zfit.Parameter 'p3770_s' floating=True>, {'error': 0.20792441648302717}), (<zfit.Parameter 'p4415_s' floating=True>, {'error': 0.18559764903708653}), (<zfit.Parameter 'bplus_1' floating=True>, {'error': 0.02634389474470422}), (<zfit.Parameter 'bplus_2' floating=True>, {'error': 0.04215911524062499}), (<zfit.Parameter 'p3770_p' floating=True>, {'error': 0.12126027556157148}), (<zfit.Parameter 'bplus_0' floating=True>, {'error': 0.01325866696713196}), (<zfit.Parameter 'p4160_p' floating=True>, {'error': 0.07358429037598313}), (<zfit.Parameter 'p4160_s' floating=True>, {'error': 0.20004363541048886}), (<zfit.Parameter 'Dbar_s' floating=True>, {'error': 0.14726695003307771}), (<zfit.Parameter 'omega_s' floating=True>, {'error': 0.8574758160021965}), (<zfit.Parameter 'phi_p' floating=True>, {'error': 0.1386560878188816}), (<zfit.Parameter 'Ctt' floating=True>, {'error': 0.15908203934676957}), (<zfit.Parameter 'p4040_p' floating=True>, {'error': 0.13998911285038518}), (<zfit.Parameter 'DDstar_s' floating=True>, {'error': 0.17548878917420122}), (<zfit.Parameter 'p4040_s' floating=True>, {'error': 0.19920259043511457})])\n",
"Step: 0/5\n",
"Current Ctt: 0.0\n",
"Ctt floating: True\n",
"Toy 0/1 - Fit 0/1\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"------------------------------------------------------------------\n",
"| FCN = 2.977E+05 | Ncalls=1119 (1119 total) |\n",
"| EDM = 0.000952 (Goal: 5E-06) | up = 0.5 |\n",
"------------------------------------------------------------------\n",
"| Valid Min. | Valid Param. | Above EDM | Reached call limit |\n",
"------------------------------------------------------------------\n",
"| True | True | False | False |\n",
"------------------------------------------------------------------\n",
"| Hesse failed | Has cov. | Accurate | Pos. def. | Forced |\n",
"------------------------------------------------------------------\n",
"| False | True | True | True | False |\n",
"------------------------------------------------------------------\n",
"Function minimum: 297676.8602903808\n",
"----------------------------------------------------------------------------------------------\n",
"| | Name | Value | Hesse Err | Minos Err- | Minos Err+ | Limit- | Limit+ | Fixed |\n",
"----------------------------------------------------------------------------------------------\n",
"| 0 | phi_s | 20.6 | 0.9 | | | 14.8182 | 23.5818 | |\n",
"| 1 | rho_p | -0.3 | 0.7 | | |-6.28319 | 6.28319 | |\n",
"| 2 | rho_s | 0.5 | 0.3 | | |0.0253049| 2.0747 | |\n",
"| 3 | p4415_p | -2.5 | 0.3 | | |-6.28319 | 6.28319 | |\n",
"| 4 | psi2s_p | 2.08 | 0.05 | | |-6.28319 | 6.28319 | |\n",
"| 5 | Dbar_p | -0.9 | 0.5 | | |-6.28319 | 6.28319 | |\n",
"| 6 | omega_p | 0.5 | 0.5 | | |-6.28319 | 6.28319 | |\n",
"| 7 | DDstar_p | 2.1 | 1.7 | | |-6.28319 | 6.28319 | |\n",
"| 8 | jpsi_p | -1.46 | 0.08 | | |-6.28319 | 6.28319 | |\n",
"| 9 | p3770_s | 2.7 | 0.3 | | |0.918861 | 4.08114 | |\n",
"| 10| p4415_s | 1.06 | 0.19 | | |0.126447 | 2.35355 | |\n",
"| 11| bplus_1 | -0.879 | 0.022 | | | -2 | 2 | |\n",
"| 12| bplus_2 | -0.17 | 0.07 | | | -2 | 2 | |\n",
"| 13| p3770_p | -2.78 | 0.18 | | |-6.28319 | 6.28319 | |\n",
"| 14| bplus_0 | 0.452 | 0.011 | | | -2 | 2 | |\n",
"| 15| p4160_p | -2.13 | 0.26 | | |-6.28319 | 6.28319 | |\n",
"| 16| p4160_s | 2.00 | 0.17 | | | 0.71676 | 3.68324 | |\n",
"| 17| Dbar_s | -0.30 | 0.46 | | | -0.3 | 0.3 | |\n",
"| 18| omega_s | 6.5 | 2.1 | | | 4.19232 | 9.40768 | |\n",
"| 19| phi_p | 0.75 | 0.17 | | |-6.28319 | 6.28319 | |\n",
"| 20| Ctt | 0.31 | 0.20 | | | -2.5 | 2.5 | |\n",
"| 21| p4040_p | -2.67 | 0.28 | | |-6.28319 | 6.28319 | |\n",
"| 22| DDstar_s | -0.15 | 0.34 | | | -0.3 | 0.3 | |\n",
"| 23| p4040_s | 1.31 | 0.17 | | |0.00501244| 2.01499 | |\n",
"----------------------------------------------------------------------------------------------\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"| | phi_s rho_p rho_s p4415_p psi2s_p Dbar_p omega_p DDstar_p jpsi_p p3770_s p4415_s bplus_1 bplus_2 p3770_p bplus_0 p4160_p p4160_s Dbar_s omega_s phi_p Ctt p4040_p DDstar_s p4040_s |\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"| phi_s | 1.000 -0.153 -0.029 -0.022 -0.005 0.021 -0.014 -0.001 0.005 -0.017 0.010 -0.010 -0.129 -0.010 -0.016 -0.012 -0.026 -0.001 -0.033 0.529 -0.020 -0.017 -0.006 -0.010 |\n",
"| rho_p | -0.153 1.000 -0.006 0.062 0.025 -0.021 0.089 0.018 0.001 0.040 -0.034 0.059 0.273 0.025 -0.049 0.043 0.051 0.016 0.243 -0.306 0.055 0.048 0.039 0.015 |\n",
"| rho_s | -0.029 -0.006 1.000 -0.002 0.007 -0.047 -0.088 0.046 0.027 0.004 0.032 -0.038 -0.046 0.037 0.213 0.005 0.006 0.014 -0.314 0.093 -0.013 0.006 0.002 0.008 |\n",
"| p4415_p | -0.022 0.062 -0.002 1.000 0.616 0.342 0.026 0.737 0.757 0.549 -0.278 0.035 0.084 0.612 0.072 0.806 0.085 0.541 0.024 0.023 -0.144 0.746 0.818 -0.185 |\n",
"| psi2s_p | -0.005 0.025 0.007 0.616 1.000 0.258 0.017 0.695 0.734 0.396 -0.154 0.073 -0.040 0.577 0.080 0.698 0.092 0.453 0.011 0.032 -0.111 0.642 0.778 -0.055 |\n",
"| Dbar_p | 0.021 -0.021 -0.047 0.342 0.258 1.000 0.014 0.492 0.205 0.318 -0.159 -0.142 -0.017 -0.003 -0.094 0.336 0.134 0.393 0.019 0.015 0.218 0.341 0.395 -0.027 |\n",
"| omega_p | -0.014 0.089 -0.088 0.026 0.017 0.014 1.000 0.011 0.006 0.018 -0.017 0.023 0.051 0.009 -0.052 0.023 0.009 0.011 0.884 -0.136 0.003 0.023 0.026 -0.002 |\n",
"| DDstar_p | -0.001 0.018 0.046 0.737 0.695 0.492 0.011 1.000 0.856 0.628 -0.220 0.206 -0.070 0.693 0.179 0.858 0.123 0.679 -0.004 0.048 -0.153 0.814 0.909 -0.139 |\n",
"| jpsi_p | 0.005 0.001 0.027 0.757 0.734 0.205 0.006 0.856 1.000 0.625 -0.211 0.110 0.007 0.749 0.107 0.865 0.129 0.577 -0.007 0.051 -0.188 0.819 0.923 -0.115 |\n",
"| p3770_s | -0.017 0.040 0.004 0.549 0.396 0.318 0.018 0.628 0.625 1.000 -0.138 0.070 -0.055 0.357 0.090 0.610 0.131 0.394 0.015 0.021 -0.206 0.584 0.670 0.013 |\n",
"| p4415_s | 0.010 -0.034 0.032 -0.278 -0.154 -0.159 -0.017 -0.220 -0.211 -0.138 1.000 0.116 -0.140 -0.196 0.066 -0.239 0.277 -0.151 -0.024 0.002 0.181 -0.178 -0.225 0.191 |\n",
"| bplus_1 | -0.010 0.059 -0.038 0.035 0.073 -0.142 0.023 0.206 0.110 0.070 0.116 1.000 -0.271 0.170 -0.511 0.073 0.022 0.078 0.035 -0.018 -0.160 0.078 0.067 0.013 |\n",
"| bplus_2 | -0.129 0.273 -0.046 0.084 -0.040 -0.017 0.051 -0.070 0.007 -0.055 -0.140 -0.271 1.000 -0.069 -0.167 -0.007 0.110 -0.010 0.078 -0.085 0.594 -0.010 -0.000 0.066 |\n",
"| p3770_p | -0.010 0.025 0.037 0.612 0.577 -0.003 0.009 0.693 0.749 0.357 -0.196 0.170 -0.069 1.000 0.166 0.740 0.051 0.392 -0.002 0.037 -0.363 0.703 0.724 -0.181 |\n",
"| bplus_0 | -0.016 -0.049 0.213 0.072 0.080 -0.094 -0.052 0.179 0.107 0.090 0.066 -0.511 -0.167 0.166 1.000 0.099 0.037 0.081 -0.096 0.034 -0.124 0.103 0.095 0.007 |\n",
"| p4160_p | -0.012 0.043 0.005 0.806 0.698 0.336 0.023 0.858 0.865 0.610 -0.239 0.073 -0.007 0.740 0.099 1.000 0.076 0.594 0.018 0.035 -0.314 0.825 0.913 -0.302 |\n",
"| p4160_s | -0.026 0.051 0.006 0.085 0.092 0.134 0.009 0.123 0.129 0.131 0.277 0.022 0.110 0.051 0.037 0.076 1.000 0.122 0.010 -0.007 0.235 0.283 0.175 -0.018 |\n",
"| Dbar_s | -0.001 0.016 0.014 0.541 0.453 0.393 0.011 0.679 0.577 0.394 -0.151 0.078 -0.010 0.392 0.081 0.594 0.122 1.000 0.005 0.030 -0.195 0.553 0.732 -0.046 |\n",
"| omega_s | -0.033 0.243 -0.314 0.024 0.011 0.019 0.884 -0.004 -0.007 0.015 -0.024 0.035 0.078 -0.002 -0.096 0.018 0.010 0.005 1.000 -0.201 0.010 0.018 0.021 -0.002 |\n",
"| phi_p | 0.529 -0.306 0.093 0.023 0.032 0.015 -0.136 0.048 0.051 0.021 0.002 -0.018 -0.085 0.037 0.034 0.035 -0.007 0.030 -0.201 1.000 -0.021 0.031 0.041 -0.009 |\n",
"| Ctt | -0.020 0.055 -0.013 -0.144 -0.111 0.218 0.003 -0.153 -0.188 -0.206 0.181 -0.160 0.594 -0.363 -0.124 -0.314 0.235 -0.195 0.010 -0.021 1.000 -0.263 -0.217 0.342 |\n",
"| p4040_p | -0.017 0.048 0.006 0.746 0.642 0.341 0.023 0.814 0.819 0.584 -0.178 0.078 -0.010 0.703 0.103 0.825 0.283 0.553 0.018 0.031 -0.263 1.000 0.877 -0.210 |\n",
"| DDstar_s | -0.006 0.039 0.002 0.818 0.778 0.395 0.026 0.909 0.923 0.670 -0.225 0.067 -0.000 0.724 0.095 0.913 0.175 0.732 0.021 0.041 -0.217 0.877 1.000 -0.105 |\n",
"| p4040_s | -0.010 0.015 0.008 -0.185 -0.055 -0.027 -0.002 -0.139 -0.115 0.013 0.191 0.013 0.066 -0.181 0.007 -0.302 -0.018 -0.046 -0.002 -0.009 0.342 -0.210 -0.105 1.000 |\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"Hesse errors: OrderedDict([(<zfit.Parameter 'phi_s' floating=True>, {'error': 0.9142957081972742}), (<zfit.Parameter 'rho_p' floating=True>, {'error': 0.6936535176060628}), (<zfit.Parameter 'rho_s' floating=True>, {'error': 0.31992987771310377}), (<zfit.Parameter 'p4415_p' floating=True>, {'error': 0.3242479743998188}), (<zfit.Parameter 'psi2s_p' floating=True>, {'error': 0.04877084961929512}), (<zfit.Parameter 'Dbar_p' floating=True>, {'error': 0.4693624571295927}), (<zfit.Parameter 'omega_p' floating=True>, {'error': 0.5448127134133323}), (<zfit.Parameter 'DDstar_p' floating=True>, {'error': 1.6782936194122264}), (<zfit.Parameter 'jpsi_p' floating=True>, {'error': 0.076502815887741}), (<zfit.Parameter 'p3770_s' floating=True>, {'error': 0.3269991972653048}), (<zfit.Parameter 'p4415_s' floating=True>, {'error': 0.19331473430113788}), (<zfit.Parameter 'bplus_1' floating=True>, {'error': 0.022453161599432336}), (<zfit.Parameter 'bplus_2' floating=True>, {'error': 0.07164712471903711}), (<zfit.Parameter 'p3770_p' floating=True>, {'error': 0.18127978100197017}), (<zfit.Parameter 'bplus_0' floating=True>, {'error': 0.010823487498148543}), (<zfit.Parameter 'p4160_p' floating=True>, {'error': 0.26287545614952745}), (<zfit.Parameter 'p4160_s' floating=True>, {'error': 0.16605504500407542}), (<zfit.Parameter 'Dbar_s' floating=True>, {'error': 0.4621450517273268}), (<zfit.Parameter 'omega_s' floating=True>, {'error': 2.074388629496891}), (<zfit.Parameter 'phi_p' floating=True>, {'error': 0.16789257624839893}), (<zfit.Parameter 'Ctt' floating=True>, {'error': 0.2011430528044098}), (<zfit.Parameter 'p4040_p' floating=True>, {'error': 0.28336427935455255}), (<zfit.Parameter 'DDstar_s' floating=True>, {'error': 0.3358643253844139}), (<zfit.Parameter 'p4040_s' floating=True>, {'error': 0.1701201938999789})])\n",
"\n",
"Time taken: 3 min, 37 s\n",
"Estimated time left: 18 min, 5 s\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEICAYAAAC9E5gJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAYU0lEQVR4nO3dfZBldX3n8fcngBAFwwANxQLaYCEbTOnoThFTRlfFB5CUYGKyQyxlE3ZnXaXWp6rNGHaN2SpThI0aLY3UGCmxFEQFVirgwxSLca31IQ3CMOygDASXkclMi4oYXBLwu3/c03Jpbnff7vvQ9+H9qrp1z/2d3zn3d37d9/c559yHk6pCkjTdfmm9GyBJWn+GgSTJMJAkGQaSJAwDSRKGgSQJw0AiybuTfHK92yGtJ8NA6yLJ7yeZS/LTJHuTfCHJbzbznjA4J/lKkn+3Pq3tnySnJ7kjyUNJbkzy9GXqzjZ1HmqWedmi+W9L8g9JHkhyaZKDB78FmlSGgYYuyduBvwT+DDgGeBrwV8DZ69muQUtyFHA18F+BI4A54MplFrkC+DZwJHAh8LkkM826XglsBU4HZoGTgD8dVNs1BarKm7eh3YBfAX4K/O4S888A/gn456bercB7gEeB/9eUfajDcl8ELlhUdivw2830B4B7gZ8ANwEvbKv3buCTzfSLgT2L1nMP8LJm+pdoDcJ3AfcDnwGO6HLbtwD/u+3xU4CfAf+yQ91nAg8Dh7WV/S/gjc305cCftc07HfiH9f77ehvfm0cGGrbfAA4Bruk0s6q+SOuI4cqqOrSqnlNVF9IaCC9oyi7osOjlwLkLD5KcCjwduK4p+jtgI6098suBzyY5ZA3t/0/AOcC/Bv4F8CPgw23PuyPJ7y+x7LNoBdTCtv4jrVB51hJ1766qB9vKbm2r+7h1NdPHJDlyVVsjNQwDDduRwA+q6pE+r/caYGPbOfjXAVdX1cMAVfXJqrq/qh6pqvcCBwOnrOF5/gNwYVXtadb9buC1SQ5snufZVXX5EsseCjywqOwB4LA11F08f2G607qkFRkGGrb7gaMWBs9+afagrwM2N0WbgU8tzE/yjiS7mjdbf0zrdNVRa3iqpwPXJPlxs55dtE5hHdPFsj8Fnrqo7KnAg2uou3j+wnSndUkrMgw0bF+nde7/nGXqdPop3W5+XvcK4NwkvwH8MnAjQJIXAn8E/B6woaoOp7UnnQ7r+EfgyQsPkhwAzLTNvxc4s6oOb7sdUlXf76J9twPPaVv3U4BnNOWd6p6UpH1P/zltdR+3rmZ6X1Xd30U7pCcwDDRUVfUA8C7gw0nOSfLkJAclOTPJxU21fcBskvb/z320PjGznOtp7bn/N1rvOfy8KT8MeASYBw5M8i6euNe94LvAIUnOSnIQ8F9onVJacAnwnoXTUUlmknT7KahrgF9L8jvN+xXvAnZU1R2LK1bVd4FbgD9JckiS1wDPBq5qqnwCOD/JqUk2NO38eJftkJ7AMNDQVdX7gLfTGsDmae1tXwD8j6bKZ5v7+5Pc3Ex/gNa5+R8l+eAS632Y1kc3X0brTeIFXwK+QGug/x6tI5N7l1jHA8CbgL8Gvk/rSGFPW5UPANcCX07yIPAN4NcXZia5Pcnrllj3PPA7tD4d9aNmuc1ty16S5JK2RTYDm5q6FwGvbdax8Eb7xbSOfr7X3P6k0/NK3UiVF7eRpGnnkYEkyTCQJBkGkiQMA0kS0Ncv/qzVUUcdVbOzs+vdDEkaKzfddNMPqmpm5ZorG4kwmJ2dZW5ubr2bIUljJcn3+rUuTxNJkgwDSZJhIEnCMJAkYRhIkjAMJEkYBpIkDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJAwDaazMbr1uvZugCbViGCQ5IcmNSXYluT3JW5ryI5JsT3Jnc7+hKU+SDybZnWRHkucNeiMkSb3p5sjgEeAdVfWrwPOBNyc5FdgK3FBVJwM3NI8BzgRObm5bgI/0vdWSpL5aMQyqam9V3dxMPwjsAo4DzgYua6pdBpzTTJ8NfKJavgEcnuTYvrdcktQ3q3rPIMks8Fzgm8AxVbUXWoEBHN1UOw64t22xPU2ZJGlEdR0GSQ4FrgLeWlU/Wa5qh7LqsL4tSeaSzM3Pz3fbDEkTwDfCR09XYZDkIFpB8Kmqurop3rdw+qe539+U7wFOaFv8eOC+xeusqm1VtamqNs3MzKy1/ZKkPujm00QBPgbsqqr3tc26FjivmT4P+Hxb+RuaTxU9H3hg4XSSJGk0HdhFnRcArwduS3JLU/bHwEXAZ5KcD/xf4HebedcDrwJ2Aw8Bf9DXFkuS+m7FMKiqr9H5fQCA0zvUL+DNPbZLkjREfgNZkmQYSJIMA0kShoEkCcNAkoRhIEnCMJAkYRhImgL+FtLKDANJkmEgSTIMJEkYBpIkDANJEoaBJAnDQJKEYSBJorvLXl6aZH+SnW1lVya5pbnds3AFtCSzSX7WNu+SQTZektQf3Vz28uPAh4BPLBRU1b9ZmE7yXuCBtvp3VdXGfjVQkjR43Vz28qtJZjvNSxLg94CX9rdZkqRh6vU9gxcC+6rqzrayE5N8O8nfJnnhUgsm2ZJkLsnc/Px8j82QJPWi1zA4F7ii7fFe4GlV9Vzg7cDlSZ7aacGq2lZVm6pq08zMTI/NkCT1Ys1hkORA4LeBKxfKqurhqrq/mb4JuAt4Zq+NlCQNVi9HBi8D7qiqPQsFSWaSHNBMnwScDNzdWxMlSYPWzUdLrwC+DpySZE+S85tZm3n8KSKAFwE7ktwKfA54Y1X9sJ8NliT1XzefJjp3ifJ/26HsKuCq3pslSRomv4EsSTIMJEmGgSQJw0CShGEgScIwkCRhGEgastmt1613E9SBYSBJMgwkSYaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJdHdxm0uT7E+ys63s3Um+n+SW5vaqtnnvTLI7yXeSvHJQDZck9U83RwYfB87oUP7+qtrY3K4HSHIqrSugPatZ5q8WLoMpSRpdK4ZBVX0V6PbSlWcDn66qh6vq74HdwGk9tE+SNAS9vGdwQZIdzWmkDU3ZccC9bXX2NGWSeuRv+miQ1hoGHwGeAWwE9gLvbcrToW51WkGSLUnmkszNz8+vsRmStLyFEDVMl7emMKiqfVX1aFX9HPgoj50K2gOc0Fb1eOC+Jdaxrao2VdWmmZmZtTRDkpZlAHRvTWGQ5Ni2h68BFj5pdC2wOcnBSU4ETga+1VsTJak/DIeldfPR0iuArwOnJNmT5Hzg4iS3JdkBvAR4G0BV3Q58Bvg/wBeBN1fVowNrvaSOpmHQm4ZtHKYDV6pQVed2KP7YMvXfA7ynl0ZJmnyzW6/jnovOWu9mqOE3kCVJhoEkyTCQJGEYSJIwDCRJGAaSJAwDSRKGgaQx55fP+sMwkLSuHMxHg2EgSTIMJEmGgaQp42mpzgwDSZJhIE2yUdsLHrX26DGGgaSRZ4gMnmEgSerqSmeXJtmfZGdb2X9PckeSHUmuSXJ4Uz6b5GdJbmlulwyy8ZLGx7D37j2aWJ1ujgw+DpyxqGw78GtV9Wzgu8A72+bdVVUbm9sb+9NMSVqeg39vVgyDqvoq8MNFZV+uqkeah98Ajh9A2yRJQ9KP9wz+EPhC2+MTk3w7yd8meeFSCyXZkmQuydz8/HwfmiFp3Lg3Pzp6CoMkFwKPAJ9qivYCT6uq5wJvBy5P8tROy1bVtqraVFWbZmZmemmGpCllmPTPmsMgyXnAbwGvq6oCqKqHq+r+Zvom4C7gmf1oqCQtxVDo3ZrCIMkZwB8Br66qh9rKZ5Ic0EyfBJwM3N2PhkrSAgf//uvmo6VXAF8HTkmyJ8n5wIeAw4Dtiz5C+iJgR5Jbgc8Bb6yqH3ZcsSR10MtAv9ZlDRc4cKUKVXVuh+KPLVH3KuCqXhslSRouv4EsTahx2tsdp7ZOKsNAkpYxLUFlGEiaOtMywK+GYSBJMgwkSYaBpBE0zNM4yz3XNJ1OMgykMdCPQWkcBrZxaOOkMgykEecAqWEwDCStyaiG1CDaNarb2k+GgTTCpmEQGieT/PcwDCRJhoE06SZ5b1b9YxhIEoamYSBJMgwkjYZu9synfe99kLoKgySXJtmfZGdb2RFJtie5s7nf0JQnyQeT7E6yI8nzBtV4Sb0Z5cF1qbaNcpvHWbdHBh8HzlhUthW4oapOBm5oHgOcSetylycDW4CP9N5MSdPIgX94VrzSGUBVfTXJ7KLis4EXN9OXAV+hdV3ks4FPVFUB30hyeJJjq2pvPxosSSvxlNPq9fKewTELA3xzf3RTfhxwb1u9PU2ZJGlEDeIN5HQoqydUSrYkmUsyNz8/P4BmSOqFe87TpZcw2JfkWIDmfn9Tvgc4oa3e8cB9ixeuqm1VtamqNs3MzPTQDEnDMg0/97zSdkzKdi7WSxhcC5zXTJ8HfL6t/A3Np4qeDzzg+wXSaFjtQDapA5+eqNuPll4BfB04JcmeJOcDFwEvT3In8PLmMcD1wN3AbuCjwJv63mppgvmrm6Nh2vqs208TnbvErNM71C3gzb00StLwzG69jnsuOmvdlu9m/cMybQHQzm8gS1Nq3Aa+cWvvuDEMJHXFwfgxk9gXhoEkYDIHOHXPMJD0C90GgsExeQwDSatiEEwmw0CSZBhIo2qte+BrWc69fRkGktSFSQ9Mw0CaMtP62ztanmEgDYEDrEadYSBpWQbZdDAMpHU0LgPtuLRzmCatTwwDSZJhIK2H9d6r7Nfzr/d2jIJJ6QPDQBpB4zbAjFt7+2WSttswkCR1d3GbTpKcAlzZVnQS8C7gcODfAwtXuf/jqrp+zS2UJA3cmo8Mquo7VbWxqjYC/wp4CLimmf3+hXkGgabZJJ1GmCbT+Hfr12mi04G7qup7fVqfNJGmcZDReOhXGGwGrmh7fEGSHUkuTbKh0wJJtiSZSzI3Pz/fqYo0NQwJrbeewyDJk4BXA59tij4CPAPYCOwF3ttpuaraVlWbqmrTzMxMr82QRoYDu8ZRP44MzgRurqp9AFW1r6oeraqfAx8FTuvDc0hjz5DQKOtHGJxL2ymiJMe2zXsNsLMPzyFJI2sSgr6nMEjyZODlwNVtxRcnuS3JDuAlwNt6eQ5pmLwGsKbVmr9nAFBVDwFHLip7fU8tkiacQaJR5DeQpRExSiExSm3RcBgG0jpxwNUoMQykPmof4Nfy/oMBofViGEiSDANJkmEgdW2pUzidyj3do3FjGEir4CCvSWUYSFIfjesOg2EgSTIMpNVa2PMb1z1AqRPDQGPJgVjqL8NA6hMDSuPMMJAWWc2g7ikjTQrDQFNjvQbsaQqKadrWSWMYaCo4SEnLMwykJRggmiY9h0GSe5orm92SZK4pOyLJ9iR3Nvcbem+qJsmon7JZXM9g0ErG/X+kX0cGL6mqjVW1qXm8Fbihqk4Gbmgea4oM84Ux7i9CTY5x/l8c1Gmis4HLmunLgHMG9DxaJ+P8T78as1uvm5ptVX+N2/9NP8KggC8nuSnJlqbsmKraC9DcH714oSRbkswlmZufn+9DMzSpBnXKZtxerBo/4/Q/dmAf1vGCqrovydHA9iR3dLNQVW0DtgFs2rSp+tAOac3G6UWr0TeO/089HxlU1X3N/X7gGuA0YF+SYwGa+/29Po+0knF8AUqjoqcwSPKUJIctTAOvAHYC1wLnNdXOAz7fy/NoejigS+uj1yODY4CvJbkV+BZwXVV9EbgIeHmSO4GXN4/VJQfE7thPUv/09J5BVd0NPKdD+f3A6b2sW+qH2a3Xcc9FZ3Us7zQtTSu/gayR5kAtDYdhoJHRzcC/Uh33+KW1MQy0bpb7QtdqB/Ll6hsK0soMA/VVPwZeB29p+AwD9X3wXe3efr+ODtaynMEjtRgG6smo/nbPKLZJGmWGgQauH78t1M/B3aCQnsgw0JJ6vb7vMD7Z48Au9YdhMKUcRCW1MwzUkWEhTRfDYMSMwyA8Dm2UtDqGgZ6g07n+fnw7eFDLSuqdYTAlBjHYOoBLk8MwmECDHKT7+TMR0rQZ5deDYaCxMcovJGkp4/J/u+YwSHJCkhuT7Epye5K3NOXvTvL9JLc0t1f1r7mTq9d/mLWc0x+Xf1JJg9fLxW0eAd5RVTc3l768Kcn2Zt77q+ovem+e+mGUBv1Raos0TKP+v7/mI4Oq2ltVNzfTDwK7gOP61bBxNsw/+nLf8h3WD9BJWt44vHb68p5BklngucA3m6ILkuxIcmmSDUsssyXJXJK5+fn5fjRjbKxXWKw0fxz+YSUNRs9hkORQ4CrgrVX1E+AjwDOAjcBe4L2dlquqbVW1qao2zczM9NqMiTUue/cGiTTeegqDJAfRCoJPVdXVAFW1r6oeraqfAx8FTuu9mdNtGL/P3+uP0kkab718mijAx4BdVfW+tvJj26q9Bti59uZNtkF+gsjLQEpajV6ODF4AvB546aKPkV6c5LYkO4CXAG/rR0NH0bC/gOUgLmlQ1vzR0qr6GpAOs65fe3Mmx+zW67jnorMe97ibesutT5IGxW8gj4BOA/1y4eGXxyT1m2GwCqsZtCVpnBgGKxjVwX5U2yVpPBkGazSIL2s5wEtaL4ZBox8D8TC+DyBJg2AYLMP3CCRNC8NgDVbzhS4vBSlpHExVGAxyYPYbv5LG2VSFQSfLfZ5/pTprXbek6TWq48LUhYE/2SxJTzR1YQCGgKT1NYpjUC+XvRwra7lGsCQNSre/SzYsE31k0O3gPkpXHpOk9TDRYQDdvUEsSdNuosLAgV+S1maiwmCBg78krc7AwiDJGUm+k2R3kq2Deh5JUu8GEgZJDgA+DJwJnAqcm+TUQTzXUjw6kKTuDerI4DRgd1XdXVX/BHwaOHtAz/U4hoAkrd6gvmdwHHBv2+M9wK+3V0iyBdjSPPxpkvuBHwyoPePmKOyLBfZFyy/6IX++zi1ZwRDaNzH/Ez321VHA0/vTksGFQTqU1eMeVG0Dtv1igWSuqjYNqD1jxb54jH3RYj88xr5oafphtl/rG9Rpoj3ACW2PjwfuG9BzSZJ6NKgw+Dvg5CQnJnkSsBm4dkDPJUnq0UBOE1XVI0kuAL4EHABcWlW3r7DYthXmTxP74jH2RYv98Bj7oqWv/ZCqWrmWJGmiTeQ3kCVJq2MYSJJGIwym4acrklyaZH+SnW1lRyTZnuTO5n5DU54kH2z6Y0eS57Utc15T/84k563HtvQiyQlJbkyyK8ntSd7SlE9VXyQ5JMm3ktza9MOfNuUnJvlms01XNh/AIMnBzePdzfzZtnW9syn/TpJXrs8W9S7JAUm+neRvmsdT2RdJ7klyW5Jbksw1ZYN/fVTVut5ovcF8F3AS8CTgVuDU9W7XALbzRcDzgJ1tZRcDW5vprcCfN9OvAr5A6/sazwe+2ZQfAdzd3G9opjes97atsh+OBZ7XTB8GfJfWT5ZMVV8023NoM30Q8M1m+z4DbG7KLwH+YzP9JuCSZnozcGUzfWrzmjkYOLF5LR2w3tu3xj55O3A58DfN46nsC+Ae4KhFZQN/fYzCkcG6/XTFMFXVV4EfLio+G7ismb4MOKet/BPV8g3g8CTHAq8EtlfVD6vqR8B24IzBt75/qmpvVd3cTD8I7KL1jfWp6otme37aPDyouRXwUuBzTfnifljon88BpydJU/7pqnq4qv4e2E3rNTVWkhwPnAX8dfM4TGlfLGHgr49RCINOP11x3Dq1ZdiOqaq90BokgaOb8qX6ZKL6qjm8fy6tveKp64vmtMgtwH5aL9a7gB9X1SNNlfZt+sX2NvMfAI5kAvqh8ZfAfwZ+3jw+kuntiwK+nOSmtH62B4bw+hiFayCv+NMVU2ipPpmYvkpyKHAV8Naq+klrx65z1Q5lE9EXVfUosDHJ4cA1wK92qtbcT2w/JPktYH9V3ZTkxQvFHapOfF80XlBV9yU5Gtie5I5l6vatL0bhyGCaf7piX3NIR3O/vylfqk8moq+SHEQrCD5VVVc3xVPZFwBV9WPgK7TO+R6eZGEnrX2bfrG9zfxfoXXacRL64QXAq5PcQ+s08UtpHSlMY19QVfc19/tp7SScxhBeH6MQBtP80xXXAgvv8p8HfL6t/A3NJwWeDzzQHBp+CXhFkg3Npwle0ZSNjebc7seAXVX1vrZZU9UXSWaaIwKS/DLwMlrvn9wIvLaptrgfFvrntcD/rNY7hdcCm5tP2JwInAx8azhb0R9V9c6qOr5aP7q2mda2vY4p7IskT0ly2MI0rf/rnQzj9bHe75y3vSP+XVrnTC9c7/YMaBuvAPYC/0wrtc+ndZ7zBuDO5v6Ipm5oXRzoLuA2YFPbev6Q1htju4E/WO/tWkM//Catw9UdwC3N7VXT1hfAs4FvN/2wE3hXU34SrQFsN/BZ4OCm/JDm8e5m/klt67qw6Z/vAGeu97b12C8v5rFPE01dXzTbfGtzu31hPBzG68Ofo5AkjcRpIknSOjMMJEmGgSTJMJAkYRhIkjAMJEkYBpIk4P8DqE+/Gh539BkAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# zfit.run.numeric_checks = False \n",
"\n",
"load = False\n",
"\n",
"bo = True\n",
"\n",
"D_contribs = True\n",
"\n",
"if not D_contribs:\n",
" Dbar_s.floating = False\n",
" Dbar_p.floating = False\n",
" DDstar_s.floating = False\n",
" DDstar_p.floating = False\n",
"\n",
"bo_set = 1\n",
"\n",
"fitting_range = 'cut'\n",
"total_BR = 1.7e-10 + 4.9e-10 + 2.5e-9 + 6.02e-5 + 4.97e-6 + 1.38e-9 + 4.2e-10 + 2.6e-9 + 6.1e-10 + 4.37e-7\n",
"cut_BR = 1.0 - (6.02e-5 + 4.97e-6)/total_BR\n",
"\n",
"Ctt_list = []\n",
"Ctt_error_list = []\n",
"\n",
"nr_of_toys = 1\n",
"nevents = int(pdg[\"number_of_decays\"]*cut_BR)\n",
"# nevents = pdg[\"number_of_decays\"]\n",
"event_stack = 1000000\n",
"# nevents *= 41\n",
"# zfit.settings.set_verbosity(10)\n",
"\n",
"# mi = 1e-4\n",
"# ma = 3e-3\n",
"# ste = 13\n",
"\n",
"BR_steps = np.linspace(mi, ma, ste)\n",
"\n",
"BR_steps[0] = 0.0\n",
"\n",
"print(BR_steps)\n",
"\n",
"Ctt_steps = np.sqrt(BR_steps/4.2*1000)\n",
"\n",
"print(Ctt_steps)\n",
"\n",
"# total_samp = []\n",
"\n",
"start = time.time()\n",
"\n",
"Nll_list = []\n",
"\n",
"sampler = total_f.create_sampler(n=nevents, fixed_params = False)\n",
"sampler.set_data_range(obs_fit)\n",
"\n",
"__ = -1\n",
"\n",
"newset = True\n",
"\n",
"#-----------------------------------------------------\n",
"\n",
"if not load:\n",
" for Ctt_step in Ctt_steps:\n",
"\n",
" if not newset:\n",
" break\n",
" \n",
" __ += 1\n",
" \n",
" for i in range(2):\n",
" Ctt_list.append([])\n",
" Ctt_error_list.append([])\n",
" Nll_list.append([])\n",
"\n",
" for param in total_f_fit.get_dependents():\n",
" if param.floating:\n",
" pull_dic[param.name].append([])\n",
" \n",
" for toy in range(nr_of_toys): \n",
" \n",
" if not newset:\n",
" break\n",
" \n",
" newset = True\n",
" \n",
" while newset:\n",
" \n",
" for floaty in [True, False]:\n",
" Ctt.floating = floaty\n",
" \n",
" if not floaty:\n",
" break\n",
" \n",
" for bo_step in range(bo_set):\n",
"\n",
" print('Step: {0}/{1}'.format(int(__), ste))\n",
" print('Current Ctt: {0}'.format(Ctt_step))\n",
" print('Ctt floating: {0}'.format(floaty))\n",
" \n",
" reset_param_values(variation = 0.0)\n",
"\n",
" if floaty:\n",
" print('Toy {0}/{1} - Fit {2}/{3}'.format(toy, nr_of_toys, bo_step, bo_set))\n",
" Ctt.set_value(Ctt_step)\n",
"\n",
" else:\n",
" Ctt.set_value(0.0)\n",
" print('Toy {0}/{1} - Fit {2}/{3}'.format(toy, nr_of_toys, bo_step, bo_set))\n",
"\n",
" if newset:\n",
" sampler.resample(n=nevents)\n",
" data = sampler\n",
" newset = False\n",
"\n",
" ### Fit data\n",
" \n",
" if floaty:\n",
" plt.clf()\n",
" plt.title('Ctt value: {:.2f}'.format(Ctt_step))\n",
" plt.hist(zfit.run(data), bins = int((x_max-x_min)/7), range = (x_min, x_max))\n",
" plt.savefig('data/CLs/plots/set_histo{}.png'.format(__))\n",
" _step = 2*__\n",
" else:\n",
" _step = 2*__+1\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",
" \n",
" save_pulls(step = _step)\n",
"\n",
" if floaty:\n",
" Nll_list[-2].append(result.fmin)\n",
" Ctt_list[-2].append(params[Ctt]['value'])\n",
" Ctt_error_list[-2].append(params[Ctt]['minuit_hesse']['error'])\n",
"\n",
" else:\n",
" Nll_list[-1].append(result.fmin)\n",
" Ctt_list[-1].append(0.0)\n",
" Ctt_error_list[-1].append(0.0)\n",
" \n",
"\n",
" else:\n",
" for _ in [1,2]:\n",
" del Nll_list[-_][toy*bo_set:]\n",
"# print(np.shape(Nll_list[-_]))\n",
" del Ctt_list[-_][toy*bo_set:]\n",
" del Ctt_error_list[-_][toy*bo_set:]\n",
" for param in total_f_fit.get_dependents():\n",
" if param.floating:\n",
" del pull_dic[param.name][_step+1-_][toy*bo_set:]\n",
" newset = True\n",
" break\n",
" \n",
" if not result.converged:\n",
" break\n",
" \n",
" print()\n",
" print('Time taken: {}'.format(display_time(int(time.time()-start))))\n",
" print('Estimated time left: {}'.format(display_time(int((time.time()-start)/(__+(toy+1)/nr_of_toys)*(ste-__-(nr_of_toys-toy-1)/nr_of_toys)))))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"if load:\n",
" \n",
" phase_combi = '-+'\n",
" \n",
" if D_contribs:\n",
" D_dir = 'D-True'\n",
" else:\n",
" D_dir = 'D-False'\n",
"\n",
" _dir = 'data/CLs/finished/f3d41/{}/{}/'.format(phase_combi, D_dir)\n",
" \n",
" jobs = os.listdir(_dir)\n",
" \n",
" First = True\n",
" \n",
" print('Number of jobs: {}'.format(len(jobs)))\n",
" \n",
" for job in jobs:\n",
" \n",
" dirName = _dir + str(job) + '/data/CLs'\n",
" \n",
" if not os.path.exists(\"{}/{}-{}_{}s{}b{}t--CLs_Nll_list.pkl\".format(dirName, mi,ma,ste,bo_set,nr_of_toys)):\n",
" print(job)\n",
" continue\n",
" \n",
" with open(r\"{}/variab.pkl\".format(dirName), \"rb\") as input_file:\n",
" variab = pkl.load(input_file)\n",
"# print(variab)\n",
" \n",
" ### sanity check:\n",
" if variab['mi'] != mi or variab['ma'] != ma or variab['ste'] != ste or bo_set != bo_set:\n",
" print('Fitting parameters of data dont equal the ones given -- Job {} skipped!'.format(job))\n",
" continue\n",
" \n",
" with open(r\"{}/{}-{}_{}s{}b{}t--CLs_Nll_list.pkl\".format(dirName, mi,ma,ste,bo_set,nr_of_toys), \"rb\") as input_file:\n",
" _Nll_list = pkl.load(input_file)\n",
" \n",
" with open(r\"{}/{}-{}_{}s{}b{}t--Ctt_list.pkl\".format(dirName, mi,ma,ste,bo_set,nr_of_toys), \"rb\") as input_file:\n",
" _Ctt_list = pkl.load(input_file)\n",
" \n",
" with open(r\"{}/{}-{}_{}s{}b{}t--Ctt_error_list.pkl\".format(dirName, mi,ma,ste,bo_set,nr_of_toys), \"rb\") as input_file:\n",
" _Ctt_error_list = pkl.load(input_file)\n",
" \n",
" with open(r\"{}/{}-{}_{}s{}b{}t--pull_dic.pkl\".format(dirName, mi,ma,ste,bo_set,nr_of_toys), \"rb\") as input_file:\n",
" _pull_dic = pkl.load(input_file)\n",
" \n",
" with open(r\"{}/{}-{}_{}s--CLs_list.pkl\".format(dirName, mi,ma,ste), \"rb\") as input_file:\n",
" _CLs_list = pkl.load(input_file)\n",
" \n",
" \n",
" if First:\n",
" Nll_list = _Nll_list\n",
" Ctt_list = _Ctt_list\n",
" Ctt_error_list = _Ctt_error_list\n",
" pull_dic = _pull_dic\n",
"# print(_pull_dic)\n",
" CLs_list = _CLs_list\n",
" First = False\n",
" else:\n",
" for step in range(2*ste):\n",
"# print(Nll_list[step], step)\n",
" Nll_list[step].extend(_Nll_list[step])\n",
" Ctt_list[step].extend(_Ctt_list[step])\n",
" Ctt_error_list[step].extend(_Ctt_error_list[step])\n",
" for key in pull_dic.keys():\n",
"# print(key, np.shape(pull_dic[key]))\n",
" pull_dic[key][step].extend(_pull_dic[key][step])\n",
" for step in range(ste):\n",
" CLs_list[step].extend(_CLs_list[step])\n",
"\n",
"# print('----------------------')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dirName = 'data/CLs'\n",
"\n",
"# if bo and not load:\n",
"# for s in range(2*ste):\n",
"# Nll_list[s] = [np.min(Nll_list[s])]\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{}b{}t--CLs_Nll_list.pkl\".format(dirName, mi,ma,ste,bo_set,nr_of_toys), \"wb\") as f:\n",
" pkl.dump(Nll_list, f, pkl.HIGHEST_PROTOCOL)\n",
" \n",
" with open(\"{}/{}-{}_{}s{}b{}t--Ctt_list.pkl\".format(dirName, mi,ma,ste,bo_set,nr_of_toys), \"wb\") as f:\n",
" pkl.dump(Ctt_list, f, pkl.HIGHEST_PROTOCOL)\n",
" \n",
" with open(\"{}/{}-{}_{}s{}b{}t--Ctt_error_list.pkl\".format(dirName, mi,ma,ste,bo_set,nr_of_toys), \"wb\") as f:\n",
" pkl.dump(Ctt_error_list, f, pkl.HIGHEST_PROTOCOL)\n",
" \n",
" with open(\"{}/{}-{}_{}s{}b{}t--pull_dic.pkl\".format(dirName, mi,ma,ste,bo_set,nr_of_toys), \"wb\") as f:\n",
" pkl.dump(pull_dic, f, pkl.HIGHEST_PROTOCOL)\n",
" \n",
" variab = {'mi': mi,\n",
" 'ma': ma,\n",
" 'ste': ste,\n",
" 'bo_set': bo_set,\n",
" 'nr_of_toys': nr_of_toys}\n",
" \n",
" with open(\"{}/variab.pkl\".format(dirName), \"wb\") as f:\n",
" pkl.dump(variab, f, pkl.HIGHEST_PROTOCOL)\n",
" \n",
" CLs_values = []\n",
" \n",
" toy_size = bo_set\n",
"\n",
" print(np.shape(Nll_list))\n",
" print(Nll_list[0:1])\n",
" \n",
" for step in range(ste):\n",
" CLs_values.append([])\n",
" for toy in range(nr_of_toys):\n",
" float_min = np.min(Nll_list[2*step][toy*bo_set:(toy+1)*bo_set])\n",
" fix_min = np.min(Nll_list[2*step+1][toy*bo_set:(toy+1)*bo_set])\n",
" CLs_values[step].append(float_min-fix_min)\n",
" \n",
" \n",
" print(np.shape(CLs_values))\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(variab['mi'] != mi)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Plot"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"l = []\n",
"\n",
"if load:\n",
" CLs_values = -1*np.array(CLs_list)\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",
"print(np.shape(CLs_values))\n",
"\n",
"figure(num=None, figsize=(2.7, 5), dpi=80, facecolor='w', edgecolor='k')\n",
"\n",
"for step in range(1,ste):\n",
" plt.clf()\n",
" plt.title('Branching ratio: {:.1E}'.format(BR_steps[step]))\n",
" plt.hist(CLs_values[0], bins = 40, range = (-5, 15), label = 'Ctt = 0', alpha = 0.8)\n",
" plt.hist(CLs_values[step], bins = 40, range = (-5, 15), label = 'Ctt = {:.2f}'.format(Ctt_steps[step]), alpha = 0.7)\n",
"# plt.hist(CLs_values[0][np.where(np.array(CLs_values[0]) > np.mean(CLs_values[0]))[0]], bins = 40, range = (-5, 15), alpha = 0.9)\n",
"# plt.hist(CLs_values[step][np.where(np.array(CLs_values[step]) < np.mean(CLs_values[0]))[0]], bins = 40, range = (-5, 15), alpha = 0.9)\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[step]))\n",
" \n",
" l.append(len(np.where(np.array(CLs_values[step]) < np.mean(CLs_values[0]))[0])) \n",
"\n",
"figure() \n",
" \n",
"for step in range(2*ste):\n",
" if step%2 == 0:\n",
" floaty = True\n",
" else:\n",
" floaty = False\n",
" for key in pull_dic.keys():\n",
" if not os.path.exists('data/CLs/plots/{}'.format(key)):\n",
" os.mkdir('data/CLs/plots/{}'.format(key))\n",
" plt.clf()\n",
" plt.title('Pull {} - Ctt value {:.2f} - floating {}'.format(key, Ctt_steps[int(step/2)], floaty))\n",
" plt.hist(pull_dic[key][step], bins = 50, range = (-5,5))\n",
" plt.xlabel('Pull')\n",
" plt.savefig('data/CLs/plots/{}/{:.2f}Ctt{}s{}f.png'.format(key, Ctt_steps[int(step/2)], step, floaty))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"for s in range(len(l)):\n",
" print('BR: {:.4f}'.format(BR_steps[s+1]))\n",
" print(l[s]/len(np.where(np.array(CLs_values[0]) < np.mean(CLs_values[0]))[0]))\n",
" print()\n",
"# print(l[s], len(CLs_values[s]))\n",
"# print()\n",
"# print(len(CLs_values[0])/2)\n",
"# print(len(np.where(np.array(CLs_values[0]) < np.mean(CLs_values[0]))[0]))\n",
"# plt.clf()\n",
"# # plt.title('Pull {} - Ctt value {:.2f} - floating {}'.format(key, Ctt_steps[int(step/2)], floaty))\n",
"# plt.hist(CLs_values[0], bins = 150, range = (-25, 50), label = 'Ctt fixed to 0')\n",
"# plt.axvline(x=np.mean(CLs_values[0]),color='red', linewidth=1.0, linestyle = 'dotted')\n",
"# plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# for step in range(2*ste):\n",
"# for key in pull_dic.keys():\n",
"# print(pull_dic[key][step])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# for param in total_f_fit.get_dependents():\n",
"# if param.floating:\n",
"# print(params[param]['value'])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(display_time(int(time.time()-start)))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"probs = total_f_fit.pdf(test_q, norm_range=False)\n",
"\n",
"calcs_test = zfit.run(probs)"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAaUAAAEMCAYAAABkwamIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO29eXhc1ZWv/a6q0mjJsq3BxpZtybY8MxszT4aAyWQ6IYkh6XDTdOivL3QnnaQDdIbOTaC/ON03dHIhyUcHbid0EkOTdOMQhhBsSEjAAzYYj1iW50my5lmqqvX9cU5J5VKNsqQqSet9Hj9Vtc8+a+86ls5Pa+111hZVxTAMwzAyAU+6J2AYhmEYIUyUDMMwjIzBRMkwDMPIGEyUDMMwjIzBRMkwDMPIGEyUDMMwjIzBl+4JZBolJSVaUVGR7mkYRsaz+0QLhblZlE/OO6O9qaOXI40dzJ9aSI4vvX/3vnusmbLCHKZOzKWhvYdjTZ0snDaRLK+kdV5jkbfeeuu0qpaerR0TpQgqKirYsmVLuqdhGBnPsgd/x/sWT+X//ci5Z7Q/t/049/58G7/8u2uYP7UwTbNzqHzgN9x7/Ty+eNMC1m46zP2/epcXH1jBOUV5iU82UkJEDg2FHQvfGYYxKFQVb5Q7iFccLyQQTP+D+aog7nw8GTQvIzYmSoZhDIqAat+NPhyPJzNu/qFqNaEZhuZlRWwyGxMlwzAGRTAYXZR87s0/mOa7f2h4T5+n5HxO97yM+JgoGYYxKIJKXE/Jn2ZPKSQ+oSla+G50YKJkGMagCCZYUwqmO3znvkaG70yTMhsTJcMwBkUgRvjOmzFrSs5rSIwsfDc6MFEyDGNQqPbf8MPpC5Ol+eYfKT59HpyJUkZjomQYxqBwsu8Gtvu8ofDdCE8oBiGRFMmseRnRMVEyDGNQBILa532EExIBf5rv/gMTHc5sNzKTpERJRFaKyF4RqRaR+6MczxGRp9zjG0WkIuzYA277XhG5OZFNEal0bexzbWbHG0NEikVkg4i0icgjMea/TkR2JHdJDMNIROgZoGjhO2+GpYSHZpgp8zLik1CURMQLPArcAiwGbheRxRHd7gIaVXUe8DCwxj13MbAaWAKsBH4gIt4ENtcAD6tqFdDo2o45BtAFfA34Uoz5fwRoS/Q9DcNInlASQ9REh77U6xGd0gBC0uOJqOhg2XeZTTKe0nKgWlVrVLUHWAusiuizCviJ+/4Z4AZxArirgLWq2q2qB4Bq115Um+45K1wbuDZvjTeGqrar6us44nQGIlIAfAF4MInvaRhGkoRu7N44nlK6s+8iw3eh13TPy4hPMqI0AzgS9vmo2xa1j6r6gWagOM65sdqLgSbXRuRYscaIx7eA/w10xOskIneLyBYR2VJXV5fApGEYkTf8cDJFlPrCd+4kvX1lhkyUMplkRClajffI/9VYfYaqPdl59E9I5AJgnqr+V6w+fUZUH1PVZaq6rLT0rCuvG8aYJyQ40RIdQg/UpjslfEDtO6voMCpIRpSOAjPDPpcDx2P1EREfUAQ0xDk3VvtpYJJrI3KsWGPE4nLgYhE5CLwOzBeRV+N+U8MwkiIkONHCd55MqejQ5yk5r7amNDpIRpQ2A1VuVlw2TuLCuog+64A73fe3AevV+TNlHbDazZyrBKqATbFsuudscG3g2nw2wRhRUdUfqup0Va0ArgLeU9Xrkvi+hmEkIBCILUoZE75zXyMLslr4LrNJuMmfqvpF5F7gJcALPKGqO0Xkm8AWVV0HPA48KSLVON7LavfcnSLyNLAL8AP3qGoAIJpNd8j7gLUi8iCwzbVNrDFcWweBiUC2iNwK3KSquwZ7UQzDiE88T6lPlDKkokOfp5Qh8zLik9TOs6r6PPB8RNvXw953AR+Lce5DwEPJ2HTba3Cy8yLb441RkWD+B4Gl8foYhpE8fWtKmewpRSQ6WPhudGAVHQzDSJm4iQ4ZklAwMNHBebWHZzMbEyXDMFImnqfkyZDKCX1pu5GJDuYqZTQmSoZhpEw8UfJlWPjOE/GckmlSZmOiZBhGysRNCc8QUQpGhO+sosPowETJMIyUiZvokClrSu5rpKdkKeGZjYmSYRgpExIcXyanhIdE0R6eHVWYKBmGkTLxqoRnWkJBZPZdusXSiI+JkmEYKdPnKXnjJTqM6JQGEJnoEHq18F1mY6JkGEbK+ON5Sn2ilGk7z2bGWpcRHxMlwzBSJnTD93mi30K8Hkl7mMw2+RudmCgZhpEy/kBoO/Tox70iaQ/fDax9d2a7kZmYKBmGkTKJPCWPJ/03/8jhMy0Bw4iOiZJhGCnj73tOKfpxn8fT502ljzPXvayiw+jARMkwjJQJ9olSDE9J0u8p9T2mJGe+pnteRnxMlAzDSBl/nCrh4CY6pLuiQ4yUcBOlzMZEyTCMlIlXZshp9/QJV7qIrH3ntTWlUYGJkmEYKZNIlLK8gj/N6XcaEb7re07JNCmjMVEyDCNl4lUJB6fSQ7o9JSWUEu7MUdy7nVV0yGxMlAzDSJlQtYaYnpLHQ2+meEruZ6+tKY0KTJQMw0iZkN5EqxIOrqeU5jhZrESHdD/Ua8QnKVESkZUisldEqkXk/ijHc0TkKff4RhGpCDv2gNu+V0RuTmRTRCpdG/tcm9nxxhCRYhHZICJtIvJImJ18EfmNiOwRkZ0i8u3UL49hGNEIeUqeWKLk8eDPtNp3VtFhVJBQlETECzwK3AIsBm4XkcUR3e4CGlV1HvAwsMY9dzGwGlgCrAR+ICLeBDbXAA+rahXQ6NqOOQbQBXwN+FKU6f+Lqi4ELgSuFJFbEn1fwzASk8hTyvIKven2lNzXkCiFqk+kO1XdiE8yntJyoFpVa1S1B1gLrIroswr4ifv+GeAGcVYXVwFrVbVbVQ8A1a69qDbdc1a4NnBt3hpvDFVtV9XXccSpD1XtUNUN7vseYCtQnsT3NQwjAX2eUoznlHze9HtKqmcmOoT0M90JGEZ8khGlGcCRsM9H3baofVTVDzQDxXHOjdVeDDS5NiLHijVGQkRkEvAh4JUYx+8WkS0isqWuri4Zk4Yxrom382yovdef7ueUnNfQDEUEn0fSvqWGEZ9kRCnaT13kT1usPkPVnuw8BiAiPuAXwPdVtSZaH1V9TFWXqeqy0tLSRCYNY9zTt59SDFHK9nnoTfvNf+CeT15P+hMwjPgkI0pHgZlhn8uB47H6uCJQBDTEOTdW+2lgkmsjcqxYYyTiMWCfqv5rEn0Nw0iC/irhsT2ldN/8I2vfgTsvC99lNMmI0magys2Ky8ZJXFgX0WcdcKf7/jZgvToB3XXAajdzrhKoAjbFsumes8G1gWvz2QRjxEREHsQRr88n8T0Nw0gSf4KKDj5v+p9TCkapz+fzeizRIcPxJeqgqn4RuRd4CfACT6jqThH5JrBFVdcBjwNPikg1jvey2j13p4g8DewC/MA9qhoAiGbTHfI+YK0rKNtc28Qaw7V1EJgIZIvIrcBNQAvwFWAPsNVd7HxEVX+c+mUyDCOcYDJlhtJe+855lXBR8kjaEzCM+CQUJQBVfR54PqLt62Hvu4CPxTj3IeChZGy67TU42XmR7fHGqIgx9ei/MYZhnBWJqoQ7+yllxnNK4bqZCdXLjfhYRQfDMFImmCDRwZcBzykFo9Tny4S1LiM+JkqGYaSMP6gxQ3fg1L5Ld5gs5BGFh++8XvOUMh0TJcMwUiYQ1JiZd5BZte/O9JTSv8+TER8TJcMwUqY3oGR5Y98+sjIg+y7kEdma0ujCRMkwjJTpDQTJ8sbxlDLgeaD+RAfLvhtNmCgZhpEy/mAQXxxPyef1pD18F02UrKJD5mOiZMSkvdvPW4ca0z0NIwPpDShZ8RIdvJL2MkPBqGtK6ffgjPiYKBkx+fxTb/PRH/6JhvaedE/FyDD8gQSekseDanq3iYi2pmQVHTIfEyUjJjuONQPQ1RtI80yMTKM3qPjirSm5x9KZ7NAXvvNEhO9sTSmjMVEyDCNl/IEgWZ542XeOEKQzVBa5HTrgbl1hnlImY6JkGEbK9AaULF+87Dvn1tLrT59XEisl3NaUMhsTJcMwUqY3EOwTnmiEPKV0JjvESgkPidW2w41U3P8baura0jI/IzomSoZhpIw/oHGfUwo9WJvO9Ovoa0r9qepPvnEIgC0HLcM0kzBRMgwjZfzB+J6SLyNEyXn1xvCUcrKcOfakufKEcSYmSoZhpExvIH72XSaE76KuKXn7s+9CohqMv1eoMcKYKBkJsV9ZIxJ/MBi39l3ohp9OT0mjhO/CH54NPVSb7i02jDMxUTIMI2X8gfhVwrN9bmgsI7LvopcZComSPYeXWZgoGQmx7XuNSHoCQbJ8sW8fOSFRCqTvhp9oTSmEiVJmYaJkJMSCG0Yk/gS170KeUndv+lPCJewu5/P276cUqjbR0WOilEmYKBmGkTKJat+FPKXuNIbv+rZDH+ApOXMKhRY7zVPKKJISJRFZKSJ7RaRaRO6PcjxHRJ5yj28UkYqwYw+47XtF5OZENkWk0rWxz7WZHW8MESkWkQ0i0iYij0TM62IRedc95/sSvi+ykTR20YxIeoPxn1PK8XmBdIuS8zpgTck9EEoFt/BdZpFQlETECzwK3AIsBm4XkcUR3e4CGlV1HvAwsMY9dzGwGlgCrAR+ICLeBDbXAA+rahXQ6NqOOQbQBXwN+FKU6f8QuBuocv+tTPR9jYFY+M6IxJ+gokNf+M6fvht+aO0o/E/R8DWlUNadiVJmkYyntByoVtUaVe0B1gKrIvqsAn7ivn8GuMH1SlYBa1W1W1UPANWuvag23XNWuDZwbd4abwxVbVfV13HEqQ8ROQeYqKpvqJMb+tMwW4ZhnAX+BM8pZUL4LpQSfsZ+SmHbtPe4gtlpa0oZRTKiNAM4Evb5qNsWtY+q+oFmoDjOubHai4Em10bkWLHGiDfvownmDYCI3C0iW0RkS11dXRyT4xML3xmR9CZ4TiknI1LCndfw8F2W10NvQFHVME/JKjpkEsmIUrR7UmREJ1afoWpPdh7JzGlgo+pjqrpMVZeVlpbGMTk+sfCdEUlvgtp3mbGmNLCiQ3+qetASHTKUZETpKDAz7HM5cDxWHxHxAUVAQ5xzY7WfBia5NiLHijVGvHmXJ5i3YRgpoqoEghp3TamvrlyaRUkEwvObsr3987JEh8wkGVHaDFS5WXHZOIkL6yL6rAPudN/fBqx313HWAavdzLlKnGSDTbFsuudscG3g2nw2wRhRUdUTQKuIXOauVX06zJZhGIMkFPaK5ymFbv7pTHQIqp6RDg5hNfkC2re2ZKKUWfgSdVBVv4jcC7wEeIEnVHWniHwT2KKq64DHgSdFpBrHe1ntnrtTRJ4GdgF+4B5VDQBEs+kOeR+wVkQeBLa5tok1hmvrIDARyBaRW4GbVHUX8NfAvwN5wAvuP8MwzoK+gqZx1pQ8HiHLK2kN3wWCZ64nAWS7YcUef7BPlCx8l1kkFCUAVX0eeD6i7eth77uAj8U49yHgoWRsuu01ONl5ke3xxqiI0b4FWBrtmGEYgyMUksuJU2bIOe5Na/hOVYmMMIbX5PNbokNGYhUdDMNIiZD3k51AlLJ9nrQ/pzTQU+pPdAg9RGueUmZhomQYRkr0e0reuP1yfJ40176LEr5z15QcT6m/3FBkkVYjfZgoGQmJk09ijENC3k8iTynH50nrrq5BVSJrxoZ7SuH7KI2ER7f5YAPXfGcD1bWtwz7WaMZEyUiIaZIRTmgNJtGaUnaaPaVAUM+o5gCQ7XW8u95AsC9hA0amqsOzbx/jcEMHv911atjHGs2YKBmGkRIh7yexp+RN65qSP6gDMgTDEx0CQe1LXR+JdaVD9R0AVNe2DftYoxkTJcMwUqI7SU8p3eG7QDA4YHfcrLA1pd6AUpjrJCCPRAZeU0cvACeauhL0HN+YKBkxsbCdEY2Q0GR6+M4fLXwXnn0XCFLQJ0rD7yk1d7qi1Nw57GONZkyUjISYOBnhdLs38GSy77rSnBIe6SmFF4rtDSoFOSMnSi1dIVHqsuShOJgoGTEJZdOqlWQ1wkh2TSkv25vWbSGieUpZ3jPXlELhu+FeUwoGlZbOXrK9Hrr9QduCPQ4mSkZM7I85IxrJrinlZfnSK0pRNiIMD98FgkpBThYw/Nl37T1+ggqVJRMAaGjvGdbxRjMmSkZCTJyMcJL1lPKzvXSksVpCIDhwI8JQtl3IU+lLdBjmckih9aSQKJ1u6x7W8UYzJkqGYaREsmtK+TnetIap/FHWlEJC2tHt7CPat6Y0zPNs6XTGqyw1TykRJkpGQsxRMsJJNvsuP8uX1hI+0R6eDa0phTy4kVpT6vOUih1RqjdRiomJkpEQyxQywgmtKSUTvgPo6PEP+5yi4Q8M3IgwJKTtIU9phFLCQ5l3IU+pvs1EKRYmSoZhpERPIIhHGBAaiySvT5TSE8LzB4MDPCURITfLQ2uXI0qFOSPrKU2bmEtuloeGdltTioWJkpEQ85OMcLr9QbJ9njO2GY9GftpFaWCiA0B+to8WVySyfR6yvZ5hr+gQGm9iXhbFE3IsfBcHEyUjJvZ8khGNHn8wYZIDpD98F+3hWYC8LG9fOM3n8ZCb5Rn+8F1nLyKOZzZlQraF7+JgomQkxJaUjHC6egMJ15PA8UhgZCpwR8MfULyRW8/iiGUofOfzCnnZ3mEXzpYuP4U5Pjweobggm3oL38XERMmIiRD6K9NUyeinszfQ5wXFI93hu1ieUn62ty+c5vN4KMzN6hOp4aK5s5eifOdB3ZKCHE63mqcUCxMlIyYWvjOi0dETIC8rsSjlpTv7LhjEG2VNKS/bS4srQl6PMDHX1xfOGy5aOnuZmNsvSvXt3ZbVGoOkRElEVorIXhGpFpH7oxzPEZGn3OMbRaQi7NgDbvteEbk5kU0RqXRt7HNtZp/FGH8nIjtFZIeI/EJEclO7PAZY+M44k86eQJ/gxCMUvktrokOMNaU2NyU8J8tDUV5W38Otw0VzZy9FeSFRyqY3oMM+5mgloSiJiBd4FLgFWAzcLiKLI7rdBTSq6jzgYWCNe+5iYDWwBFgJ/EBEvAlsrgEeVtUqoNG1PZgxZgB/CyxT1aWA1+1nGMZZMFrCd9GeU4J+sQTnuaWJeVl9KdvDRUtXv6dUWpgDQJ2VGopKMp7ScqBaVWtUtQdYC6yK6LMK+In7/hngBnHyRVcBa1W1W1UPANWuvag23XNWuDZwbd46yDEAfECeiPiAfOB4Et/XiMAcJSMcJ3znS9gvY7PvwgQ1x+d1PKVhDt81d/YyMc+5ZsUTHFGy+nfRSUaUZgBHwj4fddui9lFVP9AMFMc5N1Z7MdDk2ogcK6UxVPUY8C/AYeAE0Kyqv432BUXkbhHZIiJb6urqYl6I8YqF74xwOnv8SXlKE7J9iEDbMCcRxMIf1KhrSvlniJKHiblZtHT2EhzGckgtnf7+8F1hNmCiFItkRCnaE3KR/3ux+gxVe8pjiMhkHC+qEpgOTBCRT0Xpi6o+pqrLVHVZaWlptC6GYbgkm+jg8QgFOb6+pIKRJtp26HCmp5TrrikF1dleYjjo8Qfp7A2ckegAcLrVRCkayYjSUWBm2OdyBobB+vq4obIioCHOubHaTwOTXBuRY6U6xo3AAVWtU9Ve4FfAFUl8XyMCy8IzwunsTS7RAXCTCIY3NBYL5zmlKJ5SVviakrcvrDZc60ohu6GU8Mn52XgETtsDtFFJRpQ2A1VuVlw2TrLAuog+64A73fe3AevVyXdcB6x2M+cqgSpgUyyb7jkbXBu4Np8d5BiHgctEJN9de7oB2J3cZTEMIxadPcklOgBOaGyY12ti0R0IRn3ItyjvzESHUFht2EXJHcfrEaZMyLHwXQwSrlaqql9E7gVewslge0JVd4rIN4EtqroOeBx4UkSqcbyX1e65O0XkaWAX4AfuUdUAQDSb7pD3AWtF5EFgm2ubQYyxUUSeAba67duAxwZ7ocYztqZkhOjxB/EHNanwHcDEPF9aUp9V1SmH5B0oSpMnZPe9z/Z5mOImHgzXHkfNYXXvQpQUZJunFIPEKTSAqj4PPB/R9vWw913Ax2Kc+xDwUDI23fYa+rPnwtsHM8Y/Av8Y7RwjeUyUjBChatrJhu8m5mZxuKFjOKcUld6A80MbzVOaEiZKE3J8fSnatS3D47m0RHhK4KSFm6cUHavoYBhG0oTq2IU/6xOPiWlaU4q3Zfvk/H5RyvJ6KAuJ0jAlHkSG78ARpdqWrmEZb7RjomQkxBIdjBChZ47yspO7dThrSiMfvuvxO6KUFSV8F+4pgeMtTcj2UjfMojQpTJRmTMrjZEsX/sDwbpkxGjFRMmJiYTsjklB1hmQengVnTamt2z/iN9+QKEXzlELhunDKJuZS2zo8nku0NaUZk/IIKpw0b2kAyf1kGeMaEycjRKia9sTcJEXJfTanrdvPpPzsBL2Hjj5RiuIpZXk9fGvVEqZO7C+FWVqYM2zhu6aOXiZke8/w2qZPygPgeFMX5ZPzh2Xc0YqJkhGTBBuLGuOQVje9uzA3K0FPh/B06xEVpYDj0cXa9+nPL68443NZYQ7vHmselrmEF2MNMWOyI0rHmjqAKcMy7mjFwndGTMxDMiIJeUqFSXpKkyc4N+PhSreORbfrKeUksRkhwKwp+Rxr7ByWMGNTRw9FEYI8vajfUzLOxETJSIiJkxEitOVDsqIUKj460tt/x1tTisbs4nz8QeVE89CLRF1bd1+GX4i8bC/FE7I52tg55OONdkyUDMNImlTDdyWF6amI3feckje556lmF08A4GB9+5DPpa61O2pyxYzJeRxtHPlnuDIdEyUjIZYSboRo7fKT4/Mk7YEUu+nX9SMcvhuMpwRwqH5oRSIYVE63RRelypIJ1NQNvQiOdkyUjIRY+M4I0dLlTzp0B5Cb5aUwxzdszwDFIlGiQyRTC3PJzfIMuUg0d/bSG9AB4TuAeaUFHGvq7Hsg2XAwUTIMI2lau3qTDt2FKElDSZ14KeHR8HiEBdMmsuvE4DLw/v2PB/j2C3sIROzJFEozD21XEc7csgIA9te1DWrMsYqJkpEQc5SMEK0pekrghPBGOtGhq9fNvstK/ha3dPpEdh5vQVMMDVTXtvGNX+/iR6/t56nNR844dsSt+1fupoCHM7fURCkaJkpGTEyMjEgcTyk1USopGHlPqaOvRl9yiQ4AS2cU0drl50hDahlxf9jn7FZdWpjDj1+vOWMH20OuKFW4iRThzC7OxyOw39aVzsBEyUhIqn85GmOXli5/X5WGZBnOagmxCNXoS7ZwLMDS6UUAvH20KaWxqmvbKMrL4h/ev5Caunb+uP9037FD9e0U5vqYlD/wmuVmeakonsCeEy0pjTfWMVEyYhIq6GCSZIRoaO8ZUNA0EdMn5dHc2duXTj4SdA7CU1p0TiGFOT7e2F+f0lj7atuYV1bA+889h+IJ2fz0jUN9x9471cqc0gIkRnmUpTOK2DFMlSRGKyZKRkxMjIxwAkGlqaOnL807Wcr7SuqM3IOiHb0BsrwStUp4LHxeD5fOmcKfwjydZNhf20ZVWQE5Pi+3L5/FK7tPcaShg0BQ2XGshfPLi2Kee155Ecebu0Y8OzGTMVEyEmLROwOc9OagnrlzazL01XkbweoFnT2BpHfHDeeKuSUcqu/oS1BIREN7D/XtPcxzM+nuuHQWAD/beJjtR5to6/Zz0azJMc8/d4YjWO8eSy1kOJYxUTIMIylC9etSDd+VT0qDp9TjT2k9KcSNi6YC8Py7J5LqX13rZM6F0runT8rjpsXT+Nmbh1jz4h6yvR6uX1gW8/ylM4oQgXeOWAgvhImSkQTmKhmDF6WSghyyvZ4R9ZQ6egIprSeFmFWcz/nlRTy3PTVRqnJFCeArH1hEts/DmzUN3H3NnAEVwsOZkONjwdRCNh9sSDjWgdPt1IyD9HETJSMhFr4zYPCi5PEIMybncWQE67x19gTIG4QoAXzo/Om8e6yZvSdbE/bdV9tKXpa3r+o3wMwp+bzyxWtZd++VfPGm+QltXDG3hLcONdLVG7uyQ2dPgA98/w+seuSPY3632qRESURWisheEakWkfujHM8Rkafc4xtFpCLs2ANu+14RuTmRTRGpdG3sc21mn8UYk0TkGRHZIyK7ReTy1C6PYRghBitKAHNKJvR5FSPBYD0lgI9eVE6Oz8O//+lAwr7VtW3MLZuAx3Nmdt2k/GzOK58UM+sunCvmFtPtD7LtcOx1pbePNNHRE6C128/bR8b2+lNCURIRL/AocAuwGLhdRBZHdLsLaFTVecDDwBr33MXAamAJsBL4gYh4E9hcAzysqlVAo2s75THcc74HvKiqC4Hzgd3JXhijH3OUDICGdidDbPIgNuurmlrIgdPt9I7QX/mt3amXQwoxeUI2H7monF9tPZbwod/q2jbmlRbE7ZOI5XOm4BF4I07WX3jixTtHx/b6UzKe0nKgWlVrVLUHWAusiuizCviJ+/4Z4AZx/kRYBaxV1W5VPQBUu/ai2nTPWeHawLV562DGEJGJwDXA4wCq2qOqY/tPDMMYRk62dDE5P4vcQWS1VZUV0BvQIa/CHYumjoG7vabCZ6+uxB9Uvv/Kvph9Wrp6OdHcxfxphYMeB5wt48+fOYlX36uL2edQQztej1BSkM27KT7cO9pIRpRmAOEFnY66bVH7qKofaAaK45wbq70YaHJtRI6V6hhzgDrg/4rINhH5sYgMrPUBiMjdIrJFRLbU1cX+wRiv2JqSAXCyuYupE3MHde78qc6Ne9+pxOs0Q0G0LchTYU5pAbcvn8nPNx6OWZsu9F3ml52dKAG8b/FUth9t5kRz9GSQww2dzJiUx3nlk9g1xitAJCNK0YKikbepWH2Gqn0wY/iAi4AfquqFQDswYD0MQFUfU9VlqrqstLQ0WpdxjZUZMsDxlM4pGpwozSsrQAT2JJE8cLYEgkprl/+sRAngczfMJy/by/2/3H5GPbsQe086YrXgLD0lgJsWTwPgd7tORT1+uL6d2cX5LDqnkP117XGTIkY7yYjSUWBm2Ody4HisPiLiA4qAhjjnxmo/DUxybUSONZgxjqrqRrf9GRyRMgxjEJxs7mLaIEUpL9vLgqmFbMoObPAAACAASURBVBuBRfqWTqec0dmKUmlhDv/4oSVsPtjIE38cmPTw7rFmCnN8zJg0sAJ4qswrK2BO6QR+G0uUGjqYOSWfRedMJBBU9p0au6nhyYjSZqDKzYrLxkkqWBfRZx1wp/v+NmC9On9erwNWu5lzlUAVsCmWTfecDa4NXJvPDmYMVT0JHBGRBe45NwC7kvi+RgTmJxnd/gCn23qYNnHwN+ALZ01m2+HGqF7HUNLsilK0Iqip8tGLZvC+xVP59gt7eLPmzJp4mw7Us6xi8oDMu8Gycsk0/rS/fkDJoZauXho7epk1JZ/F50wEYPcYDuElFCV3/eZe4CWc7LWnVXWniHxTRD7sdnscKBaRauALuGEyVd0JPI0jBi8C96hqIJZN19Z9wBdcW8Wu7ZTHcM/5G+BnIrIduAD4p1Qv0HjGonZGiNoW50Y52PAdwEWzJtHa5R/2/YMaO5zU9bP1lABEhP/98fOZVZzPX//HW1TXOuHHA6fb2V/XzhVzS856jBB/duEMAkHl2bePndF+2E0OmTUln9nFE8jL8o7pdaWk6nCo6vPA8xFtXw973wV8LMa5DwEPJWPTba/Byc6LbB/MGG8Dy6KdYySPiZMRevB1+lmEqpZVTAHgzZp6qqae/TpMLELbZJQVDl5Aw5mYm8WPP72MTzz2Jh///97kOx89j+d3nMDrET50/vQhGQOctPnzy4v45dZj/OXVc/raD7vp4LOL8/F6hAXTCse3p2QYagG8cc/B0+5mdSX5g7ZRUZzP7OJ81u+pHappReVUSxcAU4sGbkE+WOaUFvCff3U5k/Ky+MufbuFXW4/x2avnDHqNLRYfvbic3Sda2HW8X3RCafSz3Y0CF50zkd0nUt8hd7RgomQYRkIO1beT7fOcUU4nVUSEFQvL+OP++r5N+IaDUy1dzjM9E4ZOlAAqSibwwuev5pE7LuTxO5dx38oFiU9KkQ+fP50cn4f/2Ni/J9PhhnaKJ2RTkOMEthafU0hLl5/jzV1DPn4mYKJkJGZs/kFmpMCB0+3MnpJ/1ov671s0lR5/kFd2D5+3dLK5m7LCnCFLQAgnx+flg+dN54ZFU5MqIZQqk/Kz+chFM/jlW0f7yjrV1Dnp4CEWhZIdjo/NEJ6JkmEYCTlY394XPjobLptTzIxJeTy95UjizoPkSGPHkKRpp4u/uLKSbn+Qn715iGBQ2XW8hSXT+zcKXDjGM/BMlIyEmKM0vukNBDl4uoO5pWcvSh6P8PFlM/nDvtPDtg1DTV0bc8+yHl06qZpayA0Ly3jsDzW8Xn2a1m4/54XtXluQ42PWlHx2nzRRMsYdJkcG7K9roycQZPH0iUNi745LZ5GX5Y1bV26wNHX0cLqth7llZy+g6eSB9y+iqzfAp5/YhNcjAzYKXHzORHafGJmSTSONiZKRkDGa5GMkSSgTLPTg5tlSWpjDp6+YzbPvHB/ybRhCN+qqIahHl07mlRXwz7edz6wp+dy/ciElBWcmbSw6ZyIH69uHNWEkXZgoGXFwFnItJXx8s/N4C7lZHuYMYUjsnuvnMW1iLl/6z3fo7Bm6Om5vHXJ2cL1o1uQhs5kubr1wBr//8vV89po5A44tOqcQVcakt2SiZMTBxMhwarwtmDYR7xBms03MzWLNR89jf10bf7t225DtpvrH6noWTC2kaAhKDGUy55VPAmDb4cY0z2ToMVEyEmLhu/FLV2+At480ccnsofc8rplfyj9+cDEv7zrFXz35Fm3dZxeKqm3tYuOBem5aMnWIZpi5TCvKZXZxPhsPNKR7KkOOiZJhGDF550gTPf4gl84pHhb7/+PKSr5161I27K3l5od/z8u7Tg26UsHjrx9AccJe44HlFVPYfLBh2AvcjjQmSkZCxtaPvJEKGw80IAKXVAzfGs2fXzabp//qcnKzPHz2p1u45Xt/4PHXD3C8KfqGd9F4Y389j//hALdeMGNUp4OnwqVzimnq6GXvCG2cOFIkVZDVGN+M1RpbRmLW76ll6fQiJuVnD+s4yyqm8MLnrmHdO8d54vUDfOu5XXzruV2UT87jwlmTmVdawKziPM4pyqMgx0d+tpfO3gAnm7vYsLeWpzYfYVZxPt/40JJhnWcmcWmlU+B204GGvioPYwETJcMwonKqpYu3jzTxpZvmj8h42T4Pt11czm0Xl1NT18b6PbVsO9zE1kONPLf9eMy1zWyvh49cWM4D71845hMcwimfnMf0olze2F/PnVdUpHs6Q4aJkhGTUKja/KTxyW93ngTgpiXTRnzsOaUFZ6Sgd/UGONbUyYmmLtp7/HT2BMjxeSgpzGHp9CLysr0jPsd0IyJcXVXK8ztO4A8E8XnHxmqMiZIREwvbjW+e3nKUBVMLqSpL/xpNbpaXuaUF42a9KFmuXVDKU1uO8PaRpr79qkY7Y0NajWFBB7wxxgvvHm3m3WPN3HHprGGphm0MDVfOLcEj8Np7demeypBhomTERPvCd6ZK443HX68hL8s7btKrRytF+VlcOGuyiZIxPrDw3fikuraNde8c59OXz6Yob/wkDoxWrp1fyrvHmqlv6073VIYEEyUjJn2ekmnTuOLbL+whx+eNWnPNyDyunV+KKrxefTrdUxkSkhIlEVkpIntFpFpE7o9yPEdEnnKPbxSRirBjD7jte0Xk5kQ2RaTStbHPtZk92DHcY14R2SYizyV/WQywpaTxyEs7T/K73af4/I1VAypTG5nJuTOKmDIhm1f3jo0QXkJREhEv8ChwC7AYuF1EFkd0uwtoVNV5wMPAGvfcxcBqYAmwEviBKxLxbK4BHlbVKqDRtZ3yGGFz+xywO7nLYYQTCt+ZpzQ+ONncxVf+610WTivkL66qTPd0jCTxeIRrqkp47b26MVFyKBlPaTlQrao1qtoDrAVWRfRZBfzEff8McIM4KTurgLWq2q2qB4Bq115Um+45K1wbuDZvHeQYiEg58AHgx8ldDiMcjXg1xi5dvQHu+flWOnoCPHLHhWSNkWdexgvXLyyjob2H7cea0z2VsyaZn7wZwJGwz0fdtqh9VNUPNAPFcc6N1V4MNLk2IsdKdQyAfwW+DMStiy8id4vIFhHZUlc3NlzgocA8pPFBbyDIvT/fyluHGvnObecxb5RvkDceubqqFBHYsKc23VM5a5IRpWgPKUTermL1Gar2lMcQkQ8Ctar6VpTjZ3ZWfUxVl6nqstLS0kTdxw3BvvCdqdNYpbMnwD0/28rvdtfyrVVL+OB509M9JWMQTJmQzQUzJ/HqGEgNT0aUjgIzwz6XA8dj9RERH1AENMQ5N1b7aWCSayNyrFTHuBL4sIgcxAkPrhCR/0ji+xouJkVjm9rWLlb/25u8vPsU//ihxfz55RXpnpJxFly/oIztR5tGfWp4MqK0Gahys+KycZIK1kX0WQfc6b6/DVivzp/X64DVbuZcJVAFbIpl0z1ng2sD1+azgxlDVR9Q1XJVrXDtr1fVTyV5XQzoUyUTp7HH+j2nuOVf/8Deky386FMX85krLbFhtHPdAic1/Pf7Rre3lLD2nar6ReRe4CXACzyhqjtF5JvAFlVdBzwOPCki1Tjey2r33J0i8jSwC/AD96hqACCaTXfI+4C1IvIgsM21zWDGMM6OUCUHi96NHZo6eljz4l5+sekwC6cVsvb2y6iaamtIY4Gl04soKchmw546/uzC8nRPZ9AkVZBVVZ8Hno9o+3rY+y7gYzHOfQh4KBmbbnsNbvZcRHvKY4QdfxV4NdZxIzomRmOHYFD5z7eO8O0X9tDc2ctfXlXJl25eQG7W+KuuPVbxeIRr5peyfk8tgaDi9YzOmoVWJdyISVAtKXy0o6q8vOsU3335PfacbGXZ7Ml8c9VSFk8fO5vCGf1cv6CMX209xttHmrh49vDtFjycmCgZMTEpGr2oKq++V8fDL7/H9qPNVBTn873VF/Dh86db1e8xzDVVpXgEXt1ba6JkjD2s9t3oo9sfYN3bx3n89QPsOdnKjEl5fOe28/jIhTPGzCZwRmyK8rO4aNZkNuyt5Ys3LUj3dAaFiZKRENOkzKepo4efbTzMv//pIHWt3SycVsg/33Yeqy6YQbbPxGg8cf3CMv75pb3UtnZRVpib7umkjImSYYxSVJW3DjXy842H+c27J+j2B7lmfinf/XglV80rsTDdOOW6BaX880t7eW1vHR9bNjPxCRmGiZKREAvfZRbNHb38attRfrHpMO+daqMwx8fHl83kU5fNZsE0S+8e7yw+ZyJlhTm8+p6JkmEYw0QwqGw51MhTm4/w3PbjdPuDnF9exJqPnsuHzp9Ofrb9KhsOIsJ1C0p5YcdJ/IHgqFtLtJ9kIyG2HXr6OHi6nV9tO8Z/bTvKkYZOJmR7+ejF5dyxfBZLZxSle3pGhnLdgjKe3nKUrYebWF45Jd3TSQkTJSMhFr4bWZo7e/nN9hP8autRthxqRASunFvC3904n5VLp5lXZCTkqqoSvB7h1b21JkqGYaRObyDIH/bV8cutx3h51yl6/EHmlRVw38qF3HrhdM4pykv3FI1RxMTcLJbNnsyGvXV8eeXCdE8nJUyUjISYozQ8BIPK5oMN/Hr7cZ5/9yQN7T1Mzs/ijuWz+MhFMzh3RpFl0BmD5roFZax5cQ8nm7uYVjR6UsNNlAxjBFFV3jnazK/fOc5vtp/gZEsXuVkeblw0lQ+fP53rFpTZc0XGkHD9wlLWvLiH196r5ROXzEr3dJLGRMlIiG3yd3aoKntOtvLrd47z6+3HOdLQSbbXw7ULSvmH8xdxw8IyJuTYr6IxtCyYWsg5Rbls2FNnomQYBtTUtfHrd07w6+3Hqa5tw+sRrpxXwt+uqOKmJdMoystK9xSNMUwoNfzX75ygNxAka5SkhpsoGcYQoarsq23j+XdP8OKOk+w52YoIXFIxhQdvXcotS6dRXJCT7mka44hrqkr5xaYjvHOkiWUVoyMLz0TJiEp4yM6id7FRVXYca+GFHSd4cedJauraEYGLZ03mqx9YxAfOO8cy54y0cdmcYkTgjf31JkrG6CZoQhSTYFDZdqSJF3ec4IUdJzna2InXI1xaOYXPXFHBzUumUTZx9GQ7GWOXyROyWTRtIn/aX8/f3FCV7ukkhYmSEZXeQLDvvVV0gEBQ2XSggRddj+hUSzdZXmeN6G9WzON9i6cxZUJ2uqdpGAO4fG4xT755iK7ewKjYadhEyYiKP2jhu/ZuP69Xn2b97lp+t/sU9e095Pg8XDu/lFvOncaKhVMtWcHIeK6YW8zjrx9g6+FGrphbku7pJMREaYzz5JuHOHdGERfMnJTSef4wT2k8cbi+g/V7TvHKnlo21jTQEwhSmOPj2gWl3LL0HK5bUGrp28aoYnnlFLwe4Y399WNHlERkJfA9wAv8WFW/HXE8B/gpcDFQD3xCVQ+6xx4A7gICwN+q6kvxbIpIJbAWmAJsBf5cVXtSHUNEZrr9pwFB4DFV/V6qF2i087X/3gHAwW9/IKXzOnsDfe/HsqfkDwR561Aj6/fU8sqeWqpr2wCYUzqBO6+YzYqFU1lWMXnUpNMaRiSFuVksnVHEG/vr0z2VpEgoSiLiBR4F3gccBTaLyDpV3RXW7S6gUVXnichqYA3wCRFZDKwGlgDTgd+JyHz3nFg21wAPq+paEfmRa/uHgxjDD3xRVbeKSCHwloi8HDHvMU3wLLIVOnsCiTuNUhrbe3jtvTpe2VPLa3traenyk+UVLq0s5o7ls1ixsIyKkgnpnqZhDBlXzC3m335fQ3u3P+M9/WRmtxyoVtUaABFZC6wCwm/uq4BvuO+fAR4Rp2jXKmCtqnYDB0Sk2rVHNJsishtYAdzh9vmJa/eHqY6hqm8AJwBUtdW1PSNi3mOajt7BC8sZntJQTCaNqCrvnWrjlT2n2LCnlrcONRJUKCnI5qYl07hhYRlXVZVQmGvrQ8bY5PI5xfzw1f1sPtjAdQvK0j2duCQjSjOAI2GfjwKXxuqjqn4RaQaK3fY3I86d4b6PZrMYaFJVf5T+gxkDABGpAC4ENkb7giJyN3A3wKxZo6ccRyLauvyJO8Ug3FMajWWGunoDvFFTz4Y9tbyyu5ZjTZ0ALJk+kXuvn8eKRVM5b0YRHo8VPDXGPk4IWnijpn5MiFK039rIu1SsPrHaowXo4/UfzBjOSSIFwC+Bz6tqS5S+qOpjwGMAy5YtG3134Bi095yFKJ2Fl5UuTjZ3sWGvI0J/rD5NZ2+AvCwvV84r4d4V87h+QdmoqpZsGENFfraPC2ZOGhXrSsmI0lEgfKP3cuB4jD5HRcQHFAENCc6N1n4amCQiPtdbCu+f8hgikoUjSD9T1V8l8V3HFB3dgxeWjp7MD98Fg8r2Y82s3+1ky+087vzNMWNSHh9bVs6KhWVcNqd4VDybYRjDzaWVxfzwtf0Zv66UzMw2A1VuVtwxnKSCOyL6rAPuBN4AbgPWq6qKyDrg5yLyXZwkhCpgE453M8Cme84G18Za1+azgxnDXW96HNitqt9N9cKMBdq6B+8pNXX0DOFMho7Wrl5e33eaV/bU8ureWk639eARuHj2ZO5buZAVC8uYP7XA9iEyjAguqZzCIxuq2Xq4kaurStM9nZgkFCV3/eZe4CWc9O0nVHWniHwT2KKq63Bu/k+6SQYNOCKD2+9pnOQCP3CPqgYAotl0h7wPWCsiDwLbXNukOoaIXAX8OfCuiLzt2vgHVX1+cJdq9NHhhu+yB5HO3NDe2/8hza7SwdPtvLKnlvV7TrHpQAO9AWViro9rF5Rxw8Iyrp1fymSrpmAYcblo1iQ8ApsPNIxuUQJwb+TPR7R9Pex9F/CxGOc+BDyUjE23vYb+DL3w9pTGUNXXib7eNG4IeUqD2TSuMcxTGukyQ72BIJsPNrB+dy3r99ZSU9cOwLyyAv7iykpWLCzj4tmT8dmzQ4aRNIW5WSyZXsTGAw3pnkpcMjewaJw1IVHKGYQonWrpwusRAiNUmbW+rZtX99axfm8tv99bR2u3n2yvh0vnTOHTlzkPsc4qzh+RuRjGWOWSiin8bOMhuv0BcnyZudZqojSGqW9zvJ2Jg6jPdvB0OxXF+eyvax+Wig6h3VjX76nlld2n2HakCVUoLczh/eeew4pFZVw1rySjF2QNY7SxvHIyT/zxADuONXPx7MzcysJ+48cwda3dgzovGFRq6tq5dE4x+93Q2VDQ1RvgT/tP88ruWtbvqeVEcxcA55UX8bkbqlixsIyl0+3ZIcMYLkJ7Km060GiiZIw8oZt+b4rFVfeeaqW128+Fsybxu92nzmpFqba1y62yXcvr1XV09QbJz/Zy1bwSPn9jFdcvKLO9hwxjhCgpyGFu6QQ2H2zgr5mb7ulExURpDPPeqVYA/IHUZOXlXacAp14WpFaQNVTS53e7T/HyrlO8faQJcJ4d+viymdy4aCqXzpmSsfFswxjrLK+cwnPbTxAIKt4MjEqYKI1RTjZ3cbihAwB/MHlPqb3bz3+8eYir5pUwfVJy23iH1oee236c57af4FC9M+755UV88X3zuXHxVBZOK7RnhwwjA7ikYgq/2HSEvSdbWTx9YrqnMwATpTHKv/2hBoCrq0p4x/VWEhEMKl/97x3Utnbzw09d3JdPH4zhKtW1dvP0liP817ZjVNe24fUIV8wt5u5r5nDjoqlMtbCcYWQcl7jrSpsPNpgoGcOLqrLxQAPff2Uff9pfzycvnUVelpe3DjUmPLelq5f7ntnOCztO8sX3zefi2ZNpbHey9yI3/HvvVCuPrK/mhR0n6A0oyyum8K1bl3LL0mmUFOQMy3czDGNoKJ+cx/SiXDYdbODOKyrSPZ0BmCiNAVSVV/fW8eiGarYcaqSkIJuvfmARn7mykn/57d64a0qdPQF++sZBfvTaflq6/Hzl/Yv4y6srgf6Hbrv9jii1dvXy0G9289SWI0zI9vGpy2bzyUtnM6+sYNi/o2EYQ4OIcEnlFP60vx5VzbiwuonSKCYQVF7ccZIfvFrNzuMtTC/K5X99eAmfuGRmXxHSLI/QGwwO+OHr9gdYu+kIj2yopq61m2vml/L3Ny3g3PKivj45YaJ0pKGDTz2+kSMNHdx1ZSX3XD/PSvsYxijlkoopPPv2cQ7Vd2TchpYmSqOQHn+Q/952jB+9tp+a0+1UlkzgO7edx60XzBhQUsjn9aDqCJjP61Ro+OXWo3zvd/s41tTJ8sopPHrHRSyvHPjMgs/rwesRGjt6+PQTm2hs7+Gpv7q8LyZtGMboJPT7vulgg4mSMXhq6tp4avMRnnnrKPXtPSyZPpFH77iIlUunxUztzHLrw/mDytbDDXz1v9/lvVNtnF9exLc/ei5XzSuJ675nez383z8eRATWfvYyEyTDGAPMKy1gcn4WG2sa+PiymYlPGEFMlDKc2pYuXtp1iufeOc7GAw14PcINC8v41GWzuboqvqAAZHmd4//x5iH+6fndlE/O54efdIQsmVhyl9/ZV2nV+dO5dE7x2X8hwzDSjscjXDanmDdrMm9dyUQpAzlc38GLO0/w4o6TbD3spHPPKZnA39+8gI9dXJ5SBQSf60E9+JvdXDO/lB988iIKUqgnF8oG/+w1c5L/AoZhZDyXzy3mhR0nOdLQmVHFjk2UMoBQFYQXd5zkxZ0n2X3C2UF1yfSJfPF981m5dBrzyga3cV1W2BrT/7n9wpQECeBrH1zMscZOlkwvStzZMIxRw+Vu5OONmtPMKp6V5tn0Y6KUJlSVd4428+KOk7y08yQHTrcjAstmT+arH1jEzUumMXPK2f/1Mn9qIQD/87q5FA2iWvhdV1We9RwMw8g85pUVUFKQwxv76/nEJSZK45JAUHnrUCPPv3uCl3ae5ERzFz6PcPncYv7y6kret3gqZYVDWwXhkoopvPj5q5lXas8SGYbRj4hw2ZwpvJFh60omSiPAvlOtPLX5CM++c5y61m6yfR6uqSrl729ewA0Lp1KUn7oHkwoLp2VeKRHDMNLP5XOLeW77CQ6cbmdOhvzhaqI0TASDyu92n+JHr+1n6+EmfB7hhkVlfPC86Vy/sCzltR3DMIyhpn9dqd5EaSyz5WADX3t2J7tPtDBzSh5fef8i/uyiGVYXzjCMjKKyZAIzJuWxYU8tn7x0drqnA4AncRcQkZUisldEqkXk/ijHc0TkKff4RhGpCDv2gNu+V0RuTmRTRCpdG/tcm9lDPcZw4Q8E+afnd3Pbj96guaOH7378fDZ88To+e80cEyTDMDIOEeHmJdP4/b7TtHb1pns6QBKiJCJe4FHgFmAxcLuILI7odhfQqKrzgIeBNe65i4HVwBJgJfADEfEmsLkGeFhVq4BG1/ZQjzHk9AaC/PXPtvLY72v45KWzePkL1/KRi8rxeZPSfcMwjLRwy7nT6PEHeWV3bbqnAiTnKS0HqlW1RlV7gLXAqog+q4CfuO+fAW4QJ5VjFbBWVbtV9QBQ7dqLatM9Z4VrA9fmrUM5RnKXJTWCQeXLz2zn5V2n+MaHFvPQn53LBFszMgxjFHDxrMlUlkzge6/s43Rbd7qnk9Sa0gzgSNjno8Clsfqoql9EmoFit/3NiHNnuO+j2SwGmlTVH6X/UI0xABG5G7jb/dgmIvXA6Wh9E/GZNfCZwZyYmZQwyOswBrFr0Y9dC4cxdx1K/37Qp5YAQ7IolYwoRUtej9ygJ1afWO3RPLR4/YdyjIGNqo8Bj4U+i8gWVV0Wre94wq5DP3Yt+rFr4WDXoR/3WlQMha1kwndHgfAysuXA8Vh9RMQHFAENcc6N1X4amOTaiBxrqMYwDMMwMpRkRGkzUOVmxWXjJBWsi+izDrjTfX8bsF5V1W1f7WbOVQJVwKZYNt1zNrg2cG0+O5RjJHdZDMMwjHSQMHznrt/cC7wEeIEnVHWniHwT2KKq64DHgSdFpBrHe1ntnrtTRJ4GdgF+4B5VDQBEs+kOeR+wVkQeBLa5thniMRLxWOIu4wK7Dv3YtejHroWDXYd+huxaiGrUZRbDMAzDGHHsIRrDMAwjYzBRMgzDMDIGE6UwRrosUToQkSdEpFZEdoS1TRGRl93STi+LyGS3XUTk++712C4iF4Wdc6fbf5+I3BltrExGRGaKyAYR2S0iO0Xkc277eLwWuSKySUTeca/F/3Lbh6zk12jCrQizTUSecz+P1+twUETeFZG3RWSL2zb8vx+qav+cdTUvsB+YA2QD7wCL0z2vYfie1wAXATvC2r4D3O++vx9Y475/P/ACzrNglwEb3fYpQI37Otl9Pznd3y3F63AOcJH7vhB4D6cc1Xi8FgIUuO+zgI3ud3waWO22/wj4a/f9/wR+5L5fDTzlvl/s/t7kAJXu75M33d9vENfjC8DPgefcz+P1OhwESiLahv33wzylfkasLFE6UdXf42QvhhNewimytNNP1eFNnGfIzgFuBl5W1QZVbQRexqk7OGpQ1ROqutV93wrsxqkEMh6vhapqm/sxy/2nDF3Jr1GDiJQDHwB+7H4eytJnY4Fh//0wUeonWjmlGTH6jjWmquoJcG7WQJnbHuuajKlr5YZdLsTxEMbltXBDVm8DtTg3jv0kWfILCC/5Ndqvxb8CXwaC7uekS58xtq4DOH+Y/FZE3hKnFBuMwO+HVQ3tJ5lySuONVEs7jTpEpAD4JfB5VW2R2FtCj+lroc6zfReIyCTgv4BF0bq5r2PyWojIB4FaVX1LRK4LNUfpOqavQxhXqupxESkDXhaRPXH6Dtm1ME+pn/FcluiU62rjvoZq2I/pEk4ikoUjSD9T1V+5zePyWoRQ1SbgVZx1gaEq+TVauBL4sIgcxAnfr8DxnMbbdQBAVY+7r7U4f6gsZwR+P0yU+hnPZYnCSzhFlnb6tJtZcxnQ7LrsLwE3ichkN/vmJrdt1ODG/h8Hdqvqd8MOjcdrUep6SIhIHnAjzhrbUJX8GhWo6gOqWq5OYdHVON/rk4yz6wAgIhNEpDD0Hufnegcj8fuR7gyPTPqHk0HyeZ1UHwAAA3NJREFUHk48/Svpns8wfcdfACeAXpy/Yu7CiYO/AuxzX6e4fQVno8T9wLvAsjA7f4GzgFsNfCbd32sQ1+EqnDDCduBt99/7x+m1OA+npNd298bzdbd9Ds7NtBr4TyDHbc91P1e7x+eE2fqKe432Arek+7udxTW5jv7su3F3Hdzv/I77b2fofjgSvx9WZsgwDMPIGCx8ZxiGYWQMJkqGYRhGxmCiZBiGYWQMJkqGYRhGxmCiZBiGYWQMJkqGYRhGxmCiZBgjgIjcKiL/JiLPishN42Vsw0gVEyXDGAFU9b9V9bPA/wA+AU4hWBHpdAuh4rb9lYioiFwb1nav23ZjLPsiMk9E3o1oyxGRA8B7UcbOc/fJ6RGRkiH8qoZxVpgoGcbI8lWcJ99D7FfVC8I+n4dTWWERgIjk41TdqMN5Uj4WNcBMEQn/nb4beE1Vd0WOraqd7rijriabMbYxUTKMIUZEqkTkVRHZIiLfcXfjFBFZA7yg7j5OMTgXpxTUQvfz3+KUsgmq6inXfqUbitsizo6xC1Q1CBwGKtw+ecAXgW+kMLZhpB0TJcMYQkTEC/wU+IKqLgPycGqH/Q1OodPbROT/iWNiEc5OpwtFpAgn3PYnnJp0ocrmPw6z/w2cHUDBKaIaErN7gHWqejCFsQ0j7dh+SoYxtNwK7ArzSHbjbBL3feD78U4UkZlAvarWuHvYfBn4P8B8nJBeyP4S4Jfu3k8+4A9hYy0Qkd/jiNJlAMmMbRiZgomSYQwtF+JUHA9xPs5OrslwHv3rRq0420Yvx9nTJyRy5+NUbH48yvm7cfYA+hzOHlGnUpu6YaQfC98ZxtBSjxtCE5FLgU/T7+Uk4lz6RemfgXvV2RH23DAbJ4CbQwkNInKu9G+XuxtHxP7CPd8wRh0mSoYxtDwJLHPTsz+CI1LVSZ57Lu7akao+p6pvuO2LgVAG3RM4v7e73VTy+7R//5m9ro3HVLX5rL+JYaQBC98ZxhCiqqeBS6Fvjeg6NzMumXM/GaO9LOx9J/27oEb268Z+p41RjnlKhjF8nE/80F0AKAp/eHakCD08C2QBSYmmYYwEtvOsYRiGkTGYp2QYhmFkDCZKhmEYRsZgomQYhmFkDCZKhmEYRsZgomQYhmFkDCZKhmEYRsZgomQYhmFkDCZKhmEYRsZgomQYhmFkDP8/qxvmJg8u0UsAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.clf()\n",
"plt.plot(test_q, calcs_test)\n",
"# plt.axvline(x=jpsi_mass -70,color='red', linewidth=1.0, linestyle = 'dotted')\n",
"# plt.axvline(x=jpsi_mass +70,color='red', linewidth=1.0, linestyle = 'dotted')\n",
"# plt.axvline(x=psi2s_mass -50,color='red', linewidth=1.0, linestyle = 'dotted')\n",
"# plt.axvline(x=psi2s_mass +50,color='red', linewidth=1.0, linestyle = 'dotted')\n",
"plt.ylim(0.0, 1.5e-6)\n",
"plt.xlabel(r'$q^2 [MeV^2]$')\n",
"plt.savefig('test.png')"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"45702000\n",
"0 22\n",
"1 22\n",
"2 22\n",
"3 22\n",
"4 22\n",
"5 22\n",
"6 22\n",
"7 22\n",
"8 22\n",
"9 22\n",
"10 22\n",
"11 22\n",
"12 22\n",
"13 22\n",
"14 22\n",
"15 22\n",
"16 22\n",
"17 22\n",
"18 22\n",
"19 22\n",
"20 22\n",
"21 22\n",
"Full integration finished in 4 min, 19 s\n",
"1051146\n"
]
}
],
"source": [
"# total_f_fit.update_integration_options(draws_per_dim=2000000, mc_sampler=None)\n",
"\n",
"start = time.time()\n",
"\n",
"_max_size = 2000000\n",
"\n",
"step_size = 1000\n",
"\n",
"steps = np.arange(x_min, x_max, 0.1/step_size)\n",
"\n",
"l = len(steps)\n",
"\n",
"parts = int(l/_max_size)\n",
"\n",
"print(l)\n",
"\n",
"start_ = time.time()\n",
"\n",
"_list = []\n",
"\n",
"for j in range(parts):\n",
" \n",
" print(j, parts)\n",
" \n",
" _c = total_f_fit.pdf(steps[j*_max_size:(j+1)*_max_size], norm_range=False)\n",
"\n",
" inte_fl = zfit.run(_c)\n",
"\n",
" for i in range(int(l/step_size)):\n",
" _list.append(np.mean(inte_fl[int(i*step_size):int((i+1)*step_size)]))\n",
" \n",
"_c = total_f_fit.pdf(steps[(parts-1)*_max_size:], norm_range=False)\n",
"\n",
"inte_fl = zfit.run(_c)\n",
"\n",
"for i in range(int(l/step_size)):\n",
" _list.append(np.mean(steps[int(i*step_size):int((i+1)*step_size)]))\n",
" \n",
"print('Full integration finished in {}'.format(display_time(int(time.time()-start))))\n",
"print(len(_list))"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [],
"source": [
"dirName = 'data/CLs'\n",
"with open(\"{}/inte_100keV_steps.pkl\".format(dirName), \"wb\") as f:\n",
" pkl.dump(_list, f, pkl.HIGHEST_PROTOCOL)"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"45701.99999999999"
]
},
"execution_count": 40,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"(x_max-x_min)/0.1"
]
},
{
"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
}
