{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Import"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\util\\execution.py:57: UserWarning: Not running on Linux. Determining available cpus for thread can failand be overestimated. Workaround (only if too many cpus are used):`zfit.run.set_n_cpu(your_cpu_number)`\n",
" warnings.warn(\"Not running on Linux. Determining available cpus for thread can fail\"\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"WARNING: The TensorFlow contrib module will not be included in TensorFlow 2.0.\n",
"For more information, please see:\n",
" * https://github.com/tensorflow/community/blob/master/rfcs/20180907-contrib-sunset.md\n",
" * https://github.com/tensorflow/addons\n",
"If you depend on functionality not listed there, please file an issue.\n",
"\n"
]
}
],
"source": [
"import os\n",
"\n",
"# os.environ[\"CUDA_VISIBLE_DEVICES\"] = \"-1\"\n",
"\n",
"import numpy as np\n",
"from pdg_const import pdg\n",
"import matplotlib\n",
"import matplotlib.pyplot as plt\n",
"import pickle as pkl\n",
"import sys\n",
"import time\n",
"from helperfunctions import display_time, prepare_plot\n",
"import cmath as c\n",
"import scipy.integrate as integrate\n",
"from scipy.optimize import fminbound\n",
"from array import array as arr\n",
"import collections\n",
"from itertools import compress\n",
"import tensorflow as tf\n",
"import zfit\n",
"from zfit import ztf\n",
"# from IPython.display import clear_output\n",
"import os\n",
"import tensorflow_probability as tfp\n",
"tfd = tfp.distributions"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# chunksize = 10000\n",
"# zfit.run.chunking.active = True\n",
"# zfit.run.chunking.max_n_points = chunksize"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Build model and graphs\n",
"## Create graphs"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"def formfactor(q2, subscript, b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2): #returns real value\n",
" #check if subscript is viable\n",
"\n",
" if subscript != \"0\" and subscript != \"+\" and subscript != \"T\":\n",
" raise ValueError('Wrong subscript entered, choose either 0, + or T')\n",
"\n",
" #get constants\n",
"\n",
" mK = ztf.constant(pdg['Ks_M'])\n",
" mbstar0 = ztf.constant(pdg[\"mbstar0\"])\n",
" mbstar = ztf.constant(pdg[\"mbstar\"])\n",
"\n",
"\n",
" mmu = ztf.constant(pdg['muon_M'])\n",
" mb = ztf.constant(pdg['bquark_M'])\n",
" ms = ztf.constant(pdg['squark_M'])\n",
" mB = ztf.constant(pdg['Bplus_M'])\n",
"\n",
" #N comes from derivation in paper\n",
"\n",
" N = 3\n",
"\n",
" #some helperfunctions\n",
"\n",
" tpos = (mB - mK)**2\n",
" tzero = (mB + mK)*(ztf.sqrt(mB)-ztf.sqrt(mK))**2\n",
"\n",
" z_oben = ztf.sqrt(tpos - q2) - ztf.sqrt(tpos - tzero)\n",
" z_unten = ztf.sqrt(tpos - q2) + ztf.sqrt(tpos - tzero)\n",
" z = tf.divide(z_oben, z_unten)\n",
"\n",
" #calculate f0\n",
"\n",
" if subscript == \"0\":\n",
" prefactor = 1/(1 - q2/(mbstar0**2))\n",
" _sum = 0\n",
" b0 = [b0_0, b0_1, b0_2]\n",
"\n",
" for i in range(N):\n",
" _sum += b0[i]*(tf.pow(z,i))\n",
"\n",
" return ztf.to_complex(prefactor * _sum)\n",
"\n",
" #calculate f+ or fT\n",
"\n",
" else:\n",
" prefactor = 1/(1 - q2/(mbstar**2))\n",
" _sum = 0\n",
"\n",
" if subscript == \"T\":\n",
" bT = [bT_0, bT_1, bT_2]\n",
" for i in range(N):\n",
" _sum += bT[i] * (tf.pow(z, i) - ((-1)**(i-N)) * (i/N) * tf.pow(z, N))\n",
" else:\n",
" bplus = [bplus_0, bplus_1, bplus_2]\n",
" for i in range(N):\n",
" _sum += bplus[i] * (tf.pow(z, i) - ((-1)**(i-N)) * (i/N) * tf.pow(z, N))\n",
"\n",
" return ztf.to_complex(prefactor * _sum)\n",
"\n",
"def resonance(q, _mass, width, phase, scale):\n",
"\n",
" q2 = tf.pow(q, 2)\n",
"\n",
" mmu = ztf.constant(pdg['muon_M'])\n",
"\n",
" p = 0.5 * ztf.sqrt(q2 - 4*(mmu**2))\n",
"\n",
" p0 = 0.5 * ztf.sqrt(_mass**2 - 4*mmu**2)\n",
"\n",
" gamma_j = tf.divide(p, q) * _mass * width / p0\n",
"\n",
" #Calculate the resonance\n",
"\n",
" _top = tf.complex(_mass * width, ztf.constant(0.0))\n",
"\n",
" _bottom = tf.complex(_mass**2 - q2, -_mass*gamma_j)\n",
"\n",
" com = _top/_bottom\n",
"\n",
" #Rotate by the phase\n",
"\n",
" r = ztf.to_complex(scale*tf.abs(com))\n",
"\n",
" _phase = tf.angle(com)\n",
"\n",
" _phase += phase\n",
"\n",
" com = r * tf.exp(tf.complex(ztf.constant(0.0), _phase))\n",
"\n",
" return com\n",
"\n",
"\n",
"def axiv_nonres(q, b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2):\n",
"\n",
" GF = ztf.constant(pdg['GF'])\n",
" alpha_ew = ztf.constant(pdg['alpha_ew'])\n",
" Vtb = ztf.constant(pdg['Vtb'])\n",
" Vts = ztf.constant(pdg['Vts'])\n",
" C10eff = ztf.constant(pdg['C10eff'])\n",
"\n",
" mmu = ztf.constant(pdg['muon_M'])\n",
" mb = ztf.constant(pdg['bquark_M'])\n",
" ms = ztf.constant(pdg['squark_M'])\n",
" mK = ztf.constant(pdg['Ks_M'])\n",
" mB = ztf.constant(pdg['Bplus_M'])\n",
"\n",
" q2 = tf.pow(q, 2)\n",
"\n",
" #Some helperfunctions\n",
"\n",
" beta = 1. - 4. * mmu**2. / q2\n",
"\n",
" kabs = ztf.sqrt(mB**2. + tf.pow(q2, 2)/mB**2. + mK**4./mB**2. - 2. * (mB**2. * mK**2. + mK**2. * q2 + mB**2. * q2) / mB**2.)\n",
"\n",
" #prefactor in front of whole bracket\n",
"\n",
" prefactor1 = GF**2. *alpha_ew**2. * (tf.abs(Vtb*Vts))**2. * kabs * beta / (128. * np.pi**5.)\n",
"\n",
" #left term in bracket\n",
"\n",
" bracket_left = 2./3. * tf.pow(kabs,2) * tf.pow(beta,2) * tf.pow(tf.abs(ztf.to_complex(C10eff)*formfactor(q2, \"+\", b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2)),2)\n",
"\n",
" #middle term in bracket\n",
"\n",
" _top = 4. * mmu**2. * (mB**2. - mK**2.) * (mB**2. - mK**2.)\n",
"\n",
" _under = q2 * mB**2.\n",
"\n",
" bracket_middle = _top/_under *tf.pow(tf.abs(ztf.to_complex(C10eff) * formfactor(q2, \"0\", b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2)), 2)\n",
" \n",
" #Note sqrt(q2) comes from derivation as we use q2 and plot q\n",
"\n",
" return prefactor1 * (bracket_left + bracket_middle) * 2 * q\n",
"\n",
"def vec(q, funcs, b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2):\n",
" \n",
" q2 = tf.pow(q, 2)\n",
"\n",
" GF = ztf.constant(pdg['GF'])\n",
" alpha_ew = ztf.constant(pdg['alpha_ew'])\n",
" Vtb = ztf.constant(pdg['Vtb'])\n",
" Vts = ztf.constant(pdg['Vts'])\n",
" C7eff = ztf.constant(pdg['C7eff'])\n",
"\n",
" mmu = ztf.constant(pdg['muon_M'])\n",
" mb = ztf.constant(pdg['bquark_M'])\n",
" ms = ztf.constant(pdg['squark_M'])\n",
" mK = ztf.constant(pdg['Ks_M'])\n",
" mB = ztf.constant(pdg['Bplus_M'])\n",
"\n",
" #Some helperfunctions\n",
"\n",
" beta = 1. - 4. * mmu**2. / q2\n",
"\n",
" kabs = ztf.sqrt(mB**2. + tf.pow(q2, 2)/mB**2. + mK**4./mB**2. - 2 * (mB**2 * mK**2 + mK**2 * q2 + mB**2 * q2) / mB**2)\n",
" \n",
" #prefactor in front of whole bracket\n",
"\n",
" prefactor1 = GF**2. *alpha_ew**2. * (tf.abs(Vtb*Vts))**2 * kabs * beta / (128. * np.pi**5.)\n",
"\n",
" #right term in bracket\n",
"\n",
" prefactor2 = tf.pow(kabs,2) * (1. - 1./3. * beta)\n",
"\n",
" abs_bracket = tf.pow(tf.abs(c9eff(q, funcs) * formfactor(q2, \"+\", b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2) + ztf.to_complex(2.0 * C7eff * (mb + ms)/(mB + mK)) * formfactor(q2, \"T\", b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2)),2)\n",
"\n",
" bracket_right = prefactor2 * abs_bracket\n",
"\n",
" #Note sqrt(q2) comes from derivation as we use q2 and plot q\n",
"\n",
" return prefactor1 * bracket_right * 2 * q\n",
"\n",
"def c9eff(q, funcs):\n",
"\n",
" C9eff_nr = ztf.to_complex(ztf.constant(pdg['C9eff']))\n",
"\n",
" c9 = C9eff_nr + funcs\n",
"\n",
" return c9"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def G(y):\n",
" \n",
" def inner_rect_bracket(q):\n",
" return tf.log(ztf.to_complex((1+tf.sqrt(q))/(1-tf.sqrt(q)))-tf.complex(ztf.constant(0), -1*ztf.constant(np.pi))) \n",
" \n",
" def inner_right(q):\n",
" return ztf.to_complex(2 * tf.atan(1/tf.sqrt(tf.math.real(-q))))\n",
" \n",
" big_bracket = tf.where(tf.math.real(y) > ztf.constant(0.0), inner_rect_bracket(y), inner_right(y))\n",
" \n",
" return ztf.to_complex(tf.sqrt(tf.abs(y))) * big_bracket\n",
"\n",
"def h_S(m, q):\n",
" \n",
" return ztf.to_complex(2) - G(ztf.to_complex(1) - ztf.to_complex(4*tf.pow(m, 2)) / ztf.to_complex(tf.pow(q, 2)))\n",
"\n",
"def h_P(m, q):\n",
" \n",
" return ztf.to_complex(2/3) + (ztf.to_complex(1) - ztf.to_complex(4*tf.pow(m, 2)) / ztf.to_complex(tf.pow(q, 2))) * h_S(m,q)\n",
"\n",
"def two_p_ccbar(mD, m_D_bar, m_D_star, q):\n",
" \n",
" \n",
" #Load constants\n",
" nu_D_bar = ztf.to_complex(pdg[\"nu_D_bar\"])\n",
" nu_D = ztf.to_complex(pdg[\"nu_D\"])\n",
" nu_D_star = ztf.to_complex(pdg[\"nu_D_star\"])\n",
" \n",
" phase_D_bar = ztf.to_complex(pdg[\"phase_D_bar\"])\n",
" phase_D = ztf.to_complex(pdg[\"phase_D\"])\n",
" phase_D_star = ztf.to_complex(pdg[\"phase_D_star\"])\n",
" \n",
" #Calculation\n",
" left_part = nu_D_bar * tf.exp(tf.complex(ztf.constant(0.0), phase_D_bar)) * h_S(m_D_bar, q) \n",
" \n",
" right_part_D = nu_D * tf.exp(tf.complex(ztf.constant(0.0), phase_D)) * h_P(m_D, q) \n",
" \n",
" right_part_D_star = nu_D_star * tf.exp(tf.complex(ztf.constant(0.0), phase_D_star)) * h_P(m_D_star, q) \n",
"\n",
" return left_part + right_part_D + right_part_D_star"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Build pdf"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"class total_pdf_cut(zfit.pdf.ZPDF):\n",
" _N_OBS = 1 # dimension, can be omitted\n",
" _PARAMS = ['b0_0', 'b0_1', 'b0_2', \n",
" 'bplus_0', 'bplus_1', 'bplus_2', \n",
" 'bT_0', 'bT_1', 'bT_2', \n",
" 'rho_mass', 'rho_scale', 'rho_phase', 'rho_width',\n",
" 'jpsi_mass', 'jpsi_scale', 'jpsi_phase', 'jpsi_width',\n",
" 'psi2s_mass', 'psi2s_scale', 'psi2s_phase', 'psi2s_width',\n",
" 'p3770_mass', 'p3770_scale', 'p3770_phase', 'p3770_width',\n",
" 'p4040_mass', 'p4040_scale', 'p4040_phase', 'p4040_width',\n",
" 'p4160_mass', 'p4160_scale', 'p4160_phase', 'p4160_width',\n",
" 'p4415_mass', 'p4415_scale', 'p4415_phase', 'p4415_width',\n",
" 'omega_mass', 'omega_scale', 'omega_phase', 'omega_width',\n",
" 'phi_mass', 'phi_scale', 'phi_phase', 'phi_width',\n",
" 'Dbar_mass', 'Dbar_scale', 'Dbar_phase',\n",
" 'Dstar_mass', 'DDstar_scale', 'DDstar_phase', 'D_mass',\n",
" 'tau_mass', 'C_tt']\n",
"# the name of the parameters\n",
"\n",
" def _unnormalized_pdf(self, x):\n",
" \n",
" x = x.unstack_x()\n",
" \n",
" b0 = [self.params['b0_0'], self.params['b0_1'], self.params['b0_2']]\n",
" bplus = [self.params['bplus_0'], self.params['bplus_1'], self.params['bplus_2']]\n",
" bT = [self.params['bT_0'], self.params['bT_1'], self.params['bT_2']]\n",
" \n",
" def rho_res(q):\n",
" return resonance(q, _mass = self.params['rho_mass'], scale = self.params['rho_scale'],\n",
" phase = self.params['rho_phase'], width = self.params['rho_width'])\n",
" \n",
" def omega_res(q):\n",
" return resonance(q, _mass = self.params['omega_mass'], scale = self.params['omega_scale'],\n",
" phase = self.params['omega_phase'], width = self.params['omega_width'])\n",
" \n",
" def phi_res(q):\n",
" return resonance(q, _mass = self.params['phi_mass'], scale = self.params['phi_scale'],\n",
" phase = self.params['phi_phase'], width = self.params['phi_width'])\n",
"\n",
" def jpsi_res(q):\n",
" return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['jpsi_mass'], 2)) * resonance(q, _mass = self.params['jpsi_mass'], \n",
" scale = self.params['jpsi_scale'],\n",
" phase = self.params['jpsi_phase'], \n",
" width = self.params['jpsi_width'])\n",
" def psi2s_res(q):\n",
" return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['psi2s_mass'], 2)) * resonance(q, _mass = self.params['psi2s_mass'], \n",
" scale = self.params['psi2s_scale'],\n",
" phase = self.params['psi2s_phase'], \n",
" width = self.params['psi2s_width'])\n",
" def p3770_res(q):\n",
" return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['p3770_mass'], 2)) * resonance(q, _mass = self.params['p3770_mass'], \n",
" scale = self.params['p3770_scale'],\n",
" phase = self.params['p3770_phase'], \n",
" width = self.params['p3770_width'])\n",
" \n",
" def p4040_res(q):\n",
" return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['p4040_mass'], 2)) * resonance(q, _mass = self.params['p4040_mass'], \n",
" scale = self.params['p4040_scale'],\n",
" phase = self.params['p4040_phase'], \n",
" width = self.params['p4040_width'])\n",
" \n",
" def p4160_res(q):\n",
" return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['p4160_mass'], 2)) * resonance(q, _mass = self.params['p4160_mass'], \n",
" scale = self.params['p4160_scale'],\n",
" phase = self.params['p4160_phase'], \n",
" width = self.params['p4160_width'])\n",
" \n",
" def p4415_res(q):\n",
" return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['p4415_mass'], 2)) * resonance(q, _mass = self.params['p4415_mass'], \n",
" scale = self.params['p4415_scale'],\n",
" phase = self.params['p4415_phase'], \n",
" width = self.params['p4415_width'])\n",
" \n",
" def P2_D(q):\n",
" Dbar_contrib = ztf.to_complex(self.params['Dbar_scale'])*tf.exp(tf.complex(ztf.constant(0.0), self.params['Dbar_phase']))*ztf.to_complex(h_S(self.params['Dbar_mass'], q))\n",
" DDstar_contrib = ztf.to_complex(self.params['DDstar_scale'])*tf.exp(tf.complex(ztf.constant(0.0), self.params['DDstar_phase']))*(ztf.to_complex(h_P(self.params['Dstar_mass'], q)) + ztf.to_complex(h_P(self.params['D_mass'], q)))\n",
" return Dbar_contrib + DDstar_contrib\n",
" \n",
" def ttau_cusp(q):\n",
" return ztf.to_complex(self.params['C_tt'])*(ztf.to_complex((h_S(self.params['tau_mass'], q))) - ztf.to_complex(h_P(self.params['tau_mass'], q)))\n",
" \n",
"\n",
" funcs = rho_res(x) + omega_res(x) + phi_res(x) + jpsi_res(x) + psi2s_res(x) + p3770_res(x) + p4040_res(x)+ p4160_res(x) + p4415_res(x) + P2_D(x) + ttau_cusp(x)\n",
"\n",
" vec_f = vec(x, funcs, b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2)\n",
"\n",
" axiv_nr = axiv_nonres(x, b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2)\n",
"\n",
" tot = vec_f + axiv_nr\n",
" \n",
" #Cut out jpsi and psi2s\n",
" \n",
" tot = tf.where(tf.math.logical_or(x < ztf.constant(jpsi_mass-60.), x > ztf.constant(jpsi_mass+70.)), tot, 0.0*tot)\n",
" \n",
" tot = tf.where(tf.math.logical_or(x < ztf.constant(psi2s_mass-50.), x > ztf.constant(psi2s_mass+50.)), tot, 0.0*tot)\n",
" \n",
" return tot\n",
" \n",
"class total_pdf_full(zfit.pdf.ZPDF):\n",
" _N_OBS = 1 # dimension, can be omitted\n",
" _PARAMS = ['b0_0', 'b0_1', 'b0_2', \n",
" 'bplus_0', 'bplus_1', 'bplus_2', \n",
" 'bT_0', 'bT_1', 'bT_2', \n",
" 'rho_mass', 'rho_scale', 'rho_phase', 'rho_width',\n",
" 'jpsi_mass', 'jpsi_scale', 'jpsi_phase', 'jpsi_width',\n",
" 'psi2s_mass', 'psi2s_scale', 'psi2s_phase', 'psi2s_width',\n",
" 'p3770_mass', 'p3770_scale', 'p3770_phase', 'p3770_width',\n",
" 'p4040_mass', 'p4040_scale', 'p4040_phase', 'p4040_width',\n",
" 'p4160_mass', 'p4160_scale', 'p4160_phase', 'p4160_width',\n",
" 'p4415_mass', 'p4415_scale', 'p4415_phase', 'p4415_width',\n",
" 'omega_mass', 'omega_scale', 'omega_phase', 'omega_width',\n",
" 'phi_mass', 'phi_scale', 'phi_phase', 'phi_width',\n",
" 'Dbar_mass', 'Dbar_scale', 'Dbar_phase',\n",
" 'Dstar_mass', 'DDstar_scale', 'DDstar_phase', 'D_mass',\n",
" 'tau_mass', 'C_tt']\n",
"# the name of the parameters\n",
"\n",
" def _unnormalized_pdf(self, x):\n",
" \n",
" x = x.unstack_x()\n",
" \n",
" b0 = [self.params['b0_0'], self.params['b0_1'], self.params['b0_2']]\n",
" bplus = [self.params['bplus_0'], self.params['bplus_1'], self.params['bplus_2']]\n",
" bT = [self.params['bT_0'], self.params['bT_1'], self.params['bT_2']]\n",
" \n",
" def rho_res(q):\n",
" return resonance(q, _mass = self.params['rho_mass'], scale = self.params['rho_scale'],\n",
" phase = self.params['rho_phase'], width = self.params['rho_width'])\n",
" \n",
" def omega_res(q):\n",
" return resonance(q, _mass = self.params['omega_mass'], scale = self.params['omega_scale'],\n",
" phase = self.params['omega_phase'], width = self.params['omega_width'])\n",
" \n",
" def phi_res(q):\n",
" return resonance(q, _mass = self.params['phi_mass'], scale = self.params['phi_scale'],\n",
" phase = self.params['phi_phase'], width = self.params['phi_width'])\n",
"\n",
" def jpsi_res(q):\n",
" return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['jpsi_mass'], 2)) * resonance(q, _mass = self.params['jpsi_mass'], \n",
" scale = self.params['jpsi_scale'],\n",
" phase = self.params['jpsi_phase'], \n",
" width = self.params['jpsi_width'])\n",
" def psi2s_res(q):\n",
" return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['psi2s_mass'], 2)) * resonance(q, _mass = self.params['psi2s_mass'], \n",
" scale = self.params['psi2s_scale'],\n",
" phase = self.params['psi2s_phase'], \n",
" width = self.params['psi2s_width'])\n",
" def p3770_res(q):\n",
" return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['p3770_mass'], 2)) * resonance(q, _mass = self.params['p3770_mass'], \n",
" scale = self.params['p3770_scale'],\n",
" phase = self.params['p3770_phase'], \n",
" width = self.params['p3770_width'])\n",
" \n",
" def p4040_res(q):\n",
" return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['p4040_mass'], 2)) * resonance(q, _mass = self.params['p4040_mass'], \n",
" scale = self.params['p4040_scale'],\n",
" phase = self.params['p4040_phase'], \n",
" width = self.params['p4040_width'])\n",
" \n",
" def p4160_res(q):\n",
" return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['p4160_mass'], 2)) * resonance(q, _mass = self.params['p4160_mass'], \n",
" scale = self.params['p4160_scale'],\n",
" phase = self.params['p4160_phase'], \n",
" width = self.params['p4160_width'])\n",
" \n",
" def p4415_res(q):\n",
" return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['p4415_mass'], 2)) * resonance(q, _mass = self.params['p4415_mass'], \n",
" scale = self.params['p4415_scale'],\n",
" phase = self.params['p4415_phase'], \n",
" width = self.params['p4415_width'])\n",
" \n",
" def P2_D(q):\n",
" Dbar_contrib = ztf.to_complex(self.params['Dbar_scale'])*tf.exp(tf.complex(ztf.constant(0.0), self.params['Dbar_phase']))*ztf.to_complex(h_S(self.params['Dbar_mass'], q))\n",
" DDstar_contrib = ztf.to_complex(self.params['DDstar_scale'])*tf.exp(tf.complex(ztf.constant(0.0), self.params['DDstar_phase']))*(ztf.to_complex(h_P(self.params['Dstar_mass'], q)) + ztf.to_complex(h_P(self.params['D_mass'], q)))\n",
" return Dbar_contrib + DDstar_contrib\n",
" \n",
" def ttau_cusp(q):\n",
" return ztf.to_complex(self.params['C_tt'])*(ztf.to_complex((h_S(self.params['tau_mass'], q))) - ztf.to_complex(h_P(self.params['tau_mass'], q)))\n",
" \n",
"\n",
" funcs = rho_res(x) + omega_res(x) + phi_res(x) + jpsi_res(x) + psi2s_res(x) + p3770_res(x) + p4040_res(x)+ p4160_res(x) + p4415_res(x) + P2_D(x) + ttau_cusp(x)\n",
"\n",
" vec_f = vec(x, funcs, b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2)\n",
"\n",
" axiv_nr = axiv_nonres(x, b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2)\n",
"\n",
" tot = vec_f + axiv_nr\n",
" \n",
" #Cut out jpsi and psi2s\n",
" \n",
"# tot = tf.where(tf.math.logical_or(x < ztf.constant(jpsi_mass-60.), x > ztf.constant(jpsi_mass+70.)), tot, 0.0*tot)\n",
" \n",
"# tot = tf.where(tf.math.logical_or(x < ztf.constant(psi2s_mass-50.), x > ztf.constant(psi2s_mass+50.)), tot, 0.0*tot)\n",
" \n",
" return tot"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup parameters"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"WARNING:tensorflow:From C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\ops\\resource_variable_ops.py:435: colocate_with (from tensorflow.python.framework.ops) is deprecated and will be removed in a future version.\n",
"Instructions for updating:\n",
"Colocations handled automatically by placer.\n"
]
}
],
"source": [
"# formfactors\n",
"\n",
"b0_0 = zfit.Parameter(\"b0_0\", ztf.constant(0.292), floating = False) #, lower_limit = -2.0, upper_limit= 2.0)\n",
"b0_1 = zfit.Parameter(\"b0_1\", ztf.constant(0.281), floating = False) #, lower_limit = -2.0, upper_limit= 2.0)\n",
"b0_2 = zfit.Parameter(\"b0_2\", ztf.constant(0.150), floating = False) #, lower_limit = -2.0, upper_limit= 2.0)\n",
"\n",
"bplus_0 = zfit.Parameter(\"bplus_0\", ztf.constant(0.466), lower_limit = -2.0, upper_limit= 2.0)\n",
"bplus_1 = zfit.Parameter(\"bplus_1\", ztf.constant(-0.885), lower_limit = -2.0, upper_limit= 2.0)\n",
"bplus_2 = zfit.Parameter(\"bplus_2\", ztf.constant(-0.213), lower_limit = -2.0, upper_limit= 2.0)\n",
"\n",
"bT_0 = zfit.Parameter(\"bT_0\", ztf.constant(0.460), floating = False) #, lower_limit = -2.0, upper_limit= 2.0)\n",
"bT_1 = zfit.Parameter(\"bT_1\", ztf.constant(-1.089), floating = False) #, lower_limit = -2.0, upper_limit= 2.0)\n",
"bT_2 = zfit.Parameter(\"bT_2\", ztf.constant(-1.114), floating = False) #, lower_limit = -2.0, upper_limit= 2.0)\n",
"\n",
"\n",
"#rho\n",
"\n",
"rho_mass, rho_width, rho_phase, rho_scale = pdg[\"rho\"]\n",
"\n",
"rho_m = zfit.Parameter(\"rho_m\", ztf.constant(rho_mass), floating = False) #lower_limit = rho_mass - rho_width, upper_limit = rho_mass + rho_width)\n",
"rho_w = zfit.Parameter(\"rho_w\", ztf.constant(rho_width), floating = False)\n",
"rho_p = zfit.Parameter(\"rho_p\", ztf.constant(rho_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"rho_s = zfit.Parameter(\"rho_s\", ztf.constant(rho_scale), lower_limit=rho_scale-np.sqrt(rho_scale), upper_limit=rho_scale+np.sqrt(rho_scale))\n",
"\n",
"#omega\n",
"\n",
"omega_mass, omega_width, omega_phase, omega_scale = pdg[\"omega\"]\n",
"\n",
"omega_m = zfit.Parameter(\"omega_m\", ztf.constant(omega_mass), floating = False)\n",
"omega_w = zfit.Parameter(\"omega_w\", ztf.constant(omega_width), floating = False)\n",
"omega_p = zfit.Parameter(\"omega_p\", ztf.constant(omega_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"omega_s = zfit.Parameter(\"omega_s\", ztf.constant(omega_scale), lower_limit=omega_scale-np.sqrt(omega_scale), upper_limit=omega_scale+np.sqrt(omega_scale))\n",
"\n",
"\n",
"#phi\n",
"\n",
"phi_mass, phi_width, phi_phase, phi_scale = pdg[\"phi\"]\n",
"\n",
"phi_m = zfit.Parameter(\"phi_m\", ztf.constant(phi_mass), floating = False)\n",
"phi_w = zfit.Parameter(\"phi_w\", ztf.constant(phi_width), floating = False)\n",
"phi_p = zfit.Parameter(\"phi_p\", ztf.constant(phi_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"phi_s = zfit.Parameter(\"phi_s\", ztf.constant(phi_scale), lower_limit=phi_scale-np.sqrt(phi_scale), upper_limit=phi_scale+np.sqrt(phi_scale))\n",
"\n",
"#jpsi\n",
"\n",
"jpsi_mass, jpsi_width, jpsi_phase, jpsi_scale = pdg[\"jpsi\"]\n",
"\n",
"jpsi_m = zfit.Parameter(\"jpsi_m\", ztf.constant(jpsi_mass), floating = False)\n",
"jpsi_w = zfit.Parameter(\"jpsi_w\", ztf.constant(jpsi_width), floating = False)\n",
"jpsi_p = zfit.Parameter(\"jpsi_p\", ztf.constant(jpsi_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"jpsi_s = zfit.Parameter(\"jpsi_s\", ztf.constant(jpsi_scale), floating = False) #, lower_limit=jpsi_scale-np.sqrt(jpsi_scale), upper_limit=jpsi_scale+np.sqrt(jpsi_scale))\n",
"\n",
"#psi2s\n",
"\n",
"psi2s_mass, psi2s_width, psi2s_phase, psi2s_scale = pdg[\"psi2s\"]\n",
"\n",
"psi2s_m = zfit.Parameter(\"psi2s_m\", ztf.constant(psi2s_mass), floating = False)\n",
"psi2s_w = zfit.Parameter(\"psi2s_w\", ztf.constant(psi2s_width), floating = False)\n",
"psi2s_p = zfit.Parameter(\"psi2s_p\", ztf.constant(psi2s_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"psi2s_s = zfit.Parameter(\"psi2s_s\", ztf.constant(psi2s_scale), floating = False) #, lower_limit=psi2s_scale-np.sqrt(psi2s_scale), upper_limit=psi2s_scale+np.sqrt(psi2s_scale))\n",
"\n",
"#psi(3770)\n",
"\n",
"p3770_mass, p3770_width, p3770_phase, p3770_scale = pdg[\"p3770\"]\n",
"\n",
"p3770_m = zfit.Parameter(\"p3770_m\", ztf.constant(p3770_mass), floating = False)\n",
"p3770_w = zfit.Parameter(\"p3770_w\", ztf.constant(p3770_width), floating = False)\n",
"p3770_p = zfit.Parameter(\"p3770_p\", ztf.constant(p3770_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"p3770_s = zfit.Parameter(\"p3770_s\", ztf.constant(p3770_scale), lower_limit=p3770_scale-np.sqrt(p3770_scale), upper_limit=p3770_scale+np.sqrt(p3770_scale))\n",
"\n",
"#psi(4040)\n",
"\n",
"p4040_mass, p4040_width, p4040_phase, p4040_scale = pdg[\"p4040\"]\n",
"\n",
"p4040_m = zfit.Parameter(\"p4040_m\", ztf.constant(p4040_mass), floating = False)\n",
"p4040_w = zfit.Parameter(\"p4040_w\", ztf.constant(p4040_width), floating = False)\n",
"p4040_p = zfit.Parameter(\"p4040_p\", ztf.constant(p4040_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"p4040_s = zfit.Parameter(\"p4040_s\", ztf.constant(p4040_scale), lower_limit=p4040_scale-np.sqrt(p4040_scale), upper_limit=p4040_scale+np.sqrt(p4040_scale))\n",
"\n",
"#psi(4160)\n",
"\n",
"p4160_mass, p4160_width, p4160_phase, p4160_scale = pdg[\"p4160\"]\n",
"\n",
"p4160_m = zfit.Parameter(\"p4160_m\", ztf.constant(p4160_mass), floating = False)\n",
"p4160_w = zfit.Parameter(\"p4160_w\", ztf.constant(p4160_width), floating = False)\n",
"p4160_p = zfit.Parameter(\"p4160_p\", ztf.constant(p4160_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"p4160_s = zfit.Parameter(\"p4160_s\", ztf.constant(p4160_scale), lower_limit=p4160_scale-np.sqrt(p4160_scale), upper_limit=p4160_scale+np.sqrt(p4160_scale))\n",
"\n",
"#psi(4415)\n",
"\n",
"p4415_mass, p4415_width, p4415_phase, p4415_scale = pdg[\"p4415\"]\n",
"\n",
"p4415_m = zfit.Parameter(\"p4415_m\", ztf.constant(p4415_mass), floating = False)\n",
"p4415_w = zfit.Parameter(\"p4415_w\", ztf.constant(p4415_width), floating = False)\n",
"p4415_p = zfit.Parameter(\"p4415_p\", ztf.constant(p4415_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"p4415_s = zfit.Parameter(\"p4415_s\", ztf.constant(p4415_scale), lower_limit=p4415_scale-np.sqrt(p4415_scale), upper_limit=p4415_scale+np.sqrt(p4415_scale))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Dynamic generation of 2 particle contribution"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"m_c = 1300\n",
"\n",
"Dbar_phase = 0.0\n",
"DDstar_phase = 0.0\n",
"Dstar_mass = pdg['Dst_M']\n",
"Dbar_mass = pdg['D0_M']\n",
"D_mass = pdg['D0_M']\n",
"\n",
"Dbar_s = zfit.Parameter(\"Dbar_s\", ztf.constant(0.0), lower_limit=-0.3, upper_limit=0.3)\n",
"Dbar_m = zfit.Parameter(\"Dbar_m\", ztf.constant(Dbar_mass), floating = False)\n",
"Dbar_p = zfit.Parameter(\"Dbar_p\", ztf.constant(Dbar_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)#, floating = False)\n",
"DDstar_s = zfit.Parameter(\"DDstar_s\", ztf.constant(0.0), lower_limit=-0.3, upper_limit=0.3)#, floating = False)\n",
"Dstar_m = zfit.Parameter(\"Dstar_m\", ztf.constant(Dstar_mass), floating = False)\n",
"D_m = zfit.Parameter(\"D_m\", ztf.constant(D_mass), floating = False)\n",
"DDstar_p = zfit.Parameter(\"DDstar_p\", ztf.constant(DDstar_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)#, floating = False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Tau parameters"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"tau_m = zfit.Parameter(\"tau_m\", ztf.constant(pdg['tau_M']), floating = False)\n",
"Ctt = zfit.Parameter(\"Ctt\", ztf.constant(0.0), lower_limit=-1.5, upper_limit=1.5)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Load data"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"x_min = 2*pdg['muon_M']\n",
"x_max = (pdg[\"Bplus_M\"]-pdg[\"Ks_M\"]-0.1)\n",
"\n",
"# # Full spectrum\n",
"\n",
"obs_toy = zfit.Space('q', limits = (x_min, x_max))\n",
"\n",
"# Jpsi and Psi2s cut out\n",
"\n",
"obs1 = zfit.Space('q', limits = (x_min, jpsi_mass - 60.))\n",
"obs2 = zfit.Space('q', limits = (jpsi_mass + 70., psi2s_mass - 50.))\n",
"obs3 = zfit.Space('q', limits = (psi2s_mass + 50., x_max))\n",
"\n",
"obs_fit = obs1 + obs2 + obs3\n",
"\n",
"# with open(r\"./data/slim_points/slim_points_toy_0_range({0}-{1}).pkl\".format(int(x_min), int(x_max)), \"rb\") as input_file:\n",
"# part_set = pkl.load(input_file)\n",
"\n",
"# x_part = part_set['x_part']\n",
"\n",
"# x_part = x_part.astype('float64')\n",
"\n",
"# data = zfit.data.Data.from_numpy(array=x_part, obs=obs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup pdf"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"total_f = total_pdf_cut(obs=obs_toy, jpsi_mass = jpsi_m, jpsi_scale = jpsi_s, jpsi_phase = jpsi_p, jpsi_width = jpsi_w,\n",
" psi2s_mass = psi2s_m, psi2s_scale = psi2s_s, psi2s_phase = psi2s_p, psi2s_width = psi2s_w,\n",
" p3770_mass = p3770_m, p3770_scale = p3770_s, p3770_phase = p3770_p, p3770_width = p3770_w,\n",
" p4040_mass = p4040_m, p4040_scale = p4040_s, p4040_phase = p4040_p, p4040_width = p4040_w,\n",
" p4160_mass = p4160_m, p4160_scale = p4160_s, p4160_phase = p4160_p, p4160_width = p4160_w,\n",
" p4415_mass = p4415_m, p4415_scale = p4415_s, p4415_phase = p4415_p, p4415_width = p4415_w,\n",
" rho_mass = rho_m, rho_scale = rho_s, rho_phase = rho_p, rho_width = rho_w,\n",
" omega_mass = omega_m, omega_scale = omega_s, omega_phase = omega_p, omega_width = omega_w,\n",
" phi_mass = phi_m, phi_scale = phi_s, phi_phase = phi_p, phi_width = phi_w,\n",
" Dstar_mass = Dstar_m, DDstar_scale = DDstar_s, DDstar_phase = DDstar_p, D_mass = D_m,\n",
" Dbar_mass = Dbar_m, Dbar_scale = Dbar_s, Dbar_phase = Dbar_p,\n",
" tau_mass = tau_m, C_tt = Ctt, b0_0 = b0_0, b0_1 = b0_1, b0_2 = b0_2,\n",
" bplus_0 = bplus_0, bplus_1 = bplus_1, bplus_2 = bplus_2,\n",
" bT_0 = bT_0, bT_1 = bT_1, bT_2 = bT_2)\n",
"\n",
"total_f_fit = total_pdf_full(obs=obs_fit, jpsi_mass = jpsi_m, jpsi_scale = jpsi_s, jpsi_phase = jpsi_p, jpsi_width = jpsi_w,\n",
" psi2s_mass = psi2s_m, psi2s_scale = psi2s_s, psi2s_phase = psi2s_p, psi2s_width = psi2s_w,\n",
" p3770_mass = p3770_m, p3770_scale = p3770_s, p3770_phase = p3770_p, p3770_width = p3770_w,\n",
" p4040_mass = p4040_m, p4040_scale = p4040_s, p4040_phase = p4040_p, p4040_width = p4040_w,\n",
" p4160_mass = p4160_m, p4160_scale = p4160_s, p4160_phase = p4160_p, p4160_width = p4160_w,\n",
" p4415_mass = p4415_m, p4415_scale = p4415_s, p4415_phase = p4415_p, p4415_width = p4415_w,\n",
" rho_mass = rho_m, rho_scale = rho_s, rho_phase = rho_p, rho_width = rho_w,\n",
" omega_mass = omega_m, omega_scale = omega_s, omega_phase = omega_p, omega_width = omega_w,\n",
" phi_mass = phi_m, phi_scale = phi_s, phi_phase = phi_p, phi_width = phi_w,\n",
" Dstar_mass = Dstar_m, DDstar_scale = DDstar_s, DDstar_phase = DDstar_p, D_mass = D_m,\n",
" Dbar_mass = Dbar_m, Dbar_scale = Dbar_s, Dbar_phase = Dbar_p,\n",
" tau_mass = tau_m, C_tt = Ctt, b0_0 = b0_0, b0_1 = b0_1, b0_2 = b0_2,\n",
" bplus_0 = bplus_0, bplus_1 = bplus_1, bplus_2 = bplus_2,\n",
" bT_0 = bT_0, bT_1 = bT_1, bT_2 = bT_2)\n",
" \n",
"# print(total_pdf.obs)\n",
"\n",
"# print(calcs_test)\n",
"\n",
"# for param in total_f.get_dependents():\n",
"# print(zfit.run(param))"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<tf.Tensor 'normalization/hook_integrate/hook_numeric_integrate/mul_1:0' shape=() dtype=float64>"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"total_f_fit.normalization(obs_toy)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Test if graphs actually work and compute values"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"# def total_test_tf(xq):\n",
"\n",
"# def jpsi_res(q):\n",
"# return resonance(q, jpsi_m, jpsi_s, jpsi_p, jpsi_w)\n",
"\n",
"# def psi2s_res(q):\n",
"# return resonance(q, psi2s_m, psi2s_s, psi2s_p, psi2s_w)\n",
"\n",
"# def cusp(q):\n",
"# return bifur_gauss(q, cusp_m, sig_L, sig_R, cusp_s)\n",
"\n",
"# funcs = jpsi_res(xq) + psi2s_res(xq) + cusp(xq)\n",
"\n",
"# vec_f = vec(xq, funcs)\n",
"\n",
"# axiv_nr = axiv_nonres(xq)\n",
"\n",
"# tot = vec_f + axiv_nr\n",
" \n",
"# return tot\n",
"\n",
"# def jpsi_res(q):\n",
"# return resonance(q, jpsi_m, jpsi_s, jpsi_p, jpsi_w)\n",
"\n",
"# calcs = zfit.run(total_test_tf(x_part))\n",
"\n",
"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",
"# res_y = zfit.run(jpsi_res(test_q))\n",
"# b0 = [b0_0, b0_1, b0_2]\n",
"# bplus = [bplus_0, bplus_1, bplus_2]\n",
"# bT = [bT_0, bT_1, bT_2]\n",
"# f0_y = zfit.run(tf.math.real(formfactor(test_q,\"0\", b0, bplus, bT)))\n",
"# fplus_y = zfit.run(tf.math.real(formfactor(test_q,\"+\", b0, bplus, bT)))\n",
"# fT_y = zfit.run(tf.math.real(formfactor(test_q,\"T\", b0, bplus, bT)))"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\ipykernel_launcher.py:12: UserWarning: Creating legend with loc=\"best\" can be slow with large amounts of data.\n",
" if sys.path[0] == '':\n",
"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\IPython\\core\\pylabtools.py:128: UserWarning: Creating legend with loc=\"best\" can be slow with large amounts of data.\n",
" fig.canvas.print_figure(bytes_io, **kw)\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAaUAAAD4CAYAAABMtfkzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO29eXzU13X//T4z2tAOWthBAgQYMN4weCExwU7AdRLc1G5I6sRP4tZ50jhtki6xnz5NnjhxWzdp3PaXza7txElaY8d1auK4xo7BSfCCwQbbiMUIxCI27QJJaJmZ8/zx/Y40jGYVkmYknffrxUszd+733DtfpPnMOffcc0VVMQzDMIx0wJPqCRiGYRhGEBMlwzAMI20wUTIMwzDSBhMlwzAMI20wUTIMwzDShoxUTyDdKC0t1YqKilRPwzDSnn2nzpKX7WXmxNzz2tvO9XK0uZOq8gJyMlP3vVcVdp9oY0phDmUF2TS2d3OyrYtFUwvxeiRl8xqrvPnmm42qWnahdkyUwqioqGDHjh2pnoZhpD1X/+NLrJxXyrdvveS89ufePcmf/+db/OIv38dFUwtTNDs41+Pnoq89z1fXLuTzq+by2KuH+frGal76+w8yKS8rZfMaq4jIkaGwY+E7wzAGhT+gET0Oj0jf66lEccZ3p0NwqqmelxEbEyXDMAZFQMETQZQy3LZAijfmB4cPTjE4VysYkN6YKBmGMSgCqkRamgl6T74UeyRBURSc+QQ9OHOU0htbUzIMY1D4A4pXIoTvgp5SysN3DgPCdyPsKfX29lJXV0dXV9eIjjtc5OTkMGPGDDIzM4fFvomSYRiDIhDQmOG7VK/dBLVHJMxTGuF51dXVUVBQQEVFRd9cRiuqSlNTE3V1dVRWVg7LGBa+MwxjUAQ0iqeULokOfeE7h+C8RnpJqauri5KSklEvSOAIfElJybB6fSZKhmEMCr9G9pSCa0ojHSYLp99Tcn563E+7VMxrLAhSkOF+LyZKhmEMikCg3/sIxZsu4Tv3Z7inlOqsQCM2JkqGYQyKgCreCJ8g3rRJCXfGD3pz/eE7E6VIvPzyy3z4wx8GoLu7mxtuuIFLL72UJ554YkTnYYkOhmEMCr9qZE/JbfP5U50S7vwM95T8gdTMZzSxc+dOent72bVr14iPnZCnJCJrRWS/iNSIyN0RXs8WkSfc17eJSEXIa/e47ftFZE08myJS6do44NrMijWGiJSIyBYRaReR70WZ/0YR2Z3YLTEMIx6qimrs8F3KPSXOX1QKenWpnlcqOHz4MAsXLuT2229n6dKl3HLLLXR2dvL888+zcOFCVq5cydNPPw1AfX09t912G7t27eLSSy/l4MGDIzrXuJ6SiHiB7wMfBOqA7SKyUVX3hHS7A2hR1Xkish64H/i4iCwC1gOLgWnAb0RkvntNNJv3Aw+o6gYR+ZFr+4fRxgC6gL8Hlrj/wuf/MaA9qbtiGEZMgutFkcoM9a8pjeiUBhJW0UHSYE3pG7+qZs+JM0Nqc9G0Qr7+kcVx++3fv59HHnmEa6+9ls9+9rN897vf5cEHH2Tz5s3MmzePj3/84wCUl5fz8MMP853vfIdnn312SOeaCIl4SsuBGlU9pKo9wAZgXVifdcBj7uOngOvF+Q1YB2xQ1W5VrQVqXHsRbbrXrHZt4Nq8OdYYqtqhqltxxOk8RCQf+ArwrQTep2EYCRIMjUUWJeenL5BaVeoP34XvU0rVjFLLzJkzufbaawG47bbb2LFjB5WVlVRVVSEi3HbbbSmeoUMia0rTgWMhz+uAFdH6qKpPRNqAErf99bBrp7uPI9ksAVpV1Rehf7QxGmPM/ZvAvwCdsd6giNwJ3Akwa9asWF0NwyCkhE/EMkOe8/qkivCCrOkQvkvEoxkuwlO529ra0jJVPRFPKdKsw/9Xo/UZqvZE59E/IZFLgXmq+stoffqMqD6kqstUdVlZ2QUfB2IYY56+8F2MRIdUh+80LNEhHcJ3qeTo0aO89tprADz++OPccMMN1NbW9q0ZPf7446mcXh+JiFIdMDPk+QzgRLQ+IpIBFAHNMa6N1t4IFLs2wseKNkY0rgauEJHDwFZgvoi8HPOdGoaREMENqBGPrgh6JGlSkDUYthvv+5QuuugiHnvsMZYuXUpzczNf/vKXeeihh7jppptYuXIls2fPTvUUgcTCd9uBKhGpBI7jJC58MqzPRuB24DXgFmCzqqqIbAT+S0S+i5PoUAW8gfPlZYBN95otro0Nrs1nYo0RbdKq+kOcBAncTL1nVXVVAu/XMIw4qOsFxcq+S3WV8PCKDt4+UUrRhFKMx+PhRz/60Xlta9euZd++fQP6rlq1ilWrVo3QzM4nrii56zd3AZsAL/CoqlaLyL3ADlXdCDwC/ExEanC8l/XutdUi8iSwB/ABX1BVP0Akm+6QXwU2iMi3gJ2ubaKN4do6DBQCWSJyM/ChsOxAwzCGEH+fFzLwtXQrM9TvKTnPU+3BGbFJaPOsqj4HPBfW9rWQx13ArVGuvQ+4LxGbbvshnOy88PZYY1TEmf9hIqSLG4YxOGKmhKeoGnc4feE7N5wYXFNKtVimgoqKCnbvHh1bNa3MkGEYSRNewieUdAnfha8peftOnh35uYyl0kbD/V5MlAzDSBp/2Ad+KN40OeSvb59SePhuhAUiJyeHpqamMSFMwfOUcnJyhm0Mq31nGEbSxEwJT5s1pfPXvfrCdyMsljNmzKCuro6GhoYRHXe4CJ48O1yYKBmGkTTBqgiRwnfpcshfICzRIVXhu8zMzGE7pXUsYuE7wzCSJtC3T2nga+lyHHogzFNKVfjOSA4TJcMwkiaRNaV0ESUJ2zyb6nkZsTFRMgwjaYJJDJFESUQQSb1HMnCf0vjePDtaMFEyDCNpYpUZAieEly4p4X2H/LmfdmMhC24sY6JkGEbS9CU6RCky7RFJm5TwoBh5xvHm2dGEiZJhGEkTXJfJ8ET+CPF6JOVrNxplTcnCd+mNiZJhGEkTPMDP643sKnk9knKPJDwlPOjVWfguvTFRMgwjafo9pRiilCaeUn9KuGXfjQZMlAzDSBpfjIKs4FR6SPWHf7TNs6ZJ6Y2JkmEYSZPImlKqU8LDj2wP/kx1AoYRGxMlwzCSJq6n5BF8/vQQJTt5dnRhomQYRtL43USHaGtKGd7U71MK3zxr4bvRgYmSYRhJE/SConlKmR4Pvf7ASE5pAOG174Lhu1RnBRqxMVEyDCNpYp08C66nlPLwnfMzuE8peMyGpYSnNyZKhmEkjS9OSniGx9O3lylVDKwSnh6HDxqxMVEyDCNp4nlKmV6hN8WeUrSKDimelhGHhERJRNaKyH4RqRGRuyO8ni0iT7ivbxORipDX7nHb94vImng2RaTStXHAtZkVawwRKRGRLSLSLiLfC7GTKyK/FpF9IlItIv+U/O0xDCMS8VLCM7yp95T6Ex2cn8HqE+YppTdxRUlEvMD3gRuBRcAnRGRRWLc7gBZVnQc8ANzvXrsIWA8sBtYCPxARbxyb9wMPqGoV0OLajjoG0AX8PfDXEab/HVVdCFwGXCsiN8Z7v4ZhxKfPU4pSZijDk3pPKXzzbDDUmOqsQCM2iXhKy4EaVT2kqj3ABmBdWJ91wGPu46eA68XxmdcBG1S1W1VrgRrXXkSb7jWrXRu4Nm+ONYaqdqjqVhxx6kNVO1V1i/u4B3gLGL6D5Q1jHBFvTSnT68GXJtl3way7/sMHUzsvIzaJiNJ04FjI8zq3LWIfVfUBbUBJjGujtZcAra6N8LGijREXESkGPgK8FOX1O0Vkh4jsaGhoSMSkYYxrgh/so2FNyROWfWeeUnqTiChF+q0L/1+N1meo2hOdxwBEJAN4HPh3VT0UqY+qPqSqy1R1WVlZWTyThjHuiZt9502HfUrOz76KDh7BI1aQNd1JRJTqgJkhz2cAJ6L1cUWgCGiOcW209kag2LURPla0MeLxEHBAVf81gb6GYSRAItl3qfZIwlPCIZiqbqKUziQiStuBKjcrLgsncWFjWJ+NwO3u41uAzer4zhuB9W7mXCVQBbwRzaZ7zRbXBq7NZ+KMERUR+RaOeH0pgfdpGEaC+OJl33nSYU3J+RlMCYf0OFLDiE1GvA6q6hORu4BNgBd4VFWrReReYIeqbgQeAX4mIjU43st699pqEXkS2AP4gC+oqh8gkk13yK8CG1xB2enaJtoYrq3DQCGQJSI3Ax8CzgB/B+wD3nJ/Mb+nqg8nf5sMwwglkYoO6bOm1N/mZAVaokM6E1eUAFT1OeC5sLavhTzuAm6Ncu19wH2J2HTbD+Fk54W3xxqjIsrUI//FGIZxQQRLCEXNvkurig79c8zw9ntKm6pP8bmfvcnbX/sQRbmZKZmjMRCr6GAYRtL4AwFEnOSBSKRF7TtXEz3nhe/615S++8J7ABxr6RzxuRnRMVEyDCNpfAGN6iWBs08p1WGy8H1K4Hh2flcsMzOcFzp7/CM+NyM6JkqGYSSNP6BR15MgPSo6aF+iQ3+b1yMDkjTO9ZoopRMmSoZhJI3jKUX/+MjMSP2akhJtTSngtjttPT5LfEgnTJQMw0iaeJ5SpusppfLsovDNs3C+pxRs7/aZp5ROmCgZhpE0vkAg5ppShtf5aEnlnqDIm2f7s++ComSeUnphomQYRtLEXVPypr7OXOTNs56+ta5gc7eJUlphomQYRtL4/PHCd85HSyoz8KJtnu1fU3LDd5bokFaYKBmGkTQJe0opzMALBCInOvStKbmffj1W4SGtMFEyDCNp4u1TCq4p9aYwAy9SokOkNaXuXhOldMJEyYjKE9uPUnH3r+my8IYRhl/jZ99Bij2l4ObZkE+50Oy74Ou2ppRemCgZUfnX3xwAoLmjJ8UzMdINnz8Qc59Sn6eU0jUl52eodGZ4PH2eUjDhwcJ36YWJkmEYSdPr174yPZHIdNeUUlnVwa8Dj9cI9ZSCR2tYSnh6YaJkGEbS9PoDZHpjVHRwX0tlVQd/WEIDnJ99FxRMC9+lFyZKRlzsSDQjnHiilJEGa0p9Zz6FV3TwB8N3jhhZRYf0wkTJMIyk6fUrWTFEKSvDeS2V6zWRDiIMTQnvtfBdWmKiZBhG0jieUvQ1paAopTLd2h9QPHJ+RYeIiQ4mSmmFiZIRFTu214hGjy92+C47wwukNjQWKW09wyN961x9iQ6WfZdWmCgZUbG1JCMa8daUsoPhuxR6IYGAnrdxFpxQnr8vFdxNdLDNs2mFiZJhGEnT69eY4bucTDd8l0JR8keoOhG6phTMwjNPKb1ISJREZK2I7BeRGhG5O8Lr2SLyhPv6NhGpCHntHrd9v4isiWdTRCpdGwdcm1mxxhCREhHZIiLtIvK9sHldISLvutf8u4hYRCoJ7GYZ0fDF8ZSyvMHwXeo+8H0BxeMZ6Cn171OyNaV0JK4oiYgX+D5wI7AI+ISILArrdgfQoqrzgAeA+91rFwHrgcXAWuAHIuKNY/N+4AFVrQJaXNtRxwC6gL8H/jrC9H8I3AlUuf/Wxnu/hmHEp8evZGbECN/1eUqpW1MKRFhTyvR66HVFKChOlhKeXiTiKS0HalT1kKr2ABuAdWF91gGPuY+fAq53vZJ1wAZV7VbVWqDGtRfRpnvNatcGrs2bY42hqh2quhVHnPoQkalAoaq+pk4N+5+G2DIM4wLo9QdipoSnw5qSP6Dn7VECJyuw2x8UJUsJT0cSEaXpwLGQ53VuW8Q+quoD2oCSGNdGay8BWl0b4WNFGyPWvOvizBsAEblTRHaIyI6GhoYYJg3DgCRSwlMtSmGeUpbXQ48vgKr2eUomSulFIqIU6TcvPDErWp+hak90HonMaWCj6kOqukxVl5WVlcUwOT4JHpZmGEF6/YG+oquRCHpRqd6nFEmUwBHL4K+1lRlKLxIRpTpgZsjzGcCJaH1EJAMoAppjXButvREodm2EjxVtjFjznhFn3kYMLC/EiISqutl3sauEZ3iEHn9q9ymFp4QHPbhzPf3zMk8pvUhElLYDVW5WXBZO4sLGsD4bgdvdx7cAm911nI3AejdzrhIn2eCNaDbda7a4NnBtPhNnjIio6kngrIhc5a5VfTrElpEA5iEZkQiGvbJihO/AWVdKpacUCGjfCbhBgqLU6Z4RJkLfGpORHmTE66CqPhG5C9gEeIFHVbVaRO4FdqjqRuAR4GciUoPjvax3r60WkSeBPYAP+IKq+gEi2XSH/CqwQUS+Bex0bRNtDNfWYaAQyBKRm4EPqeoe4PPAT4AJwP+6/wzDuACCNeNieUrgJhWkOCU8UqIDQGe3s2ydl5VBe7cPVbXIQJoQV5QAVPU54Lmwtq+FPO4Cbo1y7X3AfYnYdNsP4WTnhbfHGqMiSvsOYEmk14z42B+pEYlen+MpxROl7AxvylPCw/cpBdeUOt3wXW6Wl/ZuHz3+QF9pJCO1WEUHIyoWvjMiEayAEGufEjh7ldIxJRz6RSkv2/lePhLz3H/qLHf/9zvnrWcZAzFRMgwjKfrCd57YnnSWN7XhO3+AAdl32X2i5ITvJmQ63tFIiNK/vLCfDduP8bsDtu0kFiZKRlQsfGdEItE1pezMVItSYGBK+ABPaeTKIZ1sc/b3Hzh9dtjHGs2YKBmGkRTBc4jihu8yvKkN3ykR1pQcEQp6SrlZIxe+O33GEaXjreeGfazRjImSYRhJEfSU4qWEO+G7FCY6RKgSHs1TGu5K4apKU0cPAMdbu+L0Ht+YKBlxsXwHI5SgV5HhSe/wnS8QiJvoEPSUhns/VXu3r+/E26b27mEda7RjomQYRlIEhSYnM3YKdXaGh67eVHpKEK6b/SnhwX1KQU9peOd5psvX97ipvWdYxxrtmCgZhpEUwZBc8CC/aEzI9HIuhaIU6Tj0cE9pQtBTGmaPrq2zF4DK0jyaO3psu0UMTJSMqNgfjhGJLjfUFW+zaW52Rkr35DgFWc//iMsOq30X9JSGXZTO9YtSjz9Ae7cvzhXjFxMlIy6mTUYoiXpKuZnePo8kFTibZ89vC3pKHcE1pRHaPBsqSmAhvFiYKBlRCe5T0pgnhBjjjYQ9pSxHlAKB1Pz+xDq6or/23chsnj3T5YjSnDJXlDos2SEaJkqGYSRFwp6S64V0pSgt3BcIDNjgO8EVobNu4sGEkRKloKdUYp5SPEyUjLhY+M4IpTsJTwlIWQiv168DDiLMzvDg9Qit5xxRKMgZoUSHc72IwKySXACaO0yUomGiZETFEh2MSAQ9n+wEsu8AOrtTI0o9voFHtosIuZleWtxsuKIJmQDDnrredq6XwpxMSvOzAWjuNFGKhomSEReTJiOUfk8pTvjOTbfu7E1NppkvECAzwgbf3GwvLa6nUpybBUDHMGfDnTnXS9GETHIyveRmeWm28F1UTJSMqPQlOpjHZITQ5fOTleGJW7A3Nzv14bvMjIFzzMvK6Ds9tyA7gyyvh/ae4RWlNleUACbmZln4LgYmSkZUTIyMSHT3BsiJ4yWBkxIOqQvf9foDEUshBZMbwEkRz8v2Dvscz3T5KJzgeI4l+Vl9dfCMgZgoGXExaTJC6fb5yY5TYghCwnfD7IVEo9cf6NuXFEpeVv+B29kZXvKyM0YkfFeY43hKk/LMU4qFiZJhGEnR3RuImw4O/eG7VJUa8vl1QKID9M8LHE8pPztj2CssnO3y9WX6mSjFxkTJiItF8YxQun2BuOng0J8S3pGC8J2q4gtoxPBdcF4ZHsHrEcdTGmZv7kxXv6dUkpdlm2djkJAoichaEdkvIjUicneE17NF5An39W0iUhHy2j1u+34RWRPPpohUujYOuDazLmCML4tItYjsFpHHRSQnudszvulfyDZVMvrp6vUn5illpi58FzyIMFL4LhhWDL6Wm+WlfRiF0+cP0Nnjp6AvfJdNV28gZWHNdCfub5aIeIHvAzcCi4BPiMiisG53AC2qOg94ALjfvXYRsB5YDKwFfiAi3jg27wceUNUqoMW1PZgxpgN/ASxT1SWA1+1nJIglOhiRSNhTSmH2XfAgwvBD/gDy3UoTQXHKH+Y1pWD1iL5EhzwnDd2qOkQmEU9pOVCjqodUtQfYAKwL67MOeMx9/BRwvThfs9cBG1S1W1VrgRrXXkSb7jWrXRu4Nm8e5BgAGcAEEckAcoETCbxfIwzTJiOURD2lTK+HnEwPZ926byOJL3hku3fgPEvzHVEInpybl53RVwtvOAjWvQuG7ya6omTrSpFJRJSmA8dCnte5bRH7qKoPaANKYlwbrb0EaHVthI+V1Biqehz4DnAUOAm0qeoLkd6giNwpIjtEZEdDQ0PUG2EYhuP55IZksMWiMCeTM+dGPkwVPN48UqJDsKpCsM9wJzoEPaXQRAcwUYpGIqIUaYdc+HfnaH2Gqj3pMURkIo4XVQlMA/JE5LYIfVHVh1R1maouKysri9RlXGOOkhFKZ4+vL1kgHkUTMvs8hZHEFwiKUiRPyRGlYGWKghxHlIarmnmwGGvhhP5EB8D2KkUhEVGqA2aGPJ/BwDBYXx83VFYENMe4Nlp7I1Ds2ggfK9kxbgBqVbVBVXuBp4FrEni/hkt/RYcUT8RIKzqS8ZRSJEq9PueXNrwgK8DUYiffaVrxBMApNRTQfo9mqDkT7inlBz0ly8CLRCKitB2ocrPisnCSBTaG9dkI3O4+vgXYrM4q+UZgvZs5VwlUAW9Es+les8W1gWvzmUGOcRS4SkRy3bWn64G9id0WAyzRwYhMZ7ev7xyieBTmZPQdcDeS9Aaih+8umlLI7VfP5r4/XALApDzHgxmuIqnha0oF2RlkesU8pSjE/bqjqj4RuQvYhJPB9qiqVovIvcAOVd0IPAL8TERqcLyX9e611SLyJLAH8AFfUFU/QCSb7pBfBTaIyLeAna5tBjHGNhF5CnjLbd8JPDTYGzWesUP+jCCBgNLZ6+87KykeRRMyOdjQMcyzGkivP3r4zuMRvrFuSd/zibn9azzBk2GHkr7wnStKIsKkvKy+orDG+ST0m6WqzwHPhbV9LeRxF3BrlGvvA+5LxKbbfoj+7LnQ9sGM8XXg65GuMQwjebp8flRJ3FNK1ZqSm30XKSU8nGDiwXCJRDB8l5/T/3FrRVmjYxUdjLhYFM8IEqzOkKin5GTf9Y74kejB85FyEqjR1+cpDVP4rqWjh+LczPOOZi8ryKbhrK0pRcJEyYiLiZIRJFiFIFFPqWhCJgFl2Mv4hBM8STbemU8w/J5Sc0dP3xhBphTmcLKta1jGG+2YKBlRMS0ywunzlBLOvnP6nRmmzLZoJOMp5WZ5ycrwDFs4ramjuy8NPMi04gk0tHfTM8zHsI9GTJSMuFiigxGkz1PKTjT7zlncb+sc2XWlLncPUiKiJCKUF2Rz+szweC4tHb0DPKVpxTmoMmxjjmZMlAzDSJiOnqCnlJgopap6Qb+nlNhH3LTiCZwYpnBaU4TwXXCP1InWc8My5mjGRMmIi60pGUGCNeISDd+VuNUTGttHdlG/y5d4+A5gevGEYRGIQEBp6RwoSlOLXFFqM1EKx0TJMIyECdaIy0tQlMpSJErBEkI5CVQzByecdqqtC/8QZwm2nuvFH1BK8rIHjAdwotXCd+GYKBlRMQ/JCCdYnaHIreMWj8IJGWR5PTSO8DENQU8pO8Hw3dSiCfgCOuRp2kHvKyhCQXKzMijOzbTwXQRMlIy4mDgZQc50+RDpr+MWDxGhJD9r5MN3vYmnhANMn+iE0461dCY9lqry/S01/PtLBwaU5uoXpQkDrps9KZfDTSNf7SLdMVEyDCNhzpzrJT87A08ClRKClOZn0zTi4Ts/2RmekNOTYzOvLB+AA6fbkx7rraOtfHvTfr774nu8vP/8o29iidLcsnwOpaAEU7pjomTExVLCjSBt53oTDt0FcTylkc++SzTJAZxEhwmZXg7Un016rJf2ngacc5l+9vqR81473nqOrAzPgH1KAHPK8jjZ1mXHoodhomTExcJ3RpAz53r79h4lSml+9sgnOvgCCaeDg1OktWpy/qA8pX2nzrJwSgG3XTWb377XcN57ralvZ05pXkSPbY7rnZm3dD4mSkZUzEMywhmMpzS50KnzNtSZbbHo6PEnnCEYZF55PvtOnU36yJZDDe3MLcvnDy+bjj+gPLOr/7i5/a5gRWJOmVOR/FCjiVIoJkpGXEyajCBnupIXpenFufgCOqLVC9q7es+ryp0Il84sprG9m7qWxDPienwBjrWcY05ZHgumFHDJjCIef+MoqkpTezcn2rpYMKUw4rUVJXl4BA6cTj5kOJYxUTIMI2HazvX21bNLlBluZtvxEUx/bu/2Je0pXTF7IgBvHmlJ+JqjzR34A9rn9fzJVbOpqW/ntUNNvHKwCYCr55ZEvDYn00tVeQG7j7clNc+xjomSERc7gdYIMpjwXTDdum4Q6daDpb3bn7SntHBKIfnZGWw/3JzwNcEDDOeUOutDH71kGuUF2Xzz2b08urWW8oJslkyL7CkBLJlexLvHz9jfWAgmSkZc7M/FACejras3QHHuwEyyWEx306GPJxEWu1Dau3spSPDMpyBej7C8chK/P9CYsEgEkxSCnlJOppdv3ryE/afOsOtYK3etnkdGhNNvg1w8vZDG9m5On7GzlYIk979mjCvsy5sRSjCrrDQ/OVHKyfRSmp+V1FrNhdLe5UvaUwK44aLJbN5Xz/7TZ1kYZS0olNrGdsoKsikIyUhcs3gKG+9aSdu5Xq6JEroLcvGMYgDeqWtlStGUpOc7FjFPyYiLiZMB0OTuNQqv45YIMybmcqRpZMJ3qkp7t4/8JD0lgBsWlSMCm3afTqj/oYYO5pTmDWhfMr2Ia+eVxt28u3haIVleT9x1rG6fn49+byv3/mpPQvMazZgoGYaREH2eUkHyojSvPJ+ahuT3AA2Gbl+AXr+SNwhRKi/I4eo5JfzizWMJpbAfbGjv2280GHIyvVw6q5hX3aSIaLx3qp136tp49JVazrnHh4xVEhIlEVkrIvtFpEZE7o7weraIPOG+vk1EKkJeu8dt3y8ia+LZFJFK18YB12bWBYxRLCJPicg+EdkrIlcnd3sMB3OVjFBPKbnwHcD8yfk0nO2mtXP4KzuccYvGFg4ifAfwJytmU9dyjt+91xCzX0tHDy2dvRHxHCMAACAASURBVMwtG+gpJcPVc0qoPtEW8yDE2pAaee/UtV7QeOlOXFESES/wfeBGYBHwCRFZFNbtDqBFVecBDwD3u9cuAtYDi4G1wA9ExBvH5v3AA6paBbS4tpMew73m34DnVXUhcAmwN9EbY/Rj4TsDoKFvTSl5T6lqsrOB9L1BVExIlib3QMGSQcwT4IOLJjO5MJsfvFwTM+HhUKPzXuZcqCjNLSGgsK02urdUG1L1Ye/JMxc0XrqTiKe0HKhR1UOq2gNsANaF9VkHPOY+fgq4Xpxg6jpgg6p2q2otUOPai2jTvWa1awPX5s2DGUNECoH3A48AqGqPqo7trxhDjGmREUpjezd5WV4mJHjqbCjz+0Rp+DeKBk+5DT9YL1GyMjx8cXUV2w+3DCiwGsrBekco5l5A+A7gslnF5GV52RJjrMNNHUwryqFoQiYH6kcmDJoqEhGl6cCxkOd1blvEPqrqA9qAkhjXRmsvAVpdG+FjJTvGHKAB+LGI7BSRh0Uk4lcaEblTRHaIyI6Ghtgu+3jExMkAJ3w3mPUkgGlFORTkZLBnBL7l93lKgxQlgI9fOZOKkly+8avqqGs4e0+dYUKmlxkTcwc9DkB2hpdVC8p5cc/pqOtYhxo7qCzLY155vokSECl9JPzOReszVO2DGSMDuBz4oapeBnQAA9bDAFT1IVVdpqrLysrKInUxjHHPqbYuJhfmxO8YARHh0pnFvJVEtYTB0uyGGQfrKQFkej38w8cu5nBTJ9/etD9in13HWrl4ehHeJI7xiMaHFk+msb2bnUcj35/DjR1UlORRVZ5PjYkSdcDMkOczgBPR+ohIBlAENMe4Nlp7I1Ds2ggfazBj1KnqNrf9KRyRMpLE1pQMcMoEzYhwLlCiXD5rIvtPn+VsV/QF/aGgqaMHj5D0Jt9wrplbyqevns2jr9Tyq7fP/8g71+On+sQZLplZdEFjBPnAwnIyvcLzu08NeK2lo4e2c71UljqeUnNHz4ifTzWSJCJK24EqNysuCyepYGNYn43A7e7jW4DN6qwQbgTWu5lzlUAV8EY0m+41W1wbuDafGcwYqnoKOCYiC9xrrgfGfpL/MGAlUAyfP8CpM10RD6tLlCtmT0QV3j42vLXeTp/pojQ/e0g8mP/3pkVcWTGRv3nqbV6taexr/92BBnp8Aa6bX37BYwAU5mSyakE5/7PrBL3+wHmvBauIV5bmhazNjV1vKa4oues3dwGbcLLXnlTVahG5V0Q+6nZ7BCgRkRrgK7hhMlWtBp7EEYPngS+oqj+aTdfWV4GvuLZKXNtJj+Fe80XgP0XkHeBS4B+SvUHjGdMiI8hp9+iJYB27wXDZrGIyPMLva4Z33bau5VxfEdgLJSvDww9vu4JZk3L5zE+28/RbdQQCyqNbaynNz2bFnElDMg7AHy+bSWN7N78NS3g47IpSRYgoDeYwwtFCQon8qvoc8FxY29dCHncBt0a59j7gvkRsuu2HcLLzwtsHM8YuYFmka4zEMW0yYh3rnSgFOZmsmDOJzXvruefGi4ZqagOoaznHpTOLh8xeaX42G+68ms/9bAdfefJtvv5MNWe7ffzDH15MZoy6dsmyakEZpflZ/OLNY9ywaHJfe21jB16PMHNiLpleoSAnY0SyGFOFVXQw4mIekxEspjr9AkQJYPXCyRyob+foMJUc8geUE63nLsiji8SkvCw23Hk1/3zLUtYsmcJ3br2ETyyfGf/CJMj0evijy2fw0t76vi8B4FSNmD0pl6wMDyLCgskFvHdqHIfvDMMwahs78AgXHBZbs3gyIvDLnceHaGbnc7S5E19AqSi5sDTtSHg9wh8vm8l3br2EW66YEbeu3WD49DUVKPDo1tq+tpr680sZVU0u4L365E/IHS2YKBlxsWPRjZqGdmZOyiUnM/mNs6HMmJjLynmlPLkjsdpyyRKsdnDR1PgVvtOR6cUT+Ogl03j8jaO0dfbS1evncFMHVZP7RWn+5HxaO3v7KmyMNUyUjBiYGBkOB+vbL7hyQZBPLJ/F8dZzEdOfL5S9J8/gkf4KEqORz103h85eP//60nvsPNpKr1+5YtbEvtcXBDPwxmgIz0TJiI9p07jGH1AONXYwr3xoRGnN4inMLcvjgd+8N+Te0vbDzSycUnjBHl0qWTilkE8un8Vjrx7mG7+qJjvDw/KQLL/5UxxR2ndqbNbAM1Ey4mKaNL6pa+mkxxe44GrYQbwe4a8/tICa+vbz1k4ulHM9ft460sr7qkqHzGaq+OqNC1k0rZB9p87y2ZWVFIYcIlian83UohzeqRve/V6pwk6eNQwjJu8edz78Fk0dmuoFAGuXTOGDiybz7Rf2s2LOJJbOuPAU7hf2nKLHH+C6BaO/VFhhTiZPf/5ajrV0RjxE8NKZxbw9Ro+wME/JiMsYTfIxEuSdujayvB4WTBm6dRoR4R8/djHlBdl85sfb2X/qwvbdqCqPvXqYGRMncFVl7CPIRwtZGR7mluVHzPK7ZGYxR5o6+yqijyVMlIyomBgZAG8fa+WiaYVkZQztx0VpfjaPfXY5Ho/wRz98lY1vnxh0mvP/7DrOW0db+fyquXiGoLxQuhPcHPz2sbHnLZkoGXGxlPDxi88fYPfxNpZOH7rQXShzy/J55gvXMrcsj794fCfrH3qdzftO4wur/xYNVeXZd05wz9PvcmXFRD6+bGg3tKYrF08vwiOwcwyKkq0pGXExj2n88nZdGx09fq6aM3whsWnFE3j6z6/l568f4Qcv1/DZn+ygaEImK+eVsmhaIVXl+ZTkZ1Ocm4kq9LrFYfefOsvzu0+x61grl84s5oe3XUHGEJb9SWfysjNYOKWQ7bXNqZ7KkGOiZBhGVLYeaEQErpk7vOs0Xo9w+zUVfGL5LDbvO81v9tbz2sEmfv3uyZjXLZxSwDc+uphPLJ815OHFdOeauSX89PUjdPX6R3UKfDgmSkZczFEav2ytaWDJtCImXsCBecmQleFh7ZKprF0yFYAzXb0cbuygqaOHts5eRCDD42FyYTYVpXmU5g/uJNyxwNVzS3h4ay1vHW3hmrmjPw0+iImSERUTo/HN6TNd7DjSwl+srkrZHApzMockXXwscmXlJDwCrx1sGlOiNL78XWNQjNXCj0Zsnnv3JKrwkUumpnoqRgQKczK5eEYxrx1sSvVUhhQTJSMuJknjD1XllzuPs3BKAfPKR28dubHONXNL2HWsddiPmB9JTJQMwxjAzmOtvFPXxidXzEr1VIwYXDe/DF9AeSXkqPbRjomSER9zlcYdj2ytpSAngz+6fEaqp2LE4IrZEynIyWDLvuE9Yn4kMVEyomJrSeOT3cfb+PU7J/nUVbPJy7ZcqHQm0+vh/VVlbNlfP2b+Xk2UjLhYRYfxg6ryj/+7l+LcTD533dxUT8dIgFULyqg/2031ibFxlEVCoiQia0Vkv4jUiMjdEV7PFpEn3Ne3iUhFyGv3uO37RWRNPJsiUunaOODazBrsGO5rXhHZKSLPJn5bDOiP2o2RL2BGAjy54xiv1DTxlQ/Op2hCZvwLjJSzakE5AC/vr0/xTIaGuKIkIl7g+8CNwCLgEyKyKKzbHUCLqs4DHgDud69dBKwHFgNrgR+4IhHL5v3AA6paBbS4tpMeI2RufwnsTex2GKEEhuG4aiN9OdjQzjef3cvVc0q4bcXsVE/HSJCygmyWzihi875xIkrAcqBGVQ+pag+wAVgX1mcd8Jj7+CngenHqra8DNqhqt6rWAjWuvYg23WtWuzZwbd48yDEQkRnATcDDid0OI5Sgh2Se0tintbOHP31sB9kZHr5969JxUWl7LLFqQTk7j7WOiaMsEhGl6cCxkOd1blvEPqrqA9qAkhjXRmsvAVpdG+FjJTsGwL8CfwvELDksIneKyA4R2dHQMHayWC4U06LxQVtnL7c/+gbHW87x4KeuYMbE3FRPyUiS1QvLUYXfvjf6vaVERCnSV6bwz6tofYaqPekxROTDQL2qvhnh9fM7qz6kqstUdVlZ2eg/tXKoCLgukonT2KX+TBd/8sjr7D15lu//yeUsq5iU6ikZg2Dp9CJK87PYPAZSwxMRpTog9JCSGcCJaH1EJAMoAppjXButvREodm2Ej5XsGNcCHxWRwzjhwdUi8vME3q/h0h++M1kai7x1tIWPfG8rB+s7ePDTV/DBRZNTPSVjkHg8wqoF5fx2f33CZ1GlK4mI0nagys2Ky8JJKtgY1mcjcLv7+BZgszqfZBuB9W7mXCVQBbwRzaZ7zRbXBq7NZwYzhqreo6ozVLXCtb9ZVW9L8L4Y9HtKxtjC5w/w/S01fPzB18jK8PD0n1/DB9wMLmP0cv3Ccs50+XjzSEuqp3JBxN0Zp6o+EbkL2AR4gUdVtVpE7gV2qOpG4BHgZyJSg+O9rHevrRaRJ4E9gA/4gqr6ASLZdIf8KrBBRL4F7HRtM5gxjAujz1NK7TSMIWTPiTPc/fQ7vFPXxk0XT+VbNy8ZsWMpjOFlZVUpmV5h8/56VgzjoYzDjVho5nyWLVumO3bsSPU00oKqv3uOXr/y4KeuYM3iKamejnEB1J/t4l82vceTbx5jYm4W31y3hJuWWvXvscYn/+N1Gs528+JXrhvxsUXkTVVddqF2rIaIEZWApYSPelo6evjxK7U8vLWWXn+AO66t5IurqyjKtY2xY5HVC8v51q/3cqy5k5mTRmcWpYmSEZV+L9pUabTR2N7Nf/z+ED9/7QgdPX5uXDKFv127kMrSvFRPzRhGgqK0ZX89n766ItXTGRQmSkZUrKDD6KOm/iyPvXqEX7x5jB5fgA8vncYXPjCPBVPsTKTxwJyyfCpKcnlpr4mSMYax8F16EwgoW/bX85NXD/P7A41kZXhYd8k0Pr9qLnPK8lM9PWOEWb1wMj/fdoTOHh+5WaPvI370zdgwDACaO3p4+q06fvb6EY40dTKlMIe/WbOA9VfOpCQ/O9XTM1LE6oXlPPpKLa/WNHHDKNx7ZqJkxMUcpfQhEFC21jTyxI5jvFh9mh5/gGWzJ/I3axawZvEUMr12Gs14Z3nlJPKyvLy0r95EyRibWPgu9dS1dPKLHXU89WYdx1vPUZybyW1XzebjV8609SLjPLIyPLyvqoyX3YP/nLrVowcTJcNIU871+Hlx72l+seMYW2saAVg5r5R7/mAhH1w0mewMbxwLxnhl9cJynq8+xd6TZ1k0rTDV00kKEyUjLnby7MjhDyivHmzkf3ae4PndJ+no8TOtKIe/WF3FrctmWAVvIyFWLXQKS2/ed9pEyTCM5FBVqk+c4X92Hmfj2yeoP9tNQXYGH146jXWXTeOqyhI738hIivKCnL6D/+5aXZXq6SSFiZIRF1tTGh6ONXey8e0T/HLncWrq28n0Ch9YUM7Nl01n9cJycjItPGcMng8sKOffNx+guaOHSaOovqGJkhEX06Sho/5MF8+9e5Jn3znJDrea85UVE7nvD5dw08VTKc4dPR8eRnqzemE5//bSAX73XgM3XxZ+Lmv6YqJkGMNMw9lunt/tCNEbh5tRhQWTC/jrD81n3aXTR22NMiO9WTK9iOLcTLbWNJooGWMLqySfPE3t3TxffYpn3z7JttomAgrzyvP5y+uruOniqVRNtjRuY3jxeoRr5pbwSk3jqEoNN1EyjCGipaOH56tP8et3TvLaoSb8AWVOWR53ra7iw0unMt+EyBhhVs4r47l3T3GosYO5o6TklImSYVwADWe7eXHPaZ6vPsUrNY34A0pFSS6fv24uNy2dysIpBaPmG6ox9lg5rxSAV2oaTZSM0U0gpES4Re/O51hzJ5uqT7Gp+hQ7jrSgCrNLcvnc++dw09KpLJpaaEJkpAWzSnKZOWkCWw80jpqq4SZKRkR6A4FUTyFtUFXeO93OpupTPL/7FHtOngHgoqmFfOn6+axZMpkFk80jMtKTlfNKefadk/j8ATJGQW1EEyUjIj5/iKc0DpPCAwFlV12r4xHtPsXhpk5E4IpZE/m7P7iINYunMKvEsuaM9OfaeaU8/sYx3j3exmWzJqZ6OnExURrjBDPnkv0W3+sff55Sjy/AG7XNbKo+xQt7TnH6TDcZHuGaeaX82fvn8MFFkykvyEn1NA0jKa6Z66wrbT3QOHZESUTWAv8GeIGHVfWfwl7PBn4KXAE0AR9X1cPua/cAdwB+4C9UdVMsmyJSCWwAJgFvAZ9S1Z5kxxCRmW7/KUAAeEhV/y3ZGzTaqbznOW65YgbfufWSpK7r9Y+PNaX6s128vK+Bzfvq2VrTSHu3jwmZXq6bX8baJVP4wMJyiiZkpnqahjFoJuVlsXhaIVtrGvni9elfciiuKImIF/g+8EGgDtguIhtVdU9ItzuAFlWdJyLrgfuBj4vIImA9sBiYBvxGROa710SzeT/wgKpuEJEfubZ/OIgxfMBfqepbIlIAvCkiL4bNe0wT9JKeerMuaVHq6vWH2BnSaaWUQEB593gbm/fVs2V/Pe/UtQEwpTCHj1wyjdULy1k5r5QJWVbixxg7rJxXyqOv1HKux5/2v9uJeErLgRpVPQQgIhuAdUDoh/s64P9zHz8FfE+ceNE6YIOqdgO1IlLj2iOSTRHZC6wGPun2ecy1+8Nkx1DV14CTAKp61rU9PWzeY5pu3+BDcJ09/vidRglnu3rZeqDRFaIGGtu7EYHLZhbzN2sW8IEF5Vw01RIVjLHLVXNKePB3h9h5tIVr3DTxdCURUZoOHAt5XgesiNZHVX0i0gaUuO2vh10brHcRyWYJ0Kqqvgj9BzMGACJSAVwGbIv0BkXkTuBOgFmzZkXqMio52+WL3ykKHT39145GR+lQQ3ufN/RGbTO9fqUwJ4PrFpSzemEZ180vH1VFKg3jQlhWMRGPwOu1zWNClCJ9fQz/nIrWJ1p7pLzEWP0HM4ZzkUg+8N/Al1T1TIS+qOpDwEMAy5YtG42fwRHp6B68KHV2jy5PKZikEBSi2sYOAKrK8/nsykpWLyjnitkTR0VKrGEMNQU5mSyeVsS2Q02pnkpcEhGlOmBmyPMZwIkofepEJAMoAprjXBupvREoFpEM11sK7Z/0GCKSiSNI/6mqTyfwXscU7RcgSud5Smm6qBSapPD7Aw109PjJyvBwzdwSPnNtBR9YUG7FTg3DZUXlJH76+hG6ev1pfSxKIqK0Hahys+KO4yQVfDKsz0bgduA14BZgs6qqiGwE/ktEvouThFAFvIHj3Qyw6V6zxbWxwbX5zGDGcNebHgH2qup3k70xY4ELWRcK9bLSRZJCkxQ276vn3eNOksLUohzWXTad1QvKuWZeCblZttPBMMJZMaeEh7fW8vaxVlbMKUn1dKIS96/XXb+5C9iEk779qKpWi8i9wA5V3Yjz4f8zN8mgGUdkcPs9iZNc4AO+oKp+gEg23SG/CmwQkW8BO13bJDuGiKwEPgW8KyK7XBv/j6o+N7hbNfoICkumN/kF/Mb27qGezqAIJim8tK+el90kBY/AZbMmWpKCYSTB8opJiMC22ubRLUoA7gf5c2FtXwt53AXcGuXa+4D7ErHpth+iP0MvtD2pMVR1K5HXm8YNwRBcdkbyrnpje0//kxF2lYJJCpv31bP9sCUpGMZQUJSbyYLJBWyrbcIJKKUnFucYwwSz77Izkl/cbzjrpE2rDn+ZoW6fvz9JYV89h5s6AZg/2UlSuH7hZC6fVWxJCoZxgVw1p4QN24/S4wuQNYjPhZHARGkM03jWCcEVDqIiQW1jBxUleX1ZbENN/Zkutux3vKGtBxrPS1L47MpKS1IwjGFgReUkfvLqYd493soVsyelejoRMVEawzS460LJLreoKjX17Vw9t4Taxo4hqegQCCi7T7Txm72ONxSepHD9wnKumWuVFAxjOFle6QjR64eaTZSMkedEaxeQfHHVPSfP0N7tY+n0Il7cc3rQ45/r8bO1ppGX9p5m87566s+en6SwemG5HYJnGCNISX42VeX5bKtt5gsfSPVsImOiNIZ57/RZ4PxjKBJh0+5TALx/fhn/8uJ7Sa0oNZzt5oU9p3hpbz2v1DTS7QuQn53BdfPLuP6iclYtsCQFw0glK+ZM4pdvHU/b85VMlMYoJ1rPcbTZSRjoTUKUWjp6+Pm2o1y/sJzJhc4xDfHCd43t3Ty/+xS/fuck22qbCCjMmpTLJ1fM4vqFk1leOSltF1UNY7yxorKEn79+lN0nznDpzOJUT2cAJkpjEFXl/2yuQQTeV1XG28daE7qu2+fnS0/s4sy5Xr7yofl4PU5YzR/hFFqfP8DL+xt4/I2jbNlfT0BhTlked62u4qaLpzJ/cr6F5QwjDVkxx1lL2naoyUTJGF56/QF+9fYJHvrdIfadOsufrqxEgTcPN8e99lhzJ3+5YSdvHW3lHz92MYunFXGmqxc4v9p4jy/AL948xg+2HOR46znKCrL53HVzWXfpNDsS3DBGAeUFOcwpy+P1Q0187rq5qZ7OAEyUxgCdPT6e2H6Mh39fy/HWc8yfnM93//gSbr50Ov+8aX/M8F392S5+9PIhfr7tCNkZHr73ycv48NJpAGS58eagKO061spfPbmLgw0dXD6rmL//8CKuv6iczDSMSxuGEZ0VlSU8+/YJ/AHti4ikCyZKo5jWzh4ee/UIP3m1lpbOXq6smMi96xbzgQXleNxftEyv0BsIoKrneTHNHT08+NuDPPbaYXr9yh9dPp0v3TCfacUT+voERanHF+CF6lPc9V87Kc3P4tH/axkfWFBuXpFhjFKumjOJx984yp4TZ7h4RlGqp3MeJkqjkLqWTn7yymH+642jdPb4Wb2wnM+vmsuVFQP3HWR4PKiCP6BkeIW2c708/PtDPLq1lnO9fm6+dDpfvL6KytK8Add6PEKW10P1iTM8+LuDXDStkMc+cyXFuZY9Zxijmavc2nevH2oyUTIGR68/wG/2nObx7cf4/YEGPCJ8ZOlU/u9Vc1k4pTDqdZkZjjfjCyjP7Krj3mf30Haul5uWTuXLN1Qxr7wg5rhZGR5+s/c0xbmZ/MenrzBBMowxwOTCHCpL89hW28SfvX9OqqdzHiZKaUyPL8Brh5p4fvcpXqg+RVNHD1OLcvji6ir+eNkMZkyMX4Yn0+OE4O5/fh8/fuUwV1ZM5BsfXcKiadGFLJTgmUx/9r45lBfkDP7NGIaRVqyonMSv3z2ZdutKJkppRmePj9/ub2BT9Sle2lfP2S4feVleVi0s548un85188uT+gUKHlvx41cO88fLZvCPH1ua1PULpxSw79RZbrtqdtLvxTCM9GXFnEls2H6MvSfPsGR6+oTwTJTSgLbOXn6z9zSbqk/x2/ca6PYFmJibyY1LprBm8RSunVc66JMigzu2s7we7l23JOlvRA9+6go6e/wUDaKoq2EY6cuKSmddaVtts4mS4aRiv1DtCNFrB5vwBZSpRTl8Yvks1iyewpUVE4ekBMiCKc6a0V+vmT8oYZtdMjABwjCM0c+04gnMmpTL64eauGNlZaqn04eJ0ghyvPUcz+8+xfO7T7LjSAuqMKc0jz97/xzWLp7C0hlFQ55mfWXFJF788vuZV54/pHYNwxj9XDVnEpuqTxMIaN82klRjojTMnOvx87+7T/LkjmO8fsiprLBwSgFfun4+N148hary4S/HUzU5doadYRjjkxWVJTy5o459p84mnPw03JgoDRMtHT38+JVafvLqYc50+ZhdkstffXA+H75kWsQ9QYZhGCPNVXOddaVXDzaaKI1V/AHl568f4dub9tPe7WPN4sl85tpKlldMShv32DAMA2B68QQWTingherT/On70mO/UkIr6SKyVkT2i0iNiNwd4fVsEXnCfX2biFSEvHaP275fRNbEsykila6NA67NrKEeY7ho7ezh9kff4Osbq7lsVjGbvvR+HvzUMq6aU2KCZBhGWvKhxVPYfqSZhrPdqZ4KkIAoiYgX+D5wI7AI+ISILArrdgfQoqrzgAeA+91rFwHrgcXAWuAHIuKNY/N+4AFVrQJaXNtDPcaQ09LRw60/eo03apv5p49dzE8/u7wv880wDCNdWbt4Cqrwq7dPpHoqQGKe0nKgRlUPqWoPsAFYF9ZnHfCY+/gp4HpxVu/XARtUtVtVa4Ea115Em+41q10buDZvHsoxErstydHjC/CZn2znSHMnP/nslaxfPsuKlRqGMSq4aGoBV1ZM5P9sPsCJ1nOpnk5Ca0rTgWMhz+uAFdH6qKpPRNqAErf99bBrp7uPI9ksAVpV1Reh/1CNMQARuRO4033aLiJNQGOkvvG49r7BXJW2lDLI+zAGsXvRj90LhzF3H6Z/fdCXlgJDUvYlEVGK9JU//ICeaH2itUfy0GL1H8oxBjaqPgQ8FHwuIjtUdVmkvuMJuw/92L3ox+6Fg92Hftx7UTEUthIJ39UBM0OezwDCg499fUQkAygCmmNcG629ESh2bYSPNVRjGIZhGGlKIqK0Hahys+KycJIKNob12Qjc7j6+Bdisquq2r3cz5yqBKuCNaDbda7a4NnBtPjOUYyR2WwzDMIxUEDd8567f3AVsArzAo6paLSL3AjtUdSPwCPAzEanB8V7Wu9dWi8iTwB7AB3xBVf0AkWy6Q34V2CAi3wJ2urYZ4jHi8VD8LuMCuw/92L3ox+6Fg92HfobsXojjbBiGYRhG6rnwMtSGYRiGMUSYKBmGYRhpg4lSCCNdligViMijIlIvIrtD2iaJyItuaacXRWSi2y4i8u/u/XhHRC4PueZ2t/8BEbk90ljpjIjMFJEtIrJXRKpF5C/d9vF4L3JE5A0Redu9F99w24es5Ndowq0Is1NEnnWfj9f7cFhE3hWRXSKyw20b/r8PVbV/zrqaFzgIzAGygLeBRame1zC8z/cDlwO7Q9r+GbjbfXw3cL/7+A+A/8XZC3YVsM1tnwQccn9OdB9PTPV7S/I+TAUudx8XAO/hlKMaj/dCgHz3cSawzX2PTwLr3fYfAZ93H/858CP38XrgfK66SwAAArdJREFUCffxIvfvJhuodP+evKl+f4O4H18B/gt41n0+Xu/DYaA0rG3Y/z7MU+pnxMoSpRJV/R1O9mIooSWcwks7/VQdXsfZQzYVWAO8qKrNqtoCvIhTd3DUoKonVfUt9/FZYC9OJZDxeC9UVdvdp5nuP2XoSn6NGkRkBnAT8LD7fChLn40Fhv3vw0Spn0jllKZH6TvWmKyqJ8H5sAbK3fZo92RM3Ss37HIZjocwLu+FG7LaBdTjfHAcJMGSX0Boya/Rfi/+FfhbIOA+T7j0GWPrPoDzxeQFEXlTnFJsMAJ/H3aeUj+JlFMabyRb2mnUISL5wH8DX1LVMxK9kO6Yvhfq7O27VESKgV8CF0Xq5v4ck/dCRD4M1KvqmyKyKtgcoeuYvg8hXKuqJ0SkHHhRRPbF6Dtk98I8pX7Gc1mi066rjfuz3m0f0yWcRCQTR5D+U1WfdpvH5b0IoqqtwMs46wJDVfJrtHAt8FEROYwTvl+N4zmNt/sAgKqecH/W43xRWc4I/H2YKPUznssShZZwCi/t9Gk3s+YqoM112TcBHxKRiW72zYfctlGDG/t/BNirqt8NeWk83osy10NCRCYAN+CssQ1Vya9Rgareo6oz1Cksuh7nff0J4+w+AIhInogUBB/j/F7vZiT+PlKd4ZFO/3AySN7Diaf/XarnM0zv8XHgJNCL8y3mDpw4+EvAAffnJLev4ByUeBB4F1gWYuezOAu4NcBnUv2+BnEfVuKEEd4Bdrn//mCc3oulOCW93nE/eL7mts/B+TCtAX4BZLvtOe7zGvf1OSG2/s69R/uBG1P93i7gnqyiP/tu3N0H9z2/7f6rDn4ejsTfh5UZMgzDMNIGC98ZhmEYaYOJkmEYhpE2mCgZhmEYaYOJkmEYhpE2mCgZhmEYaYOJkmEYhpE2mCgZhmEYacP/Dzlc4oUjVCBMAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.clf()\n",
"# plt.plot(x_part, calcs, '.')\n",
"plt.plot(test_q, calcs_test, label = 'pdf')\n",
"# plt.plot(test_q, f0_y, label = '0')\n",
"# plt.plot(test_q, fT_y, label = 'T')\n",
"# plt.plot(test_q, fplus_y, label = '+')\n",
"# plt.plot(test_q, res_y, label = 'res')\n",
"plt.legend()\n",
"plt.ylim(0.0, 1.5e-6)\n",
"# plt.yscale('log')\n",
"# plt.xlim(770, 785)\n",
"plt.savefig('test.png')\n",
"# print(jpsi_width)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"\n",
"\n",
"# probs = mixture.prob(test_q)\n",
"# probs_np = zfit.run(probs)\n",
"# probs_np *= np.max(calcs_test) / np.max(probs_np)\n",
"# plt.figure()\n",
"# plt.semilogy(test_q, probs_np,label=\"importance sampling\")\n",
"# plt.semilogy(test_q, calcs_test, label = 'pdf')\n"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"# 0.213/(0.00133+0.213+0.015)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Adjust scaling of different parts"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"total_f.update_integration_options(draws_per_dim=2000000, mc_sampler=None)\n",
"# inte = total_f.integrate(limits = (950., 1050.), norm_range=False)\n",
"# inte_fl = zfit.run(inte)\n",
"# print(inte_fl/4500)\n",
"# print(pdg[\"jpsi_BR\"]/pdg[\"NR_BR\"], inte_fl*pdg[\"psi2s_auc\"]/pdg[\"NR_auc\"])"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"# # print(\"jpsi:\", inte_fl)\n",
"# # print(\"Increase am by factor:\", np.sqrt(pdg[\"jpsi_BR\"]/pdg[\"NR_BR\"]*pdg[\"NR_auc\"]/inte_fl))\n",
"# # print(\"New amp:\", pdg[\"jpsi\"][3]*np.sqrt(pdg[\"jpsi_BR\"]/pdg[\"NR_BR\"]*pdg[\"NR_auc\"]/inte_fl))\n",
"\n",
"# # print(\"psi2s:\", inte_fl)\n",
"# # print(\"Increase am by factor:\", np.sqrt(pdg[\"psi2s_BR\"]/pdg[\"NR_BR\"]*pdg[\"NR_auc\"]/inte_fl))\n",
"# # print(\"New amp:\", pdg[\"psi2s\"][3]*np.sqrt(pdg[\"psi2s_BR\"]/pdg[\"NR_BR\"]*pdg[\"NR_auc\"]/inte_fl))\n",
"\n",
"# name = \"phi\"\n",
"\n",
"# print(name+\":\", inte_fl)\n",
"# print(\"Increase am by factor:\", np.sqrt(pdg[name+\"_BR\"]/pdg[\"NR_BR\"]*pdg[\"NR_auc\"]/inte_fl))\n",
"# print(\"New amp:\", pdg[name][0]*np.sqrt(pdg[name+\"_BR\"]/pdg[\"NR_BR\"]*pdg[\"NR_auc\"]/inte_fl))\n",
"\n",
"\n",
"# print(x_min)\n",
"# print(x_max)\n",
"# # total_f.update_integration_options(draws_per_dim=2000000, mc_sampler=None)\n",
"# total_f.update_integration_options(mc_sampler=lambda dim, num_results,\n",
"# dtype: tf.random_uniform(maxval=1., shape=(num_results, dim), dtype=dtype),\n",
"# draws_per_dim=1000000)\n",
"# # _ = []\n",
"\n",
"# # for i in range(10):\n",
"\n",
"# # inte = total_f.integrate(limits = (x_min, x_max))\n",
"# # inte_fl = zfit.run(inte)\n",
"# # print(inte_fl)\n",
"# # _.append(inte_fl)\n",
"\n",
"# # print(\"mean:\", np.mean(_))\n",
"\n",
"# _ = time.time()\n",
"\n",
"# inte = total_f.integrate(limits = (x_min, x_max))\n",
"# inte_fl = zfit.run(inte)\n",
"# print(inte_fl)\n",
"# print(\"Time taken: {}\".format(display_time(int(time.time() - _))))\n",
"\n",
"# print(pdg['NR_BR']/pdg['NR_auc']*inte_fl)\n",
"# print(0.25**2*4.2/1000)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Sampling\n",
"## Mixture distribution for sampling"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"\n",
" \n",
"# print(list_of_borders[:9])\n",
"# print(list_of_borders[-9:])\n",
"\n",
"\n",
"class UniformSampleAndWeights(zfit.util.execution.SessionHolderMixin):\n",
" def __call__(self, limits, dtype, n_to_produce):\n",
" # n_to_produce = tf.cast(n_to_produce, dtype=tf.int32)\n",
" low, high = limits.limit1d\n",
" low = tf.cast(low, dtype=dtype)\n",
" high = tf.cast(high, dtype=dtype)\n",
"# uniform = tfd.Uniform(low=low, high=high)\n",
"# uniformjpsi = tfd.Uniform(low=tf.constant(3080, dtype=dtype), high=tf.constant(3112, dtype=dtype))\n",
"# uniformpsi2s = tfd.Uniform(low=tf.constant(3670, dtype=dtype), high=tf.constant(3702, dtype=dtype))\n",
"\n",
"# list_of_borders = []\n",
"# _p = []\n",
"# splits = 10\n",
"\n",
"# _ = np.linspace(x_min, x_max, splits)\n",
"\n",
"# for i in range(splits):\n",
"# list_of_borders.append(tf.constant(_[i], dtype=dtype))\n",
"# _p.append(tf.constant(1/splits, dtype=dtype))\n",
" \n",
"# mixture = tfd.MixtureSameFamily(mixture_distribution=tfd.Categorical(probs=_p[:(splits-1)]),\n",
"# components_distribution=tfd.Uniform(low=list_of_borders[:(splits-1)], \n",
"# high=list_of_borders[-(splits-1):]))\n",
" mixture = tfd.MixtureSameFamily(mixture_distribution=tfd.Categorical(probs=[tf.constant(0.05, dtype=dtype),\n",
" tf.constant(0.93, dtype=dtype),\n",
" tf.constant(0.05, dtype=dtype),\n",
" tf.constant(0.065, dtype=dtype),\n",
" tf.constant(0.04, dtype=dtype),\n",
" tf.constant(0.05, dtype=dtype)]),\n",
" components_distribution=tfd.Uniform(low=[tf.constant(x_min, dtype=dtype), \n",
" tf.constant(3090, dtype=dtype),\n",
" tf.constant(3681, dtype=dtype), \n",
" tf.constant(3070, dtype=dtype),\n",
" tf.constant(1000, dtype=dtype),\n",
" tf.constant(3660, dtype=dtype)], \n",
" high=[tf.constant(x_max, dtype=dtype),\n",
" tf.constant(3102, dtype=dtype), \n",
" tf.constant(3691, dtype=dtype),\n",
" tf.constant(3110, dtype=dtype),\n",
" tf.constant(1040, dtype=dtype),\n",
" tf.constant(3710, dtype=dtype)]))\n",
"# dtype = tf.float64\n",
"# mixture = tfd.MixtureSameFamily(mixture_distribution=tfd.Categorical(probs=[tf.constant(0.04, dtype=dtype),\n",
"# tf.constant(0.90, dtype=dtype),\n",
"# tf.constant(0.02, dtype=dtype),\n",
"# tf.constant(0.07, dtype=dtype),\n",
"# tf.constant(0.02, dtype=dtype)]),\n",
"# components_distribution=tfd.Uniform(low=[tf.constant(x_min, dtype=dtype), \n",
"# tf.constant(3089, dtype=dtype),\n",
"# tf.constant(3103, dtype=dtype), \n",
"# tf.constant(3681, dtype=dtype),\n",
"# tf.constant(3691, dtype=dtype)], \n",
"# high=[tf.constant(3089, dtype=dtype),\n",
"# tf.constant(3103, dtype=dtype), \n",
"# tf.constant(3681, dtype=dtype),\n",
"# tf.constant(3691, dtype=dtype), \n",
"# tf.constant(x_max, dtype=dtype)]))\n",
"# mixture = tfd.Uniform(tf.constant(x_min, dtype=dtype), tf.constant(x_max, dtype=dtype))\n",
"# sample = tf.random.uniform((n_to_produce, 1), dtype=dtype)\n",
" sample = mixture.sample((n_to_produce, 1))\n",
"# sample = tf.random.uniform((n_to_produce, 1), dtype=dtype)\n",
" weights = mixture.prob(sample)[:,0]\n",
"# weights = tf.broadcast_to(tf.constant(1., dtype=dtype), shape=(n_to_produce,))\n",
" # sample = tf.expand_dims(sample, axis=-1)\n",
"# print(sample, weights)\n",
" \n",
"# weights = tf.ones(shape=(n_to_produce,), dtype=dtype)\n",
" weights_max = None\n",
" thresholds = tf.random_uniform(shape=(n_to_produce,), dtype=dtype)\n",
" return sample, thresholds, weights, weights_max, n_to_produce"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"# total_f._sample_and_weights = UniformSampleAndWeights"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [],
"source": [
"# 0.00133/(0.00133+0.213+0.015)*(x_max-3750)/(x_max-x_min)"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"# zfit.settings.set_verbosity(10)"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"# # zfit.run.numeric_checks = False \n",
"\n",
"# nr_of_toys = 1\n",
"# nevents = int(pdg[\"number_of_decays\"])\n",
"# nevents = pdg[\"number_of_decays\"]\n",
"# event_stack = 1000000\n",
"# # zfit.settings.set_verbosity(10)\n",
"# calls = int(nevents/event_stack + 1)\n",
"\n",
"# total_samp = []\n",
"\n",
"# start = time.time()\n",
"\n",
"# sampler = total_f.create_sampler(n=event_stack)\n",
"\n",
"# for toy in range(nr_of_toys):\n",
" \n",
"# dirName = 'data/zfit_toys/toy_{0}'.format(toy)\n",
" \n",
"# if not os.path.exists(dirName):\n",
"# os.mkdir(dirName)\n",
"# print(\"Directory \" , dirName , \" Created \")\n",
"\n",
"# for call in range(calls):\n",
"\n",
"# sampler.resample(n=event_stack)\n",
"# s = sampler.unstack_x()\n",
"# sam = zfit.run(s)\n",
"# # clear_output(wait=True)\n",
"\n",
"# c = call + 1\n",
" \n",
"# print(\"{0}/{1} of Toy {2}/{3}\".format(c, calls, toy+1, nr_of_toys))\n",
"# print(\"Time taken: {}\".format(display_time(int(time.time() - start))))\n",
"# print(\"Projected time left: {}\".format(display_time(int((time.time() - start)/(c+calls*(toy))*((nr_of_toys-toy)*calls-c)))))\n",
"\n",
"# with open(\"data/zfit_toys/toy_{0}/{1}.pkl\".format(toy, call), \"wb\") as f:\n",
"# pkl.dump(sam, f, pkl.HIGHEST_PROTOCOL)"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"# with open(r\"data/zfit_toys/toy_0/0.pkl\", \"rb\") as input_file:\n",
"# sam = pkl.load(input_file)\n",
"# print(sam[:10])\n",
"\n",
"# with open(r\"data/zfit_toys/toy_0/1.pkl\", \"rb\") as input_file:\n",
"# sam2 = pkl.load(input_file)\n",
"# print(sam2[:10])\n",
"\n",
"# print(np.sum(sam-sam2))"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [],
"source": [
"# print(\"Time to generate full toy: {} s\".format(int(time.time()-start)))\n",
"\n",
"# total_samp = []\n",
"\n",
"# for call in range(calls):\n",
"# with open(r\"data/zfit_toys/toy_0/{}.pkl\".format(call), \"rb\") as input_file:\n",
"# sam = pkl.load(input_file)\n",
"# total_samp = np.append(total_samp, sam)\n",
"\n",
"# total_samp = total_samp.astype('float64')\n",
"\n",
"# data2 = zfit.data.Data.from_numpy(array=total_samp[:int(nevents)], obs=obs)\n",
"\n",
"# data3 = zfit.data.Data.from_numpy(array=total_samp, obs=obs)\n",
"\n",
"# print(total_samp[:nevents].shape)"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [],
"source": [
"# plt.clf()\n",
"\n",
"# bins = int((x_max-x_min)/7)\n",
"\n",
"# # calcs = zfit.run(total_test_tf(samp))\n",
"# print(total_samp[:nevents].shape)\n",
"\n",
"# plt.hist(total_samp[:nevents], bins = bins, range = (x_min,x_max), label = 'data')\n",
"# # plt.plot(test_q, calcs_test*nevents , label = 'pdf')\n",
"\n",
"# # plt.plot(sam, calcs, '.')\n",
"# # plt.plot(test_q, calcs_test)\n",
"# # plt.yscale('log')\n",
"# plt.ylim(0, 200)\n",
"# # plt.xlim(3080, 3110)\n",
"\n",
"# plt.legend()\n",
"\n",
"# plt.savefig('test2.png')"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [],
"source": [
"# sampler = total_f.create_sampler(n=nevents)\n",
"# nll = zfit.loss.UnbinnedNLL(model=total_f, data=sampler, fit_range = (x_min, x_max))\n",
"\n",
"# # for param in pdf.get_dependents():\n",
"# # param.set_value(initial_value)\n",
"\n",
"# sampler.resample(n=nevents)\n",
"\n",
"# # Randomise initial values\n",
"# # for param in pdf.get_dependents():\n",
"# # param.set_value(random value here)\n",
"\n",
"# # Minimise the NLL\n",
"# minimizer = zfit.minimize.MinuitMinimizer(verbosity = 10)\n",
"# minimum = minimizer.minimize(nll)"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [],
"source": [
"# jpsi_width"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [],
"source": [
"# plt.hist(sample, weights=1 / prob(sample))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Fitting"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [],
"source": [
"# start = time.time()\n",
"\n",
"# for param in total_f.get_dependents():\n",
"# param.randomize()\n",
" \n",
"# # for param in total_f.get_dependents():\n",
"# # print(zfit.run(param))\n",
" \n",
"# nll = zfit.loss.UnbinnedNLL(model=total_f, data=data2, fit_range = (x_min, x_max))\n",
"\n",
"# minimizer = zfit.minimize.MinuitMinimizer(verbosity = 5)\n",
"# # minimizer._use_tfgrad = False\n",
"# result = minimizer.minimize(nll)\n",
"\n",
"# # param_errors = result.error()\n",
"\n",
"# # for var, errors in param_errors.items():\n",
"# # print('{}: ^{{+{}}}_{{{}}}'.format(var.name, errors['upper'], errors['lower']))\n",
"\n",
"# print(\"Function minimum:\", result.fmin)\n",
"# # print(\"Results:\", result.params)\n",
"# print(\"Hesse errors:\", result.hesse())"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [],
"source": [
"# print(\"Time taken for fitting: {}\".format(display_time(int(time.time()-start))))\n",
"\n",
"# # probs = total_f.pdf(test_q)\n",
"\n",
"# calcs_test = zfit.run(probs)\n",
"# res_y = zfit.run(jpsi_res(test_q))"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [],
"source": [
"# plt.clf()\n",
"# # plt.plot(x_part, calcs, '.')\n",
"# plt.plot(test_q, calcs_test, label = 'pdf')\n",
"# # plt.plot(test_q, res_y, label = 'res')\n",
"# plt.legend()\n",
"# plt.ylim(0.0, 10e-6)\n",
"# # plt.yscale('log')\n",
"# # plt.xlim(3080, 3110)\n",
"# plt.savefig('test3.png')\n",
"# # print(jpsi_width)"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [],
"source": [
"# _tot = 4.37e-7+6.02e-5+4.97e-6\n",
"# _probs = []\n",
"# _probs.append(6.02e-5/_tot)\n",
"# _probs.append(4.97e-6/_tot)\n",
"# _probs.append(4.37e-7/_tot)\n",
"# print(_probs)"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [],
"source": [
"# dtype = 'float64'\n",
"# # mixture = tfd.Uniform(tf.constant(x_min, dtype=dtype), tf.constant(x_max, dtype=dtype))\n",
"# mixture = tfd.MixtureSameFamily(mixture_distribution=tfd.Categorical(probs=[tf.constant(0.007, dtype=dtype),\n",
"# tf.constant(0.917, dtype=dtype),\n",
"# tf.constant(0.076, dtype=dtype)]),\n",
"# components_distribution=tfd.Uniform(low=[tf.constant(x_min, dtype=dtype), \n",
"# tf.constant(3080, dtype=dtype),\n",
"# tf.constant(3670, dtype=dtype)], \n",
"# high=[tf.constant(x_max, dtype=dtype),\n",
"# tf.constant(3112, dtype=dtype), \n",
"# tf.constant(3702, dtype=dtype)]))\n",
"# # for i in range(10):\n",
"# # print(zfit.run(mixture.prob(mixture.sample((10, 1)))))"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [],
"source": [
"# print((zfit.run(jpsi_p)%(2*np.pi))/np.pi)\n",
"# print((zfit.run(psi2s_p)%(2*np.pi))/np.pi)"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [],
"source": [
"# def jpsi_res(q):\n",
"# return resonance(q, _mass = jpsi_mass, scale = jpsi_scale,\n",
"# phase = jpsi_phase, width = jpsi_width)\n",
"\n",
"# def psi2s_res(q):\n",
"# return resonance(q, _mass = psi2s_mass, scale = psi2s_scale,\n",
"# phase = psi2s_phase, width = psi2s_width)\n",
" \n",
"# def p3770_res(q):\n",
"# return resonance(q, _mass = p3770_mass, scale = p3770_scale,\n",
"# phase = p3770_phase, width = p3770_width)\n",
" \n",
"# def p4040_res(q):\n",
"# return resonance(q, _mass = p4040_mass, scale = p4040_scale,\n",
"# phase = p4040_phase, width = p4040_width)\n",
" \n",
"# def p4160_res(q):\n",
"# return resonance(q, _mass = p4160_mass, scale = p4160_scale,\n",
"# phase = p4160_phase, width = p4160_width)\n",
" \n",
"# def p4415_res(q):\n",
"# return resonance(q, _mass = p4415_mass, scale = p4415_scale,\n",
"# phase = p4415_phase, width = p4415_width)"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [],
"source": [
"# 0.15**2*4.2/1000\n",
"# result.hesse()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Constraints"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [],
"source": [
"# 1. Constraint - Real part of sum of Psi contrib and D contribs\n",
"\n",
"sum_list = []\n",
"\n",
"sum_list.append(ztf.to_complex(jpsi_s) * tf.exp(tf.complex(ztf.constant(0.0), jpsi_p)) * ztf.to_complex(jpsi_w / (tf.pow(jpsi_m,3))))\n",
"sum_list.append(ztf.to_complex(psi2s_s) * tf.exp(tf.complex(ztf.constant(0.0), psi2s_p)) * ztf.to_complex(psi2s_w / (tf.pow(psi2s_m,3))))\n",
"sum_list.append(ztf.to_complex(p3770_s) * tf.exp(tf.complex(ztf.constant(0.0), p3770_p)) * ztf.to_complex(p3770_w / (tf.pow(p3770_m,3))))\n",
"sum_list.append(ztf.to_complex(p4040_s) * tf.exp(tf.complex(ztf.constant(0.0), p4040_p)) * ztf.to_complex(p4040_w / (tf.pow(p4040_m,3))))\n",
"sum_list.append(ztf.to_complex(p4160_s) * tf.exp(tf.complex(ztf.constant(0.0), p4160_p)) * ztf.to_complex(p4160_w / (tf.pow(p4160_m,3))))\n",
"sum_list.append(ztf.to_complex(p4415_s) * tf.exp(tf.complex(ztf.constant(0.0), p4415_p)) * ztf.to_complex(p4415_w / (tf.pow(p4415_m,3))))\n",
"sum_list.append(ztf.to_complex(DDstar_s) * tf.exp(tf.complex(ztf.constant(0.0), DDstar_p)) * (ztf.to_complex(1.0 / (10.0*tf.pow(Dstar_m,2)) + 1.0 / (10.0*tf.pow(D_m,2)))))\n",
"sum_list.append(ztf.to_complex(Dbar_s) * tf.exp(tf.complex(ztf.constant(0.0), Dbar_p)) * ztf.to_complex(1.0 / (6.0*tf.pow(Dbar_m,2))))\n",
"\n",
"sum_ru_1 = ztf.to_complex(ztf.constant(0.0))\n",
"\n",
"for part in sum_list:\n",
" sum_ru_1 += part\n",
"\n",
"sum_1 = tf.math.real(sum_ru_1)\n",
"# constraint1 = zfit.constraint.GaussianConstraint(params = sum_1, mu = ztf.constant(1.7*10**-8), \n",
"# sigma = ztf.constant(2.2*10**-8))\n",
"\n",
"constraint1 = tf.pow((sum_1-ztf.constant(1.7*10**-8))/ztf.constant(2.2*10**-8),2)/ztf.constant(2.)\n",
"\n",
"# 2. Constraint - Abs. of sum of Psi contribs and D contribs\n",
"\n",
"sum_2 = tf.abs(sum_ru_1)\n",
"constraint2 = tf.cond(tf.greater_equal(sum_2, 5.0e-8), lambda: 100000., lambda: 0.)\n",
"\n",
"# 3. Constraint - Maximum eta of D contribs\n",
"\n",
"constraint3_0 = tf.cond(tf.greater_equal(tf.abs(Dbar_s), 0.2), lambda: 100000., lambda: 0.)\n",
"\n",
"constraint3_1 = tf.cond(tf.greater_equal(tf.abs(DDstar_s), 0.2), lambda: 100000., lambda: 0.)\n",
"\n",
"# 4. Constraint - Formfactor multivariant gaussian covariance fplus\n",
"\n",
"Cov_matrix = [[ztf.constant( 1.), ztf.constant( 0.45), ztf.constant( 0.19), ztf.constant(0.857), ztf.constant(0.598), ztf.constant(0.531), ztf.constant(0.752), ztf.constant(0.229), ztf.constant(0,117)],\n",
" [ztf.constant( 0.45), ztf.constant( 1.), ztf.constant(0.677), ztf.constant(0.708), ztf.constant(0.958), ztf.constant(0.927), ztf.constant(0.227), ztf.constant(0.443), ztf.constant(0.287)],\n",
" [ztf.constant( 0.19), ztf.constant(0.677), ztf.constant( 1.), ztf.constant(0.595), ztf.constant(0.770), ztf.constant(0.819),ztf.constant(-0.023), ztf.constant( 0.07), ztf.constant(0.196)],\n",
" [ztf.constant(0.857), ztf.constant(0.708), ztf.constant(0.595), ztf.constant( 1.), ztf.constant( 0.83), ztf.constant(0.766), ztf.constant(0.582), ztf.constant(0.237), ztf.constant(0.192)],\n",
" [ztf.constant(0.598), ztf.constant(0.958), ztf.constant(0.770), ztf.constant( 0.83), ztf.constant( 1.), ztf.constant(0.973), ztf.constant(0.324), ztf.constant(0.372), ztf.constant(0.272)],\n",
" [ztf.constant(0.531), ztf.constant(0.927), ztf.constant(0.819), ztf.constant(0.766), ztf.constant(0.973), ztf.constant( 1.), ztf.constant(0.268), ztf.constant(0.332), ztf.constant(0.269)],\n",
" [ztf.constant(0.752), ztf.constant(0.227),ztf.constant(-0.023), ztf.constant(0.582), ztf.constant(0.324), ztf.constant(0.268), ztf.constant( 1.), ztf.constant( 0.59), ztf.constant(0.515)],\n",
" [ztf.constant(0.229), ztf.constant(0.443), ztf.constant( 0.07), ztf.constant(0.237), ztf.constant(0.372), ztf.constant(0.332), ztf.constant( 0.59), ztf.constant( 1.), ztf.constant(0.897)],\n",
" [ztf.constant(0.117), ztf.constant(0.287), ztf.constant(0.196), ztf.constant(0.192), ztf.constant(0.272), ztf.constant(0.269), ztf.constant(0.515), ztf.constant(0.897), ztf.constant( 1.)]]\n",
"\n",
"def triGauss(val1,val2,val3,m = Cov_matrix):\n",
"\n",
" mean1 = ztf.constant(0.466)\n",
" mean2 = ztf.constant(-0.885)\n",
" mean3 = ztf.constant(-0.213)\n",
" sigma1 = ztf.constant(0.014/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": 38,
"metadata": {},
"outputs": [],
"source": [
"def reset_param_values(): \n",
" jpsi_m.set_value(jpsi_mass)\n",
" jpsi_s.set_value(jpsi_scale)\n",
" jpsi_p.set_value(jpsi_phase)\n",
" jpsi_w.set_value(jpsi_width)\n",
" psi2s_m.set_value(psi2s_mass)\n",
" psi2s_s.set_value(psi2s_scale)\n",
" psi2s_p.set_value(psi2s_phase)\n",
" psi2s_w.set_value(psi2s_width)\n",
" p3770_m.set_value(p3770_mass)\n",
" p3770_s.set_value(p3770_scale)\n",
" p3770_p.set_value(p3770_phase)\n",
" p3770_w.set_value(p3770_width)\n",
" p4040_m.set_value(p4040_mass)\n",
" p4040_s.set_value(p4040_scale)\n",
" p4040_p.set_value(p4040_phase)\n",
" p4040_w.set_value(p4040_width)\n",
" p4160_m.set_value(p4160_mass)\n",
" p4160_s.set_value(p4160_scale)\n",
" p4160_p.set_value(p4160_phase)\n",
" p4160_w.set_value(p4160_width)\n",
" p4415_m.set_value(p4415_mass)\n",
" p4415_s.set_value(p4415_scale)\n",
" p4415_p.set_value(p4415_phase)\n",
" p4415_w.set_value(p4415_width)\n",
" rho_m.set_value(rho_mass)\n",
" rho_s.set_value(rho_scale)\n",
" rho_p.set_value(rho_phase)\n",
" rho_w.set_value(rho_width)\n",
" omega_m.set_value(omega_mass)\n",
" omega_s.set_value(omega_scale)\n",
" omega_p.set_value(omega_phase)\n",
" omega_w.set_value(omega_width)\n",
" phi_m.set_value(phi_mass)\n",
" phi_s.set_value(phi_scale)\n",
" phi_p.set_value(phi_phase)\n",
" phi_w.set_value(phi_width)\n",
" Dstar_m.set_value(Dstar_mass)\n",
" DDstar_s.set_value(0.0)\n",
" DDstar_p.set_value(0.0)\n",
" D_m.set_value(D_mass)\n",
" Dbar_m.set_value(Dbar_mass)\n",
" Dbar_s.set_value(0.0)\n",
" Dbar_p.set_value(0.0)\n",
" tau_m.set_value(pdg['tau_M'])\n",
" Ctt.set_value(0.0)\n",
" b0_0.set_value(0.292)\n",
" b0_1.set_value(0.281)\n",
" b0_2.set_value(0.150)\n",
" bplus_0.set_value(0.466)\n",
" bplus_1.set_value(-0.885)\n",
" bplus_2.set_value(-0.213)\n",
" bT_0.set_value(0.460)\n",
" bT_1.set_value(-1.089)\n",
" bT_2.set_value(-1.114)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Analysis"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"WARNING:tensorflow:From C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\core\\sample.py:163: to_int32 (from tensorflow.python.ops.math_ops) is deprecated and will be removed in a future version.\n",
"Instructions for updating:\n",
"Use tf.cast instead.\n",
"Toy 0: Generating data...\n",
"Toy 0: Data generation finished\n",
"Toy 0: Loading data...\n",
"Toy 0: Loading data finished\n",
"Toy 0: Fitting pdf...\n",
"------------------------------------------------------------------\n",
"| FCN = 2.979E+05 | Ncalls=616 (627 total) |\n",
"| EDM = 1.12E+06 (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: 297878.03533313115\n",
"----------------------------------------------------------------------------------------------\n",
"| | Name | Value | Hesse Err | Minos Err- | Minos Err+ | Limit- | Limit+ | Fixed |\n",
"----------------------------------------------------------------------------------------------\n",
"| 0 | Ctt | -0.017 | 0.133 | | | -1.5 | 1.5 | |\n",
"| 1 | omega_p | 4 | 11 | | |-6.28319 | 6.28319 | |\n",
"| 2 | phi_s | 19.27 | 0.29 | | | 14.8182 | 23.5818 | |\n",
"| 3 | p4040_s | 1.780 | 0.021 | | |0.00501244| 2.01499 | |\n",
"| 4 | rho_p | 5.93 | 0.24 | | |-6.28319 | 6.28319 | |\n",
"| 5 | Dbar_p | -4.8 | 0.4 | | |-6.28319 | 6.28319 | |\n",
"| 6 | p3770_s | 2.54 | 0.05 | | |0.918861 | 4.08114 | |\n",
"| 7 | phi_p | 4.45 | 0.17 | | |-6.28319 | 6.28319 | |\n",
"| 8 | psi2s_p | -2.268 | 0.017 | | |-6.28319 | 6.28319 | |\n",
"| 9 | DDstar_s | 0.258 | 0.008 | | | -0.3 | 0.3 | |\n",
"| 10| bplus_2 | -0.233 | 0.016 | | | -2 | 2 | |\n",
"| 11| Dbar_s | -2.996E-1 | 0.005E-1 | | | -0.3 | 0.3 | |\n",
"| 12| DDstar_p | -0.41 | 0.20 | | |-6.28319 | 6.28319 | |\n",
"| 13| p4415_p | -3.57 | 0.17 | | |-6.28319 | 6.28319 | |\n",
"| 14| omega_s | 8.64 | 0.17 | | | 4.19232 | 9.40768 | |\n",
"| 15| bplus_1 | -0.647 | 0.003 | | | -2 | 2 | |\n",
"| 16| p4040_p | 2.94 | 0.09 | | |-6.28319 | 6.28319 | |\n",
"| 17| p4160_s | 1.04 | 0.03 | | | 0.71676 | 3.68324 | |\n",
"| 18| p3770_p | 3.72 | 0.05 | | |-6.28319 | 6.28319 | |\n",
"| 19| jpsi_p | 1.365 | 0.028 | | |-6.28319 | 6.28319 | |\n",
"| 20| bplus_0 | 0.321 | 0.014 | | | -2 | 2 | |\n",
"| 21| rho_s | 0.5 | 0.3 | | |0.0253049| 2.0747 | |\n",
"| 22| p4415_s | 0.515 | 0.026 | | |0.126447 | 2.35355 | |\n",
"| 23| p4160_p | 3.95 | 0.07 | | |-6.28319 | 6.28319 | |\n",
"----------------------------------------------------------------------------------------------\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"| | Ctt omega_p phi_s p4040_s rho_p Dbar_p p3770_s phi_p psi2s_p DDstar_s bplus_2 Dbar_s DDstar_p p4415_p omega_s bplus_1 p4040_p p4160_s p3770_p jpsi_p bplus_0 rho_s p4415_s p4160_p |\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"| Ctt | 1.000 0.965 0.862 0.207 0.933 0.950 0.074 0.882 0.631 0.797 -0.593 -0.035 0.707 0.703 -0.904 -0.323 0.781 0.492 0.179 0.878 0.961 -0.962 -0.116 0.384 |\n",
"| omega_p | 0.965 1.000 0.894 0.210 0.967 0.988 0.083 0.914 0.654 0.826 -0.638 -0.036 0.723 0.723 -0.937 -0.337 0.815 0.506 0.200 0.913 0.996 -0.997 -0.123 0.408 |\n",
"| phi_s | 0.862 0.894 1.000 0.187 0.864 0.883 0.074 0.817 0.585 0.738 -0.576 -0.032 0.646 0.646 -0.837 -0.305 0.728 0.452 0.179 0.816 0.890 -0.891 -0.110 0.364 |\n",
"| p4040_s | 0.207 0.210 0.187 1.000 0.203 0.207 0.024 0.192 0.135 0.171 -0.137 -0.008 0.151 0.147 -0.196 -0.056 0.155 0.095 0.050 0.192 0.210 -0.209 -0.023 0.070 |\n",
"| rho_p | 0.933 0.967 0.864 0.203 1.000 0.955 0.080 0.884 0.632 0.799 -0.615 -0.035 0.699 0.700 -0.906 -0.324 0.788 0.489 0.193 0.883 0.964 -0.964 -0.119 0.394 |\n",
"| Dbar_p | 0.950 0.988 0.883 0.207 0.955 1.000 0.077 0.903 0.657 0.813 -0.617 -0.036 0.706 0.713 -0.925 -0.313 0.798 0.498 0.195 0.898 0.986 -0.985 -0.120 0.397 |\n",
"| p3770_s | 0.074 0.083 0.074 0.024 0.080 0.077 1.000 0.076 0.077 0.068 -0.039 -0.004 0.048 0.064 -0.078 -0.021 0.072 0.047 0.013 0.070 0.083 -0.083 -0.010 0.030 |\n",
"| phi_p | 0.882 0.914 0.817 0.192 0.884 0.903 0.076 1.000 0.598 0.755 -0.582 -0.033 0.662 0.662 -0.857 -0.306 0.745 0.463 0.183 0.835 0.912 -0.912 -0.112 0.373 |\n",
"| psi2s_p | 0.631 0.654 0.585 0.135 0.632 0.657 0.077 0.598 1.000 0.544 -0.418 -0.020 0.542 0.472 -0.613 -0.238 0.547 0.331 0.172 0.597 0.650 -0.652 -0.085 0.287 |\n",
"| DDstar_s | 0.797 0.826 0.738 0.171 0.799 0.813 0.068 0.755 0.544 1.000 -0.521 -0.030 0.593 0.593 -0.774 -0.258 0.669 0.414 0.168 0.751 0.825 -0.823 -0.101 0.336 |\n",
"| bplus_2 | -0.593 -0.638 -0.576 -0.137 -0.615 -0.617 -0.039 -0.582 -0.418 -0.521 1.000 0.022 -0.482 -0.467 0.596 0.134 -0.494 -0.318 -0.111 -0.567 -0.643 0.637 0.066 -0.231 |\n",
"| Dbar_s | -0.035 -0.036 -0.032 -0.008 -0.035 -0.036 -0.004 -0.033 -0.020 -0.030 0.022 1.000 -0.032 -0.026 0.034 0.013 -0.029 -0.019 -0.006 -0.032 -0.036 0.036 0.004 -0.014 |\n",
"| DDstar_p | 0.707 0.723 0.646 0.151 0.699 0.706 0.048 0.662 0.542 0.593 -0.482 -0.032 1.000 0.512 -0.677 -0.200 0.610 0.356 0.170 0.700 0.723 -0.721 -0.094 0.340 |\n",
"| p4415_p | 0.703 0.723 0.646 0.147 0.700 0.713 0.064 0.662 0.472 0.593 -0.467 -0.026 0.512 1.000 -0.678 -0.210 0.598 0.366 0.152 0.657 0.724 -0.721 -0.053 0.318 |\n",
"| omega_s | -0.904 -0.937 -0.837 -0.196 -0.906 -0.925 -0.078 -0.857 -0.613 -0.774 0.596 0.034 -0.677 -0.678 1.000 0.315 -0.763 -0.474 -0.187 -0.856 -0.934 0.934 0.115 -0.382 |\n",
"| bplus_1 | -0.323 -0.337 -0.305 -0.056 -0.324 -0.313 -0.021 -0.306 -0.238 -0.258 0.134 0.013 -0.200 -0.210 0.315 1.000 -0.250 -0.146 -0.081 -0.280 -0.352 0.336 0.042 -0.139 |\n",
"| p4040_p | 0.781 0.815 0.728 0.155 0.788 0.798 0.072 0.745 0.547 0.669 -0.494 -0.029 0.610 0.598 -0.763 -0.250 1.000 0.424 0.193 0.738 0.814 -0.812 -0.096 0.292 |\n",
"| p4160_s | 0.492 0.506 0.452 0.095 0.489 0.498 0.047 0.463 0.331 0.414 -0.318 -0.019 0.356 0.366 -0.474 -0.146 0.424 1.000 0.112 0.459 0.506 -0.504 -0.054 0.175 |\n",
"| p3770_p | 0.179 0.200 0.179 0.050 0.193 0.195 0.013 0.183 0.172 0.168 -0.111 -0.006 0.170 0.152 -0.187 -0.081 0.193 0.112 1.000 0.175 0.198 -0.199 -0.028 0.096 |\n",
"| jpsi_p | 0.878 0.913 0.816 0.192 0.883 0.898 0.070 0.835 0.597 0.751 -0.567 -0.032 0.700 0.657 -0.856 -0.280 0.738 0.459 0.175 1.000 0.912 -0.910 -0.112 0.372 |\n",
"| bplus_0 | 0.961 0.996 0.890 0.210 0.964 0.986 0.083 0.912 0.650 0.825 -0.643 -0.036 0.723 0.724 -0.934 -0.352 0.814 0.506 0.198 0.912 1.000 -0.993 -0.122 0.406 |\n",
"| rho_s | -0.962 -0.997 -0.891 -0.209 -0.964 -0.985 -0.083 -0.912 -0.652 -0.823 0.637 0.036 -0.721 -0.721 0.934 0.336 -0.812 -0.504 -0.199 -0.910 -0.993 1.000 0.122 -0.407 |\n",
"| p4415_s | -0.116 -0.123 -0.110 -0.023 -0.119 -0.120 -0.010 -0.112 -0.085 -0.101 0.066 0.004 -0.094 -0.053 0.115 0.042 -0.096 -0.054 -0.028 -0.112 -0.122 0.122 1.000 -0.050 |\n",
"| p4160_p | 0.384 0.408 0.364 0.070 0.394 0.397 0.030 0.373 0.287 0.336 -0.231 -0.014 0.340 0.318 -0.382 -0.139 0.292 0.175 0.096 0.372 0.406 -0.407 -0.050 1.000 |\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"Hesse errors: OrderedDict([(<zfit.Parameter 'Ctt' floating=True>, {'error': 0.13307582874104307}), (<zfit.Parameter 'omega_p' floating=True>, {'error': 11.312544525700241}), (<zfit.Parameter 'phi_s' floating=True>, {'error': 0.29435979835990267}), (<zfit.Parameter 'p4040_s' floating=True>, {'error': 0.020773255462037232}), (<zfit.Parameter 'rho_p' floating=True>, {'error': 0.2436242306495524}), (<zfit.Parameter 'Dbar_p' floating=True>, {'error': 0.42032440023747597}), (<zfit.Parameter 'p3770_s' floating=True>, {'error': 0.04542769991233664}), (<zfit.Parameter 'phi_p' floating=True>, {'error': 0.16705000022500105}), (<zfit.Parameter 'psi2s_p' floating=True>, {'error': 0.01701667878468638}), (<zfit.Parameter 'DDstar_s' floating=True>, {'error': 0.008402722883520175}), (<zfit.Parameter 'bplus_2' floating=True>, {'error': 0.015655568513703977}), (<zfit.Parameter 'Dbar_s' floating=True>, {'error': 0.00046055330791114346}), (<zfit.Parameter 'DDstar_p' floating=True>, {'error': 0.20138193109194136}), (<zfit.Parameter 'p4415_p' floating=True>, {'error': 0.17099716522674058}), (<zfit.Parameter 'omega_s' floating=True>, {'error': 0.16848322966260731}), (<zfit.Parameter 'bplus_1' floating=True>, {'error': 0.0025357515987360912}), (<zfit.Parameter 'p4040_p' floating=True>, {'error': 0.08785001542580062}), (<zfit.Parameter 'p4160_s' floating=True>, {'error': 0.03182752693261037}), (<zfit.Parameter 'p3770_p' floating=True>, {'error': 0.04623833317689208}), (<zfit.Parameter 'jpsi_p' floating=True>, {'error': 0.02804294337792035}), (<zfit.Parameter 'bplus_0' floating=True>, {'error': 0.01433694235612526}), (<zfit.Parameter 'rho_s' floating=True>, {'error': 0.34127837859046173}), (<zfit.Parameter 'p4415_s' floating=True>, {'error': 0.02591065046357091}), (<zfit.Parameter 'p4160_p' floating=True>, {'error': 0.0690782454106591})])\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\ipykernel_launcher.py:196: UserWarning: Creating legend with loc=\"best\" can be slow with large amounts of data.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Toy 1/10\n",
"Time taken: 1 min, 40 s\n",
"Projected time left: 15 min\n",
"Toy 1: Generating data...\n",
"Toy 1: Data generation finished\n",
"Toy 1: Loading data...\n",
"Toy 1: Loading data finished\n",
"Toy 1: Fitting pdf...\n",
"------------------------------------------------------------------\n",
"| FCN = 2.979E+05 | Ncalls=1287 (1287 total) |\n",
"| EDM = 6.42E-05 (Goal: 5E-06) | up = 0.5 |\n",
"------------------------------------------------------------------\n",
"| Valid Min. | Valid Param. | Above EDM | Reached call limit |\n",
"------------------------------------------------------------------\n",
"| True | True | False | False |\n",
"------------------------------------------------------------------\n",
"| Hesse failed | Has cov. | Accurate | Pos. def. | Forced |\n",
"------------------------------------------------------------------\n",
"| False | True | False | False | True |\n",
"------------------------------------------------------------------\n",
"Function minimum: 297909.6307224464\n",
"----------------------------------------------------------------------------------------------\n",
"| | Name | Value | Hesse Err | Minos Err- | Minos Err+ | Limit- | Limit+ | Fixed |\n",
"----------------------------------------------------------------------------------------------\n",
"| 0 | Ctt | 0.62 | 0.19 | | | -1.5 | 1.5 | |\n",
"| 1 | omega_p | -5.6 | 0.4 | | |-6.28319 | 6.28319 | |\n",
"| 2 | phi_s | 22.0 | 1.0 | | | 14.8182 | 23.5818 | |\n",
"| 3 | p4040_s | 1.09 | 0.18 | | |0.00501244| 2.01499 | |\n",
"| 4 | rho_p | -0.30 | 0.71 | | |-6.28319 | 6.28319 | |\n",
"| 5 | Dbar_p | 0.9 | 0.3 | | |-6.28319 | 6.28319 | |\n",
"| 6 | p3770_s | 1.89 | 0.24 | | |0.918861 | 4.08114 | |\n",
"| 7 | phi_p | -5.53 | 0.15 | | |-6.28319 | 6.28319 | |\n",
"| 8 | psi2s_p | -2.08 | 0.03 | | |-6.28319 | 6.28319 | |\n",
"| 9 | DDstar_s | 0.300 | 0.024 | | | -0.3 | 0.3 | |\n",
"| 10| bplus_2 | 0.12 | 0.07 | | | -2 | 2 | |\n",
"| 11| Dbar_s | -0.300 | 0.015 | | | -0.3 | 0.3 | |\n",
"| 12| DDstar_p | 4.0 | 0.4 | | |-6.28319 | 6.28319 | |\n",
"| 13| p4415_p | -2.54 | 0.13 | | |-6.28319 | 6.28319 | |\n",
"| 14| omega_s | 8 | 4 | | | 4.19232 | 9.40768 | |\n",
"| 15| bplus_1 | 0.73 | 0.06 | | | -2 | 2 | |\n",
"| 16| p4040_p | -2.70 | 0.17 | | |-6.28319 | 6.28319 | |\n",
"| 17| p4160_s | 2.60 | 0.17 | | | 0.71676 | 3.68324 | |\n",
"| 18| p3770_p | 4.24 | 0.14 | | |-6.28319 | 6.28319 | |\n",
"| 19| jpsi_p | 4.580 | 0.025 | | |-6.28319 | 6.28319 | |\n",
"| 20| bplus_0 | -0.376 | 0.028 | | | -2 | 2 | |\n",
"| 21| rho_s | 0.6 | 0.4 | | |0.0253049| 2.0747 | |\n",
"| 22| p4415_s | 1.73 | 0.19 | | |0.126447 | 2.35355 | |\n",
"| 23| p4160_p | -2.13 | 0.09 | | |-6.28319 | 6.28319 | |\n",
"----------------------------------------------------------------------------------------------\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"| | Ctt omega_p phi_s p4040_s rho_p Dbar_p p3770_s phi_p psi2s_p DDstar_s bplus_2 Dbar_s DDstar_p p4415_p omega_s bplus_1 p4040_p p4160_s p3770_p jpsi_p bplus_0 rho_s p4415_s p4160_p |\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"| Ctt | 1.000 -0.014 -0.053 0.288 0.090 -0.352 -0.037 -0.015 -0.436 0.012 -0.639 0.002 -0.141 0.135 -0.012 0.102 -0.020 0.228 -0.274 0.467 0.026 -0.005 0.075 -0.120 |\n",
"| omega_p | -0.014 1.000 -0.012 -0.000 -0.141 0.005 0.004 -0.005 -0.007 0.002 0.022 0.001 0.018 -0.001 0.827 0.177 0.004 0.001 -0.006 -0.041 -0.198 0.419 0.001 -0.006 |\n",
"| phi_s | -0.053 -0.012 1.000 -0.006 -0.108 -0.014 -0.032 0.457 0.032 -0.001 0.158 -0.003 -0.007 -0.041 0.012 -0.146 -0.036 -0.039 0.016 0.056 0.161 -0.124 0.020 -0.017 |\n",
"| p4040_s | 0.288 -0.000 -0.006 1.000 0.005 0.025 0.156 -0.001 -0.200 -0.000 -0.001 0.011 -0.134 -0.143 -0.002 -0.021 -0.166 -0.070 -0.159 0.066 0.004 0.004 0.139 -0.417 |\n",
"| rho_p | 0.090 -0.141 -0.108 0.005 1.000 0.007 0.037 -0.294 -0.044 0.000 -0.253 0.003 -0.020 0.046 0.019 -0.216 0.031 0.047 -0.012 -0.049 0.235 -0.045 -0.027 0.024 |\n",
"| Dbar_p | -0.352 0.005 -0.014 0.025 0.007 1.000 -0.057 -0.005 0.023 -0.006 0.026 0.012 -0.164 -0.152 -0.002 -0.074 -0.174 -0.002 -0.061 -0.232 0.002 0.023 0.058 -0.220 |\n",
"| p3770_s | -0.037 0.004 -0.032 0.156 0.037 -0.057 1.000 -0.011 0.265 -0.004 -0.000 -0.016 0.076 0.005 -0.002 -0.020 0.033 0.090 -0.131 -0.082 -0.022 0.020 0.070 -0.048 |\n",
"| phi_p | -0.015 -0.005 0.457 -0.001 -0.294 -0.005 -0.011 1.000 0.009 0.001 0.063 -0.000 -0.001 -0.010 -0.081 0.053 -0.009 -0.012 0.002 0.015 -0.052 0.036 0.006 -0.006 |\n",
"| psi2s_p | -0.436 -0.007 0.032 -0.200 -0.044 0.023 0.265 0.009 1.000 0.019 0.156 0.026 -0.044 0.091 0.001 0.019 0.183 -0.094 0.238 -0.094 0.006 -0.026 -0.136 0.264 |\n",
"| DDstar_s | 0.012 0.002 -0.001 -0.000 0.000 -0.006 -0.004 0.001 0.019 1.000 0.016 -0.001 0.032 -0.000 0.000 -0.032 -0.003 0.000 0.013 0.037 0.006 0.006 0.010 0.012 |\n",
"| bplus_2 | -0.639 0.022 0.158 -0.001 -0.253 0.026 -0.000 0.063 0.156 0.016 1.000 -0.006 0.262 -0.234 0.023 -0.231 -0.096 -0.080 0.082 -0.264 0.032 -0.026 0.231 -0.163 |\n",
"| Dbar_s | 0.002 0.001 -0.003 0.011 0.003 0.012 -0.016 -0.000 0.026 -0.001 -0.006 1.000 -0.002 -0.001 -0.000 -0.011 -0.008 0.004 0.004 0.030 0.002 0.003 0.011 -0.003 |\n",
"| DDstar_p | -0.141 0.018 -0.007 -0.134 -0.020 -0.164 0.076 -0.001 -0.044 0.032 0.262 -0.002 1.000 -0.347 -0.001 -0.271 -0.165 -0.078 0.079 -0.296 0.030 0.056 0.049 -0.240 |\n",
"| p4415_p | 0.135 -0.001 -0.041 -0.143 0.046 -0.152 0.005 -0.010 0.091 -0.000 -0.234 -0.001 -0.347 1.000 -0.006 0.169 0.188 -0.050 0.066 0.149 -0.069 0.012 -0.183 0.383 |\n",
"| omega_s | -0.012 0.827 0.012 -0.002 0.019 -0.002 -0.002 -0.081 0.001 0.000 0.023 -0.000 -0.001 -0.006 1.000 -0.000 -0.004 -0.005 0.001 -0.017 0.000 0.070 0.002 -0.003 |\n",
"| bplus_1 | 0.102 0.177 -0.146 -0.021 -0.216 -0.074 -0.020 0.053 0.019 -0.032 -0.231 -0.011 -0.271 0.169 -0.000 1.000 0.114 -0.027 -0.016 0.047 -0.942 0.437 -0.147 0.172 |\n",
"| p4040_p | -0.020 0.004 -0.036 -0.166 0.031 -0.174 0.033 -0.009 0.183 -0.003 -0.096 -0.008 -0.165 0.188 -0.004 0.114 1.000 0.372 0.193 0.027 -0.072 0.022 0.074 0.090 |\n",
"| p4160_s | 0.228 0.001 -0.039 -0.070 0.047 -0.002 0.090 -0.012 -0.094 0.000 -0.080 0.004 -0.078 -0.050 -0.005 -0.027 0.372 1.000 -0.027 0.019 -0.011 0.020 0.330 -0.147 |\n",
"| p3770_p | -0.274 -0.006 0.016 -0.159 -0.012 -0.061 -0.131 0.002 0.238 0.013 0.082 0.004 0.079 0.066 0.001 -0.016 0.193 -0.027 1.000 -0.069 0.020 -0.020 -0.089 0.192 |\n",
"| jpsi_p | 0.467 -0.041 0.056 0.066 -0.049 -0.232 -0.082 0.015 -0.094 0.037 -0.264 0.030 -0.296 0.149 -0.017 0.047 0.027 0.019 -0.069 1.000 0.124 -0.089 -0.060 0.092 |\n",
"| bplus_0 | 0.026 -0.198 0.161 0.004 0.235 0.002 -0.022 -0.052 0.006 0.006 0.032 0.002 0.030 -0.069 0.000 -0.942 -0.072 -0.011 0.020 0.124 1.000 -0.496 0.044 -0.050 |\n",
"| rho_s | -0.005 0.419 -0.124 0.004 -0.045 0.023 0.020 0.036 -0.026 0.006 -0.026 0.003 0.056 0.012 0.070 0.437 0.022 0.020 -0.020 -0.089 -0.496 1.000 -0.004 -0.010 |\n",
"| p4415_s | 0.075 0.001 0.020 0.139 -0.027 0.058 0.070 0.006 -0.136 0.010 0.231 0.011 0.049 -0.183 0.002 -0.147 0.074 0.330 -0.089 -0.060 0.044 -0.004 1.000 -0.107 |\n",
"| p4160_p | -0.120 -0.006 -0.017 -0.417 0.024 -0.220 -0.048 -0.006 0.264 0.012 -0.163 -0.003 -0.240 0.383 -0.003 0.172 0.090 -0.147 0.192 0.092 -0.050 -0.010 -0.107 1.000 |\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"Hesse errors: OrderedDict([(<zfit.Parameter 'Ctt' floating=True>, {'error': 0.18748768755829082}), (<zfit.Parameter 'omega_p' floating=True>, {'error': 0.40323634554077215}), (<zfit.Parameter 'phi_s' floating=True>, {'error': 0.985024000483012}), (<zfit.Parameter 'p4040_s' floating=True>, {'error': 0.17586720794456895}), (<zfit.Parameter 'rho_p' floating=True>, {'error': 0.7105697644993554}), (<zfit.Parameter 'Dbar_p' floating=True>, {'error': 0.32118284401931696}), (<zfit.Parameter 'p3770_s' floating=True>, {'error': 0.2372102719369873}), (<zfit.Parameter 'phi_p' floating=True>, {'error': 0.14503339498061774}), (<zfit.Parameter 'psi2s_p' floating=True>, {'error': 0.03185662332362904}), (<zfit.Parameter 'DDstar_s' floating=True>, {'error': 0.02421642931361845}), (<zfit.Parameter 'bplus_2' floating=True>, {'error': 0.07483773848310715}), (<zfit.Parameter 'Dbar_s' floating=True>, {'error': 0.01520721320902993}), (<zfit.Parameter 'DDstar_p' floating=True>, {'error': 0.39795186227992385}), (<zfit.Parameter 'p4415_p' floating=True>, {'error': 0.13390670241756997}), (<zfit.Parameter 'omega_s' floating=True>, {'error': 3.8706357513686984}), (<zfit.Parameter 'bplus_1' floating=True>, {'error': 0.057866643089661274}), (<zfit.Parameter 'p4040_p' floating=True>, {'error': 0.16823464853939285}), (<zfit.Parameter 'p4160_s' floating=True>, {'error': 0.171806968358724}), (<zfit.Parameter 'p3770_p' floating=True>, {'error': 0.138889598592427}), (<zfit.Parameter 'jpsi_p' floating=True>, {'error': 0.02515049586844853}), (<zfit.Parameter 'bplus_0' floating=True>, {'error': 0.02826194125557635}), (<zfit.Parameter 'rho_s' floating=True>, {'error': 0.39800274915520395}), (<zfit.Parameter 'p4415_s' floating=True>, {'error': 0.19429670305431568}), (<zfit.Parameter 'p4160_p' floating=True>, {'error': 0.08939828920188164})])\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Toy 2/10\n",
"Time taken: 3 min, 54 s\n",
"Projected time left: 15 min, 36 s\n",
"Toy 2: Generating data...\n",
"Toy 2: Data generation finished\n",
"Toy 2: Loading data...\n",
"Toy 2: Loading data finished\n",
"Toy 2: Fitting pdf...\n",
"------------------------------------------------------------------\n",
"| FCN = 2.977E+05 | Ncalls=1481 (1481 total) |\n",
"| EDM = 5.87E-05 (Goal: 5E-06) | up = 0.5 |\n",
"------------------------------------------------------------------\n",
"| Valid Min. | Valid Param. | Above EDM | Reached call limit |\n",
"------------------------------------------------------------------\n",
"| True | True | False | False |\n",
"------------------------------------------------------------------\n",
"| Hesse failed | Has cov. | Accurate | Pos. def. | Forced |\n",
"------------------------------------------------------------------\n",
"| False | True | True | True | False |\n",
"------------------------------------------------------------------\n",
"Function minimum: 297742.96509176795\n",
"----------------------------------------------------------------------------------------------\n",
"| | Name | Value | Hesse Err | Minos Err- | Minos Err+ | Limit- | Limit+ | Fixed |\n",
"----------------------------------------------------------------------------------------------\n",
"| 0 | Ctt | -0.017 | 0.193 | | | -1.5 | 1.5 | |\n",
"| 1 | omega_p | 6.28 | 0.09 | | |-6.28319 | 6.28319 | |\n",
"| 2 | phi_s | 23.0 | 0.9 | | | 14.8182 | 23.5818 | |\n",
"| 3 | p4040_s | 1.56 | 0.17 | | |0.00501244| 2.01499 | |\n",
"| 4 | rho_p | -0.55 | 0.27 | | |-6.28319 | 6.28319 | |\n",
"| 5 | Dbar_p | 1.5 | 0.6 | | |-6.28319 | 6.28319 | |\n",
"| 6 | p3770_s | 0.919 | 0.024 | | |0.918861 | 4.08114 | |\n",
"| 7 | phi_p | -5.27 | 0.12 | | |-6.28319 | 6.28319 | |\n",
"| 8 | psi2s_p | -2.04 | 0.03 | | |-6.28319 | 6.28319 | |\n",
"| 9 | DDstar_s | 0.30 | 0.05 | | | -0.3 | 0.3 | |\n",
"| 10| bplus_2 | 0.23 | 0.10 | | | -2 | 2 | |\n",
"| 11| Dbar_s | -0.30 | 0.08 | | | -0.3 | 0.3 | |\n",
"| 12| DDstar_p | -0.011 | 0.472 | | |-6.28319 | 6.28319 | |\n",
"| 13| p4415_p | 4.4 | 0.5 | | |-6.28319 | 6.28319 | |\n",
"| 14| omega_s | 6.0 | 0.9 | | | 4.19232 | 9.40768 | |\n",
"| 15| bplus_1 | 1.02 | 0.08 | | | -2 | 2 | |\n",
"| 16| p4040_p | 2.86 | 0.11 | | |-6.28319 | 6.28319 | |\n",
"| 17| p4160_s | 0.717 | 0.012 | | | 0.71676 | 3.68324 | |\n",
"| 18| p3770_p | -6.283 | 0.019 | | |-6.28319 | 6.28319 | |\n",
"| 19| jpsi_p | 4.613 | 0.030 | | |-6.28319 | 6.28319 | |\n",
"| 20| bplus_0 | -0.55 | 0.04 | | | -2 | 2 | |\n",
"| 21| rho_s | 1.37 | 0.31 | | |0.0253049| 2.0747 | |\n",
"| 22| p4415_s | 0.40 | 0.19 | | |0.126447 | 2.35355 | |\n",
"| 23| p4160_p | -6.283 | 0.011 | | |-6.28319 | 6.28319 | |\n",
"----------------------------------------------------------------------------------------------\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"| | Ctt omega_p phi_s p4040_s rho_p Dbar_p p3770_s phi_p psi2s_p DDstar_s bplus_2 Dbar_s DDstar_p p4415_p omega_s bplus_1 p4040_p p4160_s p3770_p jpsi_p bplus_0 rho_s p4415_s p4160_p |\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"| Ctt | 1.000 -0.001 -0.001 0.031 0.025 -0.498 -0.008 -0.001 -0.402 -0.011 -0.589 0.022 0.398 0.240 0.026 0.041 0.334 0.004 0.001 0.479 0.157 -0.079 0.022 -0.006 |\n",
"| omega_p | -0.001 1.000 -0.002 0.000 0.021 0.001 -0.000 0.005 0.000 -0.000 0.001 0.000 -0.001 -0.000 -0.081 0.003 -0.000 -0.000 0.000 0.000 -0.004 0.014 0.001 0.000 |\n",
"| phi_s | -0.001 -0.002 1.000 -0.024 -0.143 -0.082 0.000 0.331 -0.000 0.003 0.120 -0.020 0.003 -0.020 -0.004 -0.093 -0.005 0.001 -0.001 0.029 0.117 -0.034 -0.012 0.000 |\n",
"| p4040_s | 0.031 0.000 -0.024 1.000 0.037 0.373 -0.006 -0.007 -0.081 0.000 -0.015 0.075 -0.081 -0.213 -0.013 -0.067 -0.269 -0.007 0.005 0.017 -0.057 0.028 0.104 0.010 |\n",
"| rho_p | 0.025 0.021 -0.143 0.037 1.000 0.123 -0.000 -0.248 -0.013 -0.001 -0.253 0.034 0.007 0.041 0.159 0.008 0.008 -0.001 0.001 -0.051 -0.048 0.161 0.008 -0.000 |\n",
"| Dbar_p | -0.498 0.001 -0.082 0.373 0.123 1.000 0.002 -0.029 0.275 0.032 -0.054 0.078 -0.300 -0.262 -0.030 -0.082 -0.417 -0.014 0.013 -0.132 -0.192 0.070 0.184 0.007 |\n",
"| p3770_s | -0.008 -0.000 0.000 -0.006 -0.000 0.002 1.000 0.000 0.009 -0.001 0.001 0.002 0.004 0.002 0.000 -0.001 0.005 0.000 -0.002 0.002 0.001 -0.000 -0.002 -0.000 |\n",
"| phi_p | -0.001 0.005 0.331 -0.007 -0.248 -0.029 0.000 1.000 0.009 -0.001 0.073 -0.006 0.014 -0.011 -0.117 0.021 0.003 0.000 -0.000 0.023 -0.010 0.131 0.002 0.000 |\n",
"| psi2s_p | -0.402 0.000 -0.000 -0.081 -0.013 0.275 0.009 0.009 1.000 0.022 0.188 0.133 0.236 -0.040 -0.017 -0.020 0.042 0.006 -0.006 0.122 -0.073 0.041 -0.030 0.001 |\n",
"| DDstar_s | -0.011 -0.000 0.003 0.000 -0.001 0.032 -0.001 -0.001 0.022 1.000 -0.049 0.006 0.027 0.029 0.007 0.014 0.038 0.000 0.000 0.012 0.030 -0.017 0.003 -0.001 |\n",
"| bplus_2 | -0.589 0.001 0.120 -0.015 -0.253 -0.054 0.001 0.073 0.188 -0.049 1.000 -0.047 -0.071 -0.245 -0.049 0.058 -0.042 -0.003 -0.003 -0.221 -0.180 0.123 0.161 0.001 |\n",
"| Dbar_s | 0.022 0.000 -0.020 0.075 0.034 0.078 0.002 -0.006 0.133 0.006 -0.047 1.000 -0.092 -0.038 -0.003 -0.010 -0.000 -0.001 0.003 0.104 -0.030 0.006 0.055 0.001 |\n",
"| DDstar_p | 0.398 -0.001 0.003 -0.081 0.007 -0.300 0.004 0.014 0.236 0.027 -0.071 -0.092 1.000 0.139 0.006 0.015 0.483 0.015 0.000 0.635 0.033 -0.017 0.036 -0.005 |\n",
"| p4415_p | 0.240 -0.000 -0.020 -0.213 0.041 -0.262 0.002 -0.011 -0.040 0.029 -0.245 -0.038 0.139 1.000 0.013 0.068 0.076 0.011 -0.003 0.106 0.034 -0.035 -0.149 0.001 |\n",
"| omega_s | 0.026 -0.081 -0.004 -0.013 0.159 -0.030 0.000 -0.117 -0.017 0.007 -0.049 -0.003 0.006 0.013 1.000 -0.122 0.006 0.001 -0.000 -0.002 0.132 -0.460 -0.017 -0.000 |\n",
"| bplus_1 | 0.041 0.003 -0.093 -0.067 0.008 -0.082 -0.001 0.021 -0.020 0.014 0.058 -0.010 0.015 0.068 -0.122 1.000 0.065 0.002 -0.002 -0.041 -0.941 0.268 -0.050 -0.002 |\n",
"| p4040_p | 0.334 -0.000 -0.005 -0.269 0.008 -0.417 0.005 0.003 0.042 0.038 -0.042 -0.000 0.483 0.076 0.006 0.065 1.000 0.001 -0.009 0.265 0.011 -0.017 0.092 -0.004 |\n",
"| p4160_s | 0.004 -0.000 0.001 -0.007 -0.001 -0.014 0.000 0.000 0.006 0.000 -0.003 -0.001 0.015 0.011 0.001 0.002 0.001 1.000 -0.000 0.009 0.003 -0.001 -0.001 -0.002 |\n",
"| p3770_p | 0.001 0.000 -0.001 0.005 0.001 0.013 -0.002 -0.000 -0.006 0.000 -0.003 0.003 0.000 -0.003 -0.000 -0.002 -0.009 -0.000 1.000 0.004 -0.001 0.001 0.002 0.000 |\n",
"| jpsi_p | 0.479 0.000 0.029 0.017 -0.051 -0.132 0.002 0.023 0.122 0.012 -0.221 0.104 0.635 0.106 -0.002 -0.041 0.265 0.009 0.004 1.000 0.103 -0.017 0.019 -0.002 |\n",
"| bplus_0 | 0.157 -0.004 0.117 -0.057 -0.048 -0.192 0.001 -0.010 -0.073 0.030 -0.180 -0.030 0.033 0.034 0.132 -0.941 0.011 0.003 -0.001 0.103 1.000 -0.299 -0.080 0.000 |\n",
"| rho_s | -0.079 0.014 -0.034 0.028 0.161 0.070 -0.000 0.131 0.041 -0.017 0.123 0.006 -0.017 -0.035 -0.460 0.268 -0.017 -0.001 0.001 -0.017 -0.299 1.000 0.041 0.000 |\n",
"| p4415_s | 0.022 0.001 -0.012 0.104 0.008 0.184 -0.002 0.002 -0.030 0.003 0.161 0.055 0.036 -0.149 -0.017 -0.050 0.092 -0.001 0.002 0.019 -0.080 0.041 1.000 -0.005 |\n",
"| p4160_p | -0.006 0.000 0.000 0.010 -0.000 0.007 -0.000 0.000 0.001 -0.001 0.001 0.001 -0.005 0.001 -0.000 -0.002 -0.004 -0.002 0.000 -0.002 0.000 0.000 -0.005 1.000 |\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"Hesse errors: OrderedDict([(<zfit.Parameter 'Ctt' floating=True>, {'error': 0.19338376498028875}), (<zfit.Parameter 'omega_p' floating=True>, {'error': 0.09355253342183101}), (<zfit.Parameter 'phi_s' floating=True>, {'error': 0.8591335437952221}), (<zfit.Parameter 'p4040_s' floating=True>, {'error': 0.1745741995718021}), (<zfit.Parameter 'rho_p' floating=True>, {'error': 0.26975749374548696}), (<zfit.Parameter 'Dbar_p' floating=True>, {'error': 0.5675491117289342}), (<zfit.Parameter 'p3770_s' floating=True>, {'error': 0.023645611184699322}), (<zfit.Parameter 'phi_p' floating=True>, {'error': 0.11694285135085192}), (<zfit.Parameter 'psi2s_p' floating=True>, {'error': 0.034894129483006964}), (<zfit.Parameter 'DDstar_s' floating=True>, {'error': 0.049562852297974136}), (<zfit.Parameter 'bplus_2' floating=True>, {'error': 0.09992913474301979}), (<zfit.Parameter 'Dbar_s' floating=True>, {'error': 0.0804709270108857}), (<zfit.Parameter 'DDstar_p' floating=True>, {'error': 0.4721007315609622}), (<zfit.Parameter 'p4415_p' floating=True>, {'error': 0.4971009105610804}), (<zfit.Parameter 'omega_s' floating=True>, {'error': 0.9373573080054802}), (<zfit.Parameter 'bplus_1' floating=True>, {'error': 0.07869165370384468}), (<zfit.Parameter 'p4040_p' floating=True>, {'error': 0.11078527371686864}), (<zfit.Parameter 'p4160_s' floating=True>, {'error': 0.011981558188349617}), (<zfit.Parameter 'p3770_p' floating=True>, {'error': 0.01896914517553716}), (<zfit.Parameter 'jpsi_p' floating=True>, {'error': 0.029718397699545207}), (<zfit.Parameter 'bplus_0' floating=True>, {'error': 0.0430274851462954}), (<zfit.Parameter 'rho_s' floating=True>, {'error': 0.3078557703077741}), (<zfit.Parameter 'p4415_s' floating=True>, {'error': 0.1904672232500821}), (<zfit.Parameter 'p4160_p' floating=True>, {'error': 0.010943271603894189})])\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Toy 3/10\n",
"Time taken: 6 min, 17 s\n",
"Projected time left: 14 min, 35 s\n",
"Toy 3: Generating data...\n",
"Toy 3: Data generation finished\n",
"Toy 3: Loading data...\n",
"Toy 3: Loading data finished\n",
"Toy 3: Fitting pdf...\n",
"------------------------------------------------------------------\n",
"| FCN = 2.978E+05 | Ncalls=1302 (1302 total) |\n",
"| EDM = 3.73E-05 (Goal: 5E-06) | up = 0.5 |\n",
"------------------------------------------------------------------\n",
"| Valid Min. | Valid Param. | Above EDM | Reached call limit |\n",
"------------------------------------------------------------------\n",
"| True | True | False | False |\n",
"------------------------------------------------------------------\n",
"| Hesse failed | Has cov. | Accurate | Pos. def. | Forced |\n",
"------------------------------------------------------------------\n",
"| False | True | False | False | True |\n",
"------------------------------------------------------------------\n",
"Function minimum: 297798.8924556829\n",
"----------------------------------------------------------------------------------------------\n",
"| | Name | Value | Hesse Err | Minos Err- | Minos Err+ | Limit- | Limit+ | Fixed |\n",
"----------------------------------------------------------------------------------------------\n",
"| 0 | Ctt | 0.53 | 0.24 | | | -1.5 | 1.5 | |\n",
"| 1 | omega_p | -5.63 | 0.23 | | |-6.28319 | 6.28319 | |\n",
"| 2 | phi_s | 17.9 | 1.3 | | | 14.8182 | 23.5818 | |\n",
"| 3 | p4040_s | 0.98 | 0.17 | | |0.00501244| 2.01499 | |\n",
"| 4 | rho_p | -6.28 | 0.13 | | |-6.28319 | 6.28319 | |\n",
"| 5 | Dbar_p | -2.4 | 0.8 | | |-6.28319 | 6.28319 | |\n",
"| 6 | p3770_s | 0.919 | 0.023 | | |0.918861 | 4.08114 | |\n",
"| 7 | phi_p | 0.51 | 0.29 | | |-6.28319 | 6.28319 | |\n",
"| 8 | psi2s_p | -2.18 | 0.04 | | |-6.28319 | 6.28319 | |\n",
"| 9 | DDstar_s | -0.30 | 0.41 | | | -0.3 | 0.3 | |\n",
"| 10| bplus_2 | 0.04 | 0.10 | | | -2 | 2 | |\n",
"| 11| Dbar_s | -0.30 | 0.07 | | | -0.3 | 0.3 | |\n",
"| 12| DDstar_p | -1.4 | 0.6 | | |-6.28319 | 6.28319 | |\n",
"| 13| p4415_p | 4.10 | 0.21 | | |-6.28319 | 6.28319 | |\n",
"| 14| omega_s | 8.4 | 1.0 | | | 4.19232 | 9.40768 | |\n",
"| 15| bplus_1 | -0.95 | 0.06 | | | -2 | 2 | |\n",
"| 16| p4040_p | 4.07 | 0.20 | | |-6.28319 | 6.28319 | |\n",
"| 17| p4160_s | 2.34 | 0.17 | | | 0.71676 | 3.68324 | |\n",
"| 18| p3770_p | -6.283 | 0.018 | | |-6.28319 | 6.28319 | |\n",
"| 19| jpsi_p | -4.72 | 0.04 | | |-6.28319 | 6.28319 | |\n",
"| 20| bplus_0 | 0.502 | 0.023 | | | -2 | 2 | |\n",
"| 21| rho_s | 0.7 | 0.4 | | |0.0253049| 2.0747 | |\n",
"| 22| p4415_s | 1.11 | 0.20 | | |0.126447 | 2.35355 | |\n",
"| 23| p4160_p | 4.41 | 0.12 | | |-6.28319 | 6.28319 | |\n",
"----------------------------------------------------------------------------------------------\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"| | Ctt omega_p phi_s p4040_s rho_p Dbar_p p3770_s phi_p psi2s_p DDstar_s bplus_2 Dbar_s DDstar_p p4415_p omega_s bplus_1 p4040_p p4160_s p3770_p jpsi_p bplus_0 rho_s p4415_s p4160_p |\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"| Ctt | 1.000 0.013 -0.024 0.274 0.010 0.598 -0.001 0.010 -0.033 -0.039 0.439 0.141 -0.097 0.182 -0.018 -0.114 0.038 0.178 0.006 0.311 0.027 0.012 0.106 0.112 |\n",
"| omega_p | 0.013 1.000 0.036 0.003 -0.057 -0.013 0.001 0.031 0.015 -0.031 0.009 0.001 -0.077 -0.016 0.565 -0.011 -0.010 0.012 0.000 0.010 0.104 0.448 0.019 -0.009 |\n",
"| phi_s | -0.024 0.036 1.000 -0.012 -0.012 -0.009 0.000 0.803 0.009 -0.001 -0.076 0.004 -0.013 -0.025 0.017 0.015 -0.016 -0.020 0.000 -0.033 -0.011 -0.046 0.006 -0.012 |\n",
"| p4040_s | 0.274 0.003 -0.012 1.000 0.003 0.006 -0.005 -0.003 -0.222 -0.024 0.095 -0.040 0.105 -0.150 -0.008 -0.025 -0.188 0.021 -0.003 -0.036 0.021 0.009 0.176 -0.415 |\n",
"| rho_p | 0.010 -0.057 -0.012 0.003 1.000 0.006 -0.000 -0.022 -0.008 0.004 0.025 -0.002 0.014 0.011 0.011 0.008 0.006 0.004 -0.000 0.001 -0.029 -0.062 -0.006 0.005 |\n",
"| Dbar_p | 0.598 -0.013 -0.009 0.006 0.006 1.000 0.009 0.005 0.357 0.066 -0.079 0.181 -0.038 0.320 0.003 -0.135 0.384 -0.036 -0.002 0.702 -0.032 -0.055 -0.076 0.559 |\n",
"| p3770_s | -0.001 0.001 0.000 -0.005 -0.000 0.009 1.000 0.000 0.015 -0.006 -0.007 0.002 -0.014 0.003 0.000 0.007 0.010 -0.001 -0.002 0.011 0.001 0.002 -0.002 0.012 |\n",
"| phi_p | 0.010 0.031 0.803 -0.003 -0.022 0.005 0.000 1.000 0.003 -0.009 -0.011 0.001 -0.015 -0.004 -0.017 0.006 -0.003 -0.002 0.000 -0.012 -0.001 -0.028 0.001 -0.000 |\n",
"| psi2s_p | -0.033 0.015 0.009 -0.222 -0.008 0.357 0.015 0.003 1.000 -0.296 -0.322 0.002 -0.425 0.114 0.004 0.214 0.303 -0.122 -0.008 0.425 0.024 0.035 -0.103 0.457 |\n",
"| DDstar_s | -0.039 -0.031 -0.001 -0.024 0.004 0.066 -0.006 -0.009 -0.296 1.000 0.136 -0.086 0.324 -0.015 0.005 -0.226 -0.016 -0.014 -0.003 -0.280 -0.055 -0.068 -0.074 -0.142 |\n",
"| bplus_2 | 0.439 0.009 -0.076 0.095 0.025 -0.079 -0.007 -0.011 -0.322 0.136 1.000 -0.070 0.388 0.230 -0.046 -0.497 0.006 0.143 0.002 -0.263 0.100 0.013 -0.174 -0.046 |\n",
"| Dbar_s | 0.141 0.001 0.004 -0.040 -0.002 0.181 0.002 0.001 0.002 -0.086 -0.070 1.000 -0.176 0.055 0.002 0.034 0.099 -0.016 -0.002 0.119 -0.003 0.002 -0.035 0.135 |\n",
"| DDstar_p | -0.097 -0.077 -0.013 0.105 0.014 -0.038 -0.014 -0.015 -0.425 0.324 0.388 -0.176 1.000 0.241 0.010 -0.637 0.072 0.055 -0.004 -0.310 -0.134 -0.183 -0.144 -0.047 |\n",
"| p4415_p | 0.182 -0.016 -0.025 -0.150 0.011 0.320 0.003 -0.004 0.114 -0.015 0.230 0.055 0.241 1.000 -0.009 -0.246 0.260 -0.108 -0.003 0.163 -0.022 -0.049 -0.211 0.403 |\n",
"| omega_s | -0.018 0.565 0.017 -0.008 0.011 0.003 0.000 -0.017 0.004 0.005 -0.046 0.002 0.010 -0.009 1.000 0.005 -0.003 -0.014 -0.000 -0.016 -0.019 -0.109 -0.001 -0.001 |\n",
"| bplus_1 | -0.114 -0.011 0.015 -0.025 0.008 -0.135 0.007 0.006 0.214 -0.226 -0.497 0.034 -0.637 -0.246 0.005 1.000 -0.122 -0.012 0.002 0.110 -0.610 -0.032 0.197 -0.098 |\n",
"| p4040_p | 0.038 -0.010 -0.016 -0.188 0.006 0.384 0.010 -0.003 0.303 -0.016 0.006 0.099 0.072 0.260 -0.003 -0.122 1.000 0.300 -0.005 0.249 -0.014 -0.034 -0.014 0.348 |\n",
"| p4160_s | 0.178 0.012 -0.020 0.021 0.004 -0.036 -0.001 -0.002 -0.122 -0.014 0.143 -0.016 0.055 -0.108 -0.014 -0.012 0.300 1.000 -0.001 -0.077 0.051 0.028 0.336 -0.136 |\n",
"| p3770_p | 0.006 0.000 0.000 -0.003 -0.000 -0.002 -0.002 0.000 -0.008 -0.003 0.002 -0.002 -0.004 -0.003 -0.000 0.002 -0.005 -0.001 1.000 0.001 0.000 0.001 0.001 -0.003 |\n",
"| jpsi_p | 0.311 0.010 -0.033 -0.036 0.001 0.702 0.011 -0.012 0.425 -0.280 -0.263 0.119 -0.310 0.163 -0.016 0.110 0.249 -0.077 0.001 1.000 0.057 0.025 -0.018 0.459 |\n",
"| bplus_0 | 0.027 0.104 -0.011 0.021 -0.029 -0.032 0.001 -0.001 0.024 -0.055 0.100 -0.003 -0.134 -0.022 -0.019 -0.610 -0.014 0.051 0.000 0.057 1.000 0.268 0.044 -0.020 |\n",
"| rho_s | 0.012 0.448 -0.046 0.009 -0.062 -0.055 0.002 -0.028 0.035 -0.068 0.013 0.002 -0.183 -0.049 -0.109 -0.032 -0.034 0.028 0.001 0.025 0.268 1.000 0.050 -0.035 |\n",
"| p4415_s | 0.106 0.019 0.006 0.176 -0.006 -0.076 -0.002 0.001 -0.103 -0.074 -0.174 -0.035 -0.144 -0.211 -0.001 0.197 -0.014 0.336 0.001 -0.018 0.044 0.050 1.000 -0.115 |\n",
"| p4160_p | 0.112 -0.009 -0.012 -0.415 0.005 0.559 0.012 -0.000 0.457 -0.142 -0.046 0.135 -0.047 0.403 -0.001 -0.098 0.348 -0.136 -0.003 0.459 -0.020 -0.035 -0.115 1.000 |\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"Hesse errors: OrderedDict([(<zfit.Parameter 'Ctt' floating=True>, {'error': 0.23746033613438966}), (<zfit.Parameter 'omega_p' floating=True>, {'error': 0.22827854391421543}), (<zfit.Parameter 'phi_s' floating=True>, {'error': 1.3242924020583473}), (<zfit.Parameter 'p4040_s' floating=True>, {'error': 0.17019599757914744}), (<zfit.Parameter 'rho_p' floating=True>, {'error': 0.1332345481328261}), (<zfit.Parameter 'Dbar_p' floating=True>, {'error': 0.7962709805312216}), (<zfit.Parameter 'p3770_s' floating=True>, {'error': 0.02308522313498257}), (<zfit.Parameter 'phi_p' floating=True>, {'error': 0.2876950629670567}), (<zfit.Parameter 'psi2s_p' floating=True>, {'error': 0.03769146402222612}), (<zfit.Parameter 'DDstar_s' floating=True>, {'error': 0.40560761040792503}), (<zfit.Parameter 'bplus_2' floating=True>, {'error': 0.10018804096358758}), (<zfit.Parameter 'Dbar_s' floating=True>, {'error': 0.06870386529805039}), (<zfit.Parameter 'DDstar_p' floating=True>, {'error': 0.6338643163749156}), (<zfit.Parameter 'p4415_p' floating=True>, {'error': 0.2091282373871879}), (<zfit.Parameter 'omega_s' floating=True>, {'error': 1.0387903247081258}), (<zfit.Parameter 'bplus_1' floating=True>, {'error': 0.05705684069783579}), (<zfit.Parameter 'p4040_p' floating=True>, {'error': 0.20352338541949067}), (<zfit.Parameter 'p4160_s' floating=True>, {'error': 0.1706073621912625}), (<zfit.Parameter 'p3770_p' floating=True>, {'error': 0.018358184856082094}), (<zfit.Parameter 'jpsi_p' floating=True>, {'error': 0.03880632726508182}), (<zfit.Parameter 'bplus_0' floating=True>, {'error': 0.02288158646680416}), (<zfit.Parameter 'rho_s' floating=True>, {'error': 0.3649741145457665}), (<zfit.Parameter 'p4415_s' floating=True>, {'error': 0.19896385950009027}), (<zfit.Parameter 'p4160_p' floating=True>, {'error': 0.11669689762358448})])\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Toy 4/10\n",
"Time taken: 8 min, 31 s\n",
"Projected time left: 12 min, 42 s\n",
"Toy 4: Generating data...\n",
"Toy 4: Data generation finished\n",
"Toy 4: Loading data...\n",
"Toy 4: Loading data finished\n",
"Toy 4: Fitting pdf...\n",
"------------------------------------------------------------------\n",
"| FCN = 2.98E+05 | Ncalls=1319 (1319 total) |\n",
"| EDM = 1.7E-05 (Goal: 5E-06) | up = 0.5 |\n",
"------------------------------------------------------------------\n",
"| Valid Min. | Valid Param. | Above EDM | Reached call limit |\n",
"------------------------------------------------------------------\n",
"| True | True | False | False |\n",
"------------------------------------------------------------------\n",
"| Hesse failed | Has cov. | Accurate | Pos. def. | Forced |\n",
"------------------------------------------------------------------\n",
"| False | True | False | False | True |\n",
"------------------------------------------------------------------\n",
"Function minimum: 297957.0655247679\n",
"----------------------------------------------------------------------------------------------\n",
"| | Name | Value | Hesse Err | Minos Err- | Minos Err+ | Limit- | Limit+ | Fixed |\n",
"----------------------------------------------------------------------------------------------\n",
"| 0 | Ctt | -0.017 | 0.204 | | | -1.5 | 1.5 | |\n",
"| 1 | omega_p | 0.8 | 0.3 | | |-6.28319 | 6.28319 | |\n",
"| 2 | phi_s | 18 | 6 | | | 14.8182 | 23.5818 | |\n",
"| 3 | p4040_s | 1.39 | 0.16 | | |0.00501244| 2.01499 | |\n",
"| 4 | rho_p | -0.53 | 0.30 | | |-6.28319 | 6.28319 | |\n",
"| 5 | Dbar_p | 1.1 | 0.3 | | |-6.28319 | 6.28319 | |\n",
"| 6 | p3770_s | 1.62 | 0.22 | | |0.918861 | 4.08114 | |\n",
"| 7 | phi_p | -6 | 7 | | |-6.28319 | 6.28319 | |\n",
"| 8 | psi2s_p | -2.03 | 0.03 | | |-6.28319 | 6.28319 | |\n",
"| 9 | DDstar_s | -0.300 | 0.029 | | | -0.3 | 0.3 | |\n",
"| 10| bplus_2 | -0.20 | 0.08 | | | -2 | 2 | |\n",
"| 11| Dbar_s | -0.300 | 0.016 | | | -0.3 | 0.3 | |\n",
"| 12| DDstar_p | 2.6 | 0.3 | | |-6.28319 | 6.28319 | |\n",
"| 13| p4415_p | 4.6 | 0.5 | | |-6.28319 | 6.28319 | |\n",
"| 14| omega_s | 8.9 | 1.3 | | | 4.19232 | 9.40768 | |\n",
"| 15| bplus_1 | -0.93 | 0.03 | | | -2 | 2 | |\n",
"| 16| p4040_p | 3.02 | 0.11 | | |-6.28319 | 6.28319 | |\n",
"| 17| p4160_s | 0.717 | 0.011 | | | 0.71676 | 3.68324 | |\n",
"| 18| p3770_p | -1.85 | 0.15 | | |-6.28319 | 6.28319 | |\n",
"| 19| jpsi_p | -1.752 | 0.032 | | |-6.28319 | 6.28319 | |\n",
"| 20| bplus_0 | 0.480 | 0.016 | | | -2 | 2 | |\n",
"| 21| rho_s | 1.32 | 0.31 | | |0.0253049| 2.0747 | |\n",
"| 22| p4415_s | 0.40 | 0.18 | | |0.126447 | 2.35355 | |\n",
"| 23| p4160_p | -6.283 | 0.012 | | |-6.28319 | 6.28319 | |\n",
"----------------------------------------------------------------------------------------------\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"| | Ctt omega_p phi_s p4040_s rho_p Dbar_p p3770_s phi_p psi2s_p DDstar_s bplus_2 Dbar_s DDstar_p p4415_p omega_s bplus_1 p4040_p p4160_s p3770_p jpsi_p bplus_0 rho_s p4415_s p4160_p |\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"| Ctt | 1.000 -0.019 -0.070 0.243 0.092 -0.423 -0.091 -0.067 -0.344 0.026 0.652 0.010 0.348 0.098 -0.006 -0.129 0.105 -0.002 -0.209 0.574 -0.023 -0.018 0.145 -0.003 |\n",
"| omega_p | -0.019 1.000 0.306 0.001 0.031 0.007 -0.003 0.305 0.000 -0.001 -0.058 0.001 0.021 -0.017 0.720 -0.043 -0.015 -0.000 0.004 -0.040 0.058 0.489 0.002 0.000 |\n",
"| phi_s | -0.070 0.306 1.000 -0.010 -0.245 -0.016 -0.014 0.995 -0.000 0.000 -0.203 -0.001 0.001 -0.050 -0.036 -0.019 -0.053 -0.001 0.015 -0.112 -0.001 0.145 -0.011 0.001 |\n",
"| p4040_s | 0.243 0.001 -0.010 1.000 0.011 0.079 0.084 -0.010 -0.204 0.006 0.013 0.014 0.052 -0.165 -0.001 0.072 -0.247 -0.006 -0.191 0.070 0.024 0.016 0.038 0.010 |\n",
"| rho_p | 0.092 0.031 -0.245 0.011 1.000 0.037 0.039 -0.250 -0.038 0.000 0.288 0.007 0.059 0.072 0.071 0.004 0.057 0.001 -0.027 -0.020 0.045 0.172 0.009 -0.001 |\n",
"| Dbar_p | -0.423 0.007 -0.016 0.079 0.037 1.000 -0.073 -0.015 0.096 -0.000 -0.063 0.012 -0.099 -0.123 0.000 0.192 -0.316 -0.007 -0.017 -0.239 0.073 0.037 -0.006 0.004 |\n",
"| p3770_s | -0.091 -0.003 -0.014 0.084 0.039 -0.073 1.000 -0.012 0.245 -0.012 0.018 -0.020 -0.109 0.012 0.000 -0.033 0.144 -0.001 -0.153 -0.118 0.021 -0.012 0.057 -0.002 |\n",
"| phi_p | -0.067 0.305 0.995 -0.010 -0.250 -0.015 -0.012 1.000 -0.001 0.000 -0.194 -0.001 0.003 -0.048 -0.042 -0.020 -0.051 -0.001 0.014 -0.113 0.001 0.149 -0.011 0.001 |\n",
"| psi2s_p | -0.344 0.000 -0.000 -0.204 -0.038 0.096 0.245 -0.001 1.000 0.017 -0.201 0.027 0.248 0.036 -0.000 0.091 0.149 0.007 0.238 0.025 0.021 0.017 -0.098 -0.000 |\n",
"| DDstar_s | 0.026 -0.001 0.000 0.006 0.000 -0.000 -0.012 0.000 0.017 1.000 0.027 -0.001 0.039 0.015 -0.000 -0.030 0.041 0.001 -0.001 0.049 -0.009 -0.006 0.012 -0.001 |\n",
"| bplus_2 | 0.652 -0.058 -0.203 0.013 0.288 -0.063 0.018 -0.194 -0.201 0.027 1.000 0.016 0.059 0.282 -0.008 -0.237 0.067 0.002 -0.100 0.254 -0.031 -0.069 -0.100 -0.002 |\n",
"| Dbar_s | 0.010 0.001 -0.001 0.014 0.007 0.012 -0.020 -0.001 0.027 -0.001 0.016 1.000 -0.017 -0.007 -0.000 0.006 -0.003 -0.000 0.004 0.028 0.005 0.002 0.010 0.000 |\n",
"| DDstar_p | 0.348 0.021 0.001 0.052 0.059 -0.099 -0.109 0.003 0.248 0.039 0.059 -0.017 1.000 -0.042 -0.000 0.234 0.160 0.006 0.110 0.588 0.082 0.056 0.069 -0.000 |\n",
"| p4415_p | 0.098 -0.017 -0.050 -0.165 0.072 -0.123 0.012 -0.048 0.036 0.015 0.282 -0.007 -0.042 1.000 -0.002 -0.126 -0.008 0.006 0.062 0.045 -0.010 -0.028 -0.095 0.003 |\n",
"| omega_s | -0.006 0.720 -0.036 -0.001 0.071 0.000 0.000 -0.042 -0.000 -0.000 -0.008 -0.000 -0.000 -0.002 1.000 0.000 -0.002 -0.000 0.001 -0.009 0.001 0.104 -0.000 0.000 |\n",
"| bplus_1 | -0.129 -0.043 -0.019 0.072 0.004 0.192 -0.033 -0.020 0.091 -0.030 -0.237 0.006 0.234 -0.126 0.000 1.000 -0.007 -0.002 0.014 0.022 -0.822 -0.127 0.125 0.000 |\n",
"| p4040_p | 0.105 -0.015 -0.053 -0.247 0.057 -0.316 0.144 -0.051 0.149 0.041 0.067 -0.003 0.160 -0.008 -0.002 -0.007 1.000 -0.005 0.178 0.095 0.020 -0.011 0.138 -0.003 |\n",
"| p4160_s | -0.002 -0.000 -0.001 -0.006 0.001 -0.007 -0.001 -0.001 0.007 0.001 0.002 -0.000 0.006 0.006 -0.000 -0.002 -0.005 1.000 0.006 0.004 -0.000 -0.000 0.002 -0.001 |\n",
"| p3770_p | -0.209 0.004 0.015 -0.191 -0.027 -0.017 -0.153 0.014 0.238 -0.001 -0.100 0.004 0.110 0.062 0.001 0.014 0.178 0.006 1.000 0.009 0.001 0.004 -0.065 -0.002 |\n",
"| jpsi_p | 0.574 -0.040 -0.112 0.070 -0.020 -0.239 -0.118 -0.113 0.025 0.049 0.254 0.028 0.588 0.045 -0.009 0.022 0.095 0.004 0.009 1.000 -0.014 0.008 0.031 -0.001 |\n",
"| bplus_0 | -0.023 0.058 -0.001 0.024 0.045 0.073 0.021 0.001 0.021 -0.009 -0.031 0.005 0.082 -0.010 0.001 -0.822 0.020 -0.000 0.001 -0.014 1.000 0.196 0.041 -0.000 |\n",
"| rho_s | -0.018 0.489 0.145 0.016 0.172 0.037 -0.012 0.149 0.017 -0.006 -0.069 0.002 0.056 -0.028 0.104 -0.127 -0.011 -0.000 0.004 0.008 0.196 1.000 0.021 0.000 |\n",
"| p4415_s | 0.145 0.002 -0.011 0.038 0.009 -0.006 0.057 -0.011 -0.098 0.012 -0.100 0.010 0.069 -0.095 -0.000 0.125 0.138 0.002 -0.065 0.031 0.041 0.021 1.000 -0.006 |\n",
"| p4160_p | -0.003 0.000 0.001 0.010 -0.001 0.004 -0.002 0.001 -0.000 -0.001 -0.002 0.000 -0.000 0.003 0.000 0.000 -0.003 -0.001 -0.002 -0.001 -0.000 0.000 -0.006 1.000 |\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"Hesse errors: OrderedDict([(<zfit.Parameter 'Ctt' floating=True>, {'error': 0.20356743696537594}), (<zfit.Parameter 'omega_p' floating=True>, {'error': 0.3225986671314396}), (<zfit.Parameter 'phi_s' floating=True>, {'error': 5.601509675112658}), (<zfit.Parameter 'p4040_s' floating=True>, {'error': 0.1629207983331552}), (<zfit.Parameter 'rho_p' floating=True>, {'error': 0.3048294644485092}), (<zfit.Parameter 'Dbar_p' floating=True>, {'error': 0.3470474381158444}), (<zfit.Parameter 'p3770_s' floating=True>, {'error': 0.22325462827629838}), (<zfit.Parameter 'phi_p' floating=True>, {'error': 7.06425982044021}), (<zfit.Parameter 'psi2s_p' floating=True>, {'error': 0.03353790591960104}), (<zfit.Parameter 'DDstar_s' floating=True>, {'error': 0.028565280512765218}), (<zfit.Parameter 'bplus_2' floating=True>, {'error': 0.08327237332492732}), (<zfit.Parameter 'Dbar_s' floating=True>, {'error': 0.01630754606864146}), (<zfit.Parameter 'DDstar_p' floating=True>, {'error': 0.34158714572292137}), (<zfit.Parameter 'p4415_p' floating=True>, {'error': 0.47996263337494316}), (<zfit.Parameter 'omega_s' floating=True>, {'error': 1.2607557762653512}), (<zfit.Parameter 'bplus_1' floating=True>, {'error': 0.033473910114156125}), (<zfit.Parameter 'p4040_p' floating=True>, {'error': 0.11311693903164866}), (<zfit.Parameter 'p4160_s' floating=True>, {'error': 0.011131439035562851}), (<zfit.Parameter 'p3770_p' floating=True>, {'error': 0.1531198411151733}), (<zfit.Parameter 'jpsi_p' floating=True>, {'error': 0.031531811376411145}), (<zfit.Parameter 'bplus_0' floating=True>, {'error': 0.016322193507993354}), (<zfit.Parameter 'rho_s' floating=True>, {'error': 0.31173277656278864}), (<zfit.Parameter 'p4415_s' floating=True>, {'error': 0.18417991434425623}), (<zfit.Parameter 'p4160_p' floating=True>, {'error': 0.011534763945483029})])\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Toy 5/10\n",
"Time taken: 11 min, 1 \n",
"Projected time left: 11 min\n",
"Toy 5: Generating data...\n",
"Toy 5: Data generation finished\n",
"Toy 5: Loading data...\n",
"Toy 5: Loading data finished\n",
"Toy 5: Fitting pdf...\n",
"------------------------------------------------------------------\n",
"| FCN = 2.978E+05 | Ncalls=1557 (1557 total) |\n",
"| EDM = 2.49E-05 (Goal: 5E-06) | up = 0.5 |\n",
"------------------------------------------------------------------\n",
"| Valid Min. | Valid Param. | Above EDM | Reached call limit |\n",
"------------------------------------------------------------------\n",
"| True | True | False | False |\n",
"------------------------------------------------------------------\n",
"| Hesse failed | Has cov. | Accurate | Pos. def. | Forced |\n",
"------------------------------------------------------------------\n",
"| False | True | True | True | False |\n",
"------------------------------------------------------------------\n",
"Function minimum: 297819.47573215386\n",
"----------------------------------------------------------------------------------------------\n",
"| | Name | Value | Hesse Err | Minos Err- | Minos Err+ | Limit- | Limit+ | Fixed |\n",
"----------------------------------------------------------------------------------------------\n",
"| 0 | Ctt | 1.07 | 0.18 | | | -1.5 | 1.5 | |\n",
"| 1 | omega_p | -6.3 | 1.0 | | |-6.28319 | 6.28319 | |\n",
"| 2 | phi_s | 15.4 | 1.2 | | | 14.8182 | 23.5818 | |\n",
"| 3 | p4040_s | 1.39 | 0.17 | | |0.00501244| 2.01499 | |\n",
"| 4 | rho_p | 0.24 | 0.34 | | |-6.28319 | 6.28319 | |\n",
"| 5 | Dbar_p | -5.4 | 0.3 | | |-6.28319 | 6.28319 | |\n",
"| 6 | p3770_s | 2.33 | 0.24 | | |0.918861 | 4.08114 | |\n",
"| 7 | phi_p | 6.10 | 0.27 | | |-6.28319 | 6.28319 | |\n",
"| 8 | psi2s_p | -2.080 | 0.031 | | |-6.28319 | 6.28319 | |\n",
"| 9 | DDstar_s | -0.300 | 0.026 | | | -0.3 | 0.3 | |\n",
"| 10| bplus_2 | -0.16 | 0.09 | | | -2 | 2 | |\n",
"| 11| Dbar_s | -0.300 | 0.016 | | | -0.3 | 0.3 | |\n",
"| 12| DDstar_p | -5.4 | 0.5 | | |-6.28319 | 6.28319 | |\n",
"| 13| p4415_p | -2.31 | 0.16 | | |-6.28319 | 6.28319 | |\n",
"| 14| omega_s | 7.2 | 1.2 | | | 4.19232 | 9.40768 | |\n",
"| 15| bplus_1 | 0.90 | 0.07 | | | -2 | 2 | |\n",
"| 16| p4040_p | -2.87 | 0.13 | | |-6.28319 | 6.28319 | |\n",
"| 17| p4160_s | 2.10 | 0.17 | | | 0.71676 | 3.68324 | |\n",
"| 18| p3770_p | -2.10 | 0.12 | | |-6.28319 | 6.28319 | |\n",
"| 19| jpsi_p | 4.591 | 0.025 | | |-6.28319 | 6.28319 | |\n",
"| 20| bplus_0 | -0.45 | 0.03 | | | -2 | 2 | |\n",
"| 21| rho_s | 1.0 | 0.4 | | |0.0253049| 2.0747 | |\n",
"| 22| p4415_s | 1.38 | 0.19 | | |0.126447 | 2.35355 | |\n",
"| 23| p4160_p | 4.28 | 0.10 | | |-6.28319 | 6.28319 | |\n",
"----------------------------------------------------------------------------------------------\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"| | Ctt omega_p phi_s p4040_s rho_p Dbar_p p3770_s phi_p psi2s_p DDstar_s bplus_2 Dbar_s DDstar_p p4415_p omega_s bplus_1 p4040_p p4160_s p3770_p jpsi_p bplus_0 rho_s p4415_s p4160_p |\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"| Ctt | 1.000 0.012 -0.033 0.262 0.086 -0.343 -0.022 -0.008 -0.410 0.010 -0.591 0.000 -0.126 0.105 0.017 0.070 0.019 0.200 -0.262 0.446 0.062 0.020 0.105 -0.117 |\n",
"| omega_p | 0.012 1.000 0.053 -0.001 0.145 -0.004 0.003 0.024 -0.003 -0.002 -0.025 -0.000 -0.025 0.014 0.642 -0.062 0.008 0.004 -0.000 -0.006 0.073 -0.076 -0.007 0.012 |\n",
"| phi_s | -0.033 0.053 1.000 -0.001 0.139 -0.011 -0.017 0.738 0.016 -0.002 0.091 -0.002 -0.017 -0.015 0.070 -0.081 -0.012 -0.022 0.008 0.016 0.085 0.062 0.001 -0.001 |\n",
"| p4040_s | 0.262 -0.001 -0.001 1.000 -0.002 0.060 0.147 -0.001 -0.195 0.001 0.018 0.013 -0.125 -0.157 -0.002 -0.019 -0.176 -0.176 -0.171 0.054 0.004 -0.000 0.089 -0.368 |\n",
"| rho_p | 0.086 0.145 0.139 -0.002 1.000 -0.006 0.020 0.088 -0.019 -0.008 -0.166 0.000 -0.092 0.067 0.490 -0.267 0.039 0.028 -0.005 -0.003 0.311 0.071 -0.035 0.050 |\n",
"| Dbar_p | -0.343 -0.004 -0.011 0.060 -0.006 1.000 -0.062 -0.003 0.007 -0.008 -0.026 0.014 -0.210 -0.124 -0.013 -0.023 -0.191 0.002 -0.105 -0.222 -0.030 0.029 0.051 -0.201 |\n",
"| p3770_s | -0.022 0.003 -0.017 0.147 0.020 -0.062 1.000 0.001 0.260 -0.000 -0.006 -0.015 0.106 -0.023 -0.005 -0.017 0.042 0.084 -0.141 -0.089 -0.039 0.040 0.077 -0.072 |\n",
"| phi_p | -0.008 0.024 0.738 -0.001 0.088 -0.003 0.001 1.000 -0.002 -0.000 0.002 0.000 -0.005 0.001 0.007 -0.007 0.001 -0.001 -0.000 -0.020 0.005 0.074 -0.000 0.001 |\n",
"| psi2s_p | -0.410 -0.003 0.016 -0.195 -0.019 0.007 0.260 -0.002 1.000 0.019 0.125 0.026 -0.063 0.118 0.007 0.007 0.177 -0.061 0.243 -0.057 0.033 -0.045 -0.132 0.251 |\n",
"| DDstar_s | 0.010 -0.002 -0.002 0.001 -0.008 -0.008 -0.000 -0.000 0.019 1.000 0.019 -0.001 0.043 -0.005 -0.006 -0.021 -0.001 -0.003 0.015 0.031 -0.012 0.010 0.005 0.011 |\n",
"| bplus_2 | -0.591 -0.025 0.091 0.018 -0.166 -0.026 -0.006 0.002 0.125 0.019 1.000 -0.006 0.270 -0.265 -0.002 -0.370 -0.102 -0.082 0.077 -0.219 0.158 -0.164 0.189 -0.165 |\n",
"| Dbar_s | 0.000 -0.000 -0.002 0.013 0.000 0.014 -0.015 0.000 0.026 -0.001 -0.006 1.000 -0.003 -0.003 -0.002 -0.005 -0.005 0.002 0.005 0.028 -0.004 0.005 0.010 -0.004 |\n",
"| DDstar_p | -0.126 -0.025 -0.017 -0.125 -0.092 -0.210 0.106 -0.005 -0.063 0.043 0.270 -0.003 1.000 -0.338 -0.062 -0.177 -0.222 -0.098 0.073 -0.339 -0.127 0.084 0.009 -0.212 |\n",
"| p4415_p | 0.105 0.014 -0.015 -0.157 0.067 -0.124 -0.023 0.001 0.118 -0.005 -0.265 -0.003 -0.338 1.000 0.022 0.128 0.152 -0.085 0.086 0.170 0.011 0.004 -0.158 0.368 |\n",
"| omega_s | 0.017 0.642 0.070 -0.002 0.490 -0.013 -0.005 0.007 0.007 -0.006 -0.002 -0.002 -0.062 0.022 1.000 -0.185 0.010 -0.004 0.005 0.023 0.214 -0.208 -0.016 0.025 |\n",
"| bplus_1 | 0.070 -0.062 -0.081 -0.019 -0.267 -0.023 -0.017 -0.007 0.007 -0.021 -0.370 -0.005 -0.177 0.128 -0.185 1.000 0.091 0.008 -0.016 0.015 -0.924 0.227 -0.076 0.102 |\n",
"| p4040_p | 0.019 0.008 -0.012 -0.176 0.039 -0.191 0.042 0.001 0.177 -0.001 -0.102 -0.005 -0.222 0.152 0.010 0.091 1.000 0.341 0.187 0.070 -0.020 0.010 0.105 0.004 |\n",
"| p4160_s | 0.200 0.004 -0.022 -0.176 0.028 0.002 0.084 -0.001 -0.061 -0.003 -0.082 0.002 -0.098 -0.085 -0.004 0.008 0.341 1.000 -0.010 0.013 -0.035 0.040 0.331 -0.117 |\n",
"| p3770_p | -0.262 -0.000 0.008 -0.171 -0.005 -0.105 -0.141 -0.000 0.243 0.015 0.077 0.005 0.073 0.086 0.005 -0.016 0.187 -0.010 1.000 -0.048 0.024 -0.022 -0.090 0.191 |\n",
"| jpsi_p | 0.446 -0.006 0.016 0.054 -0.003 -0.222 -0.089 -0.020 -0.057 0.031 -0.219 0.028 -0.339 0.170 0.023 0.015 0.070 0.013 -0.048 1.000 0.184 -0.115 -0.054 0.117 |\n",
"| bplus_0 | 0.062 0.073 0.085 0.004 0.311 -0.030 -0.039 0.005 0.033 -0.012 0.158 -0.004 -0.127 0.011 0.214 -0.924 -0.020 -0.035 0.024 0.184 1.000 -0.275 -0.026 0.038 |\n",
"| rho_s | 0.020 -0.076 0.062 -0.000 0.071 0.029 0.040 0.074 -0.045 0.010 -0.164 0.005 0.084 0.004 -0.208 0.227 0.010 0.040 -0.022 -0.115 -0.275 1.000 0.009 -0.023 |\n",
"| p4415_s | 0.105 -0.007 0.001 0.089 -0.035 0.051 0.077 -0.000 -0.132 0.005 0.189 0.010 0.009 -0.158 -0.016 -0.076 0.105 0.331 -0.090 -0.054 -0.026 0.009 1.000 -0.067 |\n",
"| p4160_p | -0.117 0.012 -0.001 -0.368 0.050 -0.201 -0.072 0.001 0.251 0.011 -0.165 -0.004 -0.212 0.368 0.025 0.102 0.004 -0.117 0.191 0.117 0.038 -0.023 -0.067 1.000 |\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"Hesse errors: OrderedDict([(<zfit.Parameter 'Ctt' floating=True>, {'error': 0.17756068569158656}), (<zfit.Parameter 'omega_p' floating=True>, {'error': 0.974810495678307}), (<zfit.Parameter 'phi_s' floating=True>, {'error': 1.2048099140867805}), (<zfit.Parameter 'p4040_s' floating=True>, {'error': 0.1746148426443993}), (<zfit.Parameter 'rho_p' floating=True>, {'error': 0.3388255648988139}), (<zfit.Parameter 'Dbar_p' floating=True>, {'error': 0.33809702511271134}), (<zfit.Parameter 'p3770_s' floating=True>, {'error': 0.23753713150793}), (<zfit.Parameter 'phi_p' floating=True>, {'error': 0.2688131804807581}), (<zfit.Parameter 'psi2s_p' floating=True>, {'error': 0.031223719302520436}), (<zfit.Parameter 'DDstar_s' floating=True>, {'error': 0.025508419572656593}), (<zfit.Parameter 'bplus_2' floating=True>, {'error': 0.09027231656247292}), (<zfit.Parameter 'Dbar_s' floating=True>, {'error': 0.015661261261986303}), (<zfit.Parameter 'DDstar_p' floating=True>, {'error': 0.46545022942963765}), (<zfit.Parameter 'p4415_p' floating=True>, {'error': 0.1639990663169446}), (<zfit.Parameter 'omega_s' floating=True>, {'error': 1.248215198422895}), (<zfit.Parameter 'bplus_1' floating=True>, {'error': 0.06888293344375596}), (<zfit.Parameter 'p4040_p' floating=True>, {'error': 0.13379105407082625}), (<zfit.Parameter 'p4160_s' floating=True>, {'error': 0.17211479906979676}), (<zfit.Parameter 'p3770_p' floating=True>, {'error': 0.11563312352469346}), (<zfit.Parameter 'jpsi_p' floating=True>, {'error': 0.025354427208602814}), (<zfit.Parameter 'bplus_0' floating=True>, {'error': 0.03394709032162613}), (<zfit.Parameter 'rho_s' floating=True>, {'error': 0.354306489482555}), (<zfit.Parameter 'p4415_s' floating=True>, {'error': 0.19367485659297867}), (<zfit.Parameter 'p4160_p' floating=True>, {'error': 0.1046625681834934})])\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Toy 6/10\n",
"Time taken: 13 min, 46 s\n",
"Projected time left: 9 min, 8 s\n",
"Toy 6: Generating data...\n",
"Toy 6: Data generation finished\n",
"Toy 6: Loading data...\n",
"Toy 6: Loading data finished\n",
"Toy 6: Fitting pdf...\n"
]
}
],
"source": [
"# zfit.run.numeric_checks = False \n",
"\n",
"fitting_range = 'cut'\n",
"total_BR = 1.7e-10 + 4.9e-10 + 2.5e-9 + 6.02e-5 + 4.97e-6 + 1.38e-9 + 4.2e-10 + 2.6e-9 + 6.1e-10 + 4.37e-7\n",
"cut_BR = 1.0 - (6.02e-5 + 4.97e-6)/total_BR\n",
"\n",
"Ctt_list = []\n",
"Ctt_error_list = []\n",
"\n",
"nr_of_toys = 25\n",
"if fitting_range == 'cut':\n",
" nevents = int(pdg[\"number_of_decays\"]*cut_BR)\n",
"else:\n",
" nevents = int(pdg[\"number_of_decays\"])\n",
"# nevents = pdg[\"number_of_decays\"]\n",
"event_stack = 1000000\n",
"nevents *= 41\n",
"# zfit.settings.set_verbosity(10)\n",
"calls = int(nevents/event_stack + 1)\n",
"\n",
"total_samp = []\n",
"\n",
"start = time.time()\n",
"\n",
"sampler = total_f.create_sampler(n=event_stack)\n",
"\n",
"for toy in range(nr_of_toys):\n",
" \n",
" ### Generate data\n",
" \n",
"# clear_output(wait=True)\n",
" \n",
" print(\"Toy {}: Generating data...\".format(toy))\n",
" \n",
" dirName = 'data/zfit_toys/toy_{0}'.format(toy)\n",
" \n",
" if not os.path.exists(dirName):\n",
" os.mkdir(dirName)\n",
" print(\"Directory \" , dirName , \" Created \")\n",
" \n",
" reset_param_values()\n",
" \n",
" if fitting_range == 'cut':\n",
" \n",
" sampler.resample(n=nevents)\n",
" s = sampler.unstack_x()\n",
" sam = zfit.run(s)\n",
" calls = 0\n",
" c = 1\n",
" \n",
" else: \n",
" for call in range(calls):\n",
"\n",
" sampler.resample(n=event_stack)\n",
" s = sampler.unstack_x()\n",
" sam = zfit.run(s)\n",
"\n",
" c = call + 1\n",
"\n",
" with open(\"data/zfit_toys/toy_{0}/{1}.pkl\".format(toy, call), \"wb\") as f:\n",
" pkl.dump(sam, f, pkl.HIGHEST_PROTOCOL)\n",
" \n",
" print(\"Toy {}: Data generation finished\".format(toy))\n",
" \n",
" ### Load data\n",
" \n",
" print(\"Toy {}: Loading data...\".format(toy))\n",
" \n",
" if fitting_range == 'cut':\n",
" \n",
" total_samp = sam\n",
" \n",
" else:\n",
" \n",
" for call in range(calls):\n",
" with open(r\"data/zfit_toys/toy_0/{}.pkl\".format(call), \"rb\") as input_file:\n",
" sam = pkl.load(input_file)\n",
" total_samp = np.append(total_samp, sam)\n",
"\n",
" total_samp = total_samp.astype('float64')\n",
" \n",
" if fitting_range == 'full':\n",
"\n",
" data = zfit.data.Data.from_numpy(array=total_samp[:int(nevents)], obs=obs)\n",
" \n",
" print(\"Toy {}: Loading data finished\".format(toy))\n",
"\n",
" ### Fit data\n",
"\n",
" print(\"Toy {}: Fitting pdf...\".format(toy))\n",
"\n",
" for param in total_f.get_dependents():\n",
" param.randomize()\n",
"\n",
" nll = zfit.loss.UnbinnedNLL(model=total_f, data=data, fit_range = (x_min, x_max), constraints = constraints)\n",
"\n",
" minimizer = zfit.minimize.MinuitMinimizer(verbosity = 5)\n",
" # minimizer._use_tfgrad = False\n",
" result = minimizer.minimize(nll)\n",
"\n",
" print(\"Toy {}: Fitting finished\".format(toy))\n",
"\n",
" print(\"Function minimum:\", result.fmin)\n",
" print(\"Hesse errors:\", result.hesse())\n",
"\n",
" params = result.params\n",
" Ctt_list.append(params[Ctt]['value'])\n",
" Ctt_error_list.append(params[Ctt]['minuit_hesse']['error'])\n",
"\n",
" #plotting the result\n",
"\n",
" plotdirName = 'data/plots'.format(toy)\n",
"\n",
" if not os.path.exists(plotdirName):\n",
" os.mkdir(plotdirName)\n",
"# print(\"Directory \" , dirName , \" Created \")\n",
" \n",
" probs = total_f.pdf(test_q, norm_range=False)\n",
" calcs_test = zfit.run(probs)\n",
" plt.clf()\n",
" plt.plot(test_q, calcs_test, label = 'pdf')\n",
" plt.legend()\n",
" plt.ylim(0.0, 6e-6)\n",
" plt.savefig(plotdirName + '/toy_fit_full_range{}.png'.format(toy))\n",
"\n",
" print(\"Toy {0}/{1}\".format(toy+1, nr_of_toys))\n",
" print(\"Time taken: {}\".format(display_time(int(time.time() - start))))\n",
" print(\"Projected time left: {}\".format(display_time(int((time.time() - start)/(c+calls*(toy))*((nr_of_toys-toy)*calls-c)))))\n",
" \n",
" if fitting_range == 'cut':\n",
" \n",
" _1 = np.where((total_samp >= x_min) & (total_samp <= (jpsi_mass - 60.)))\n",
" \n",
" tot_sam_1 = total_samp[_1]\n",
" \n",
" _2 = np.where((total_samp >= (jpsi_mass + 70.)) & (total_samp <= (psi2s_mass - 50.)))\n",
" \n",
" tot_sam_2 = total_samp[_2]\n",
"\n",
" _3 = np.where((total_samp >= (psi2s_mass + 50.)) & (total_samp <= x_max))\n",
" \n",
" tot_sam_3 = total_samp[_3]\n",
"\n",
" tot_sam = np.append(tot_sam_1, tot_sam_2)\n",
" tot_sam = np.append(tot_sam, tot_sam_3)\n",
" \n",
" data = zfit.data.Data.from_numpy(array=tot_sam[:int(nevents)], obs=obs_fit)\n",
" \n",
" print(\"Toy {}: Loading data finished\".format(toy))\n",
" \n",
" ### Fit data\n",
"\n",
" print(\"Toy {}: Fitting pdf...\".format(toy))\n",
"\n",
" for param in total_f_fit.get_dependents():\n",
" param.randomize()\n",
"\n",
" nll = zfit.loss.UnbinnedNLL(model=total_f_fit, data=data, constraints = constraints)\n",
"\n",
" minimizer = zfit.minimize.MinuitMinimizer(verbosity = 5)\n",
" # minimizer._use_tfgrad = False\n",
" result = minimizer.minimize(nll)\n",
"\n",
" print(\"Function minimum:\", result.fmin)\n",
" print(\"Hesse errors:\", result.hesse())\n",
"\n",
" params = result.params\n",
" \n",
" if result.converged:\n",
" Ctt_list.append(params[Ctt]['value'])\n",
" Ctt_error_list.append(params[Ctt]['minuit_hesse']['error'])\n",
"\n",
" #plotting the result\n",
"\n",
" plotdirName = 'data/plots'.format(toy)\n",
"\n",
" if not os.path.exists(plotdirName):\n",
" os.mkdir(plotdirName)\n",
" # print(\"Directory \" , dirName , \" Created \")\n",
" \n",
" plt.clf()\n",
" plt.hist(tot_sam, bins = int((x_max-x_min)/7.), label = 'toy data')\n",
" plt.savefig(plotdirName + '/toy_histo_cut_region{}.png'.format(toy))\n",
"\n",
" \n",
" probs = total_f_fit.pdf(test_q, norm_range=False)\n",
" calcs_test = zfit.run(probs)\n",
" plt.clf()\n",
" plt.plot(test_q, calcs_test, label = 'pdf')\n",
" plt.axvline(x=jpsi_mass-60.,color='red', linewidth=0.7, linestyle = 'dotted')\n",
" plt.axvline(x=jpsi_mass+70.,color='red', linewidth=0.7, linestyle = 'dotted')\n",
" plt.axvline(x=psi2s_mass-50.,color='red', linewidth=0.7, linestyle = 'dotted')\n",
" plt.axvline(x=psi2s_mass+50.,color='red', linewidth=0.7, linestyle = 'dotted')\n",
" plt.legend()\n",
" plt.ylim(0.0, 1.5e-6)\n",
" plt.savefig(plotdirName + '/toy_fit_cut_region{}.png'.format(toy))\n",
" \n",
" print(\"Toy {0}/{1}\".format(toy+1, nr_of_toys))\n",
" print(\"Time taken: {}\".format(display_time(int(time.time() - start))))\n",
" print(\"Projected time left: {}\".format(display_time(int((time.time() - start)/(toy+1))*((nr_of_toys-toy-1)))))\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"with open(\"data/results/Ctt_list.pkl\", \"wb\") as f:\n",
" pkl.dump(Ctt_list, f, pkl.HIGHEST_PROTOCOL)\n",
"with open(\"data/results/Ctt_error_list.pkl\", \"wb\") as f:\n",
" pkl.dump(Ctt_error_list, f, pkl.HIGHEST_PROTOCOL)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print('{0}/{1} fits converged'.format(len(Ctt_list), nr_of_toys))\n",
"print('Mean Ctt value = {}'.format(np.mean(Ctt_list)))\n",
"print('Mean Ctt error = {}'.format(np.mean(Ctt_error_list)))\n",
"print('95 Sensitivy = {}'.format(((2*np.mean(Ctt_error_list))**2)*4.2/1000))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.hist(tot_sam, bins = int((x_max-x_min)/7.))\n",
"\n",
"plt.show()\n",
"\n",
"# _ = np.where((total_samp >= x_min) & (total_samp <= (jpsi_mass - 50.)))\n",
"\n",
"tot_sam.shape"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# sample from original values"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
