{ "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": 3, "metadata": {}, "outputs": [], "source": [ "def formfactor( q2, subscript): #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", " b0 = ztf.constant(pdg[\"b0\"])\n", " bplus = ztf.constant(pdg[\"bplus\"])\n", " bT = ztf.constant(pdg[\"bT\"])\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", "\n", " for i in range(N):\n", " _sum += b0[i]*(tf.pow(z,i))\n", "\n", " return tf.complex(prefactor * _sum, ztf.constant(0.0))\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", " b = bT\n", " else:\n", " b = bplus\n", "\n", " for i in range(N):\n", " _sum += b[i] * (tf.pow(z, i) - ((-1)**(i-N)) * (i/N) * tf.pow(z, N))\n", "\n", " return tf.complex(prefactor * _sum, ztf.constant(0.0))\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", "def bifur_gauss(q, mean, sigma_L, sigma_R, scale):\n", "\n", " _exp = tf.where(q < mean, ztf.exp(- tf.pow((q-mean),2) / (2 * sigma_L**2)), ztf.exp(- tf.pow((q-mean),2) / (2 * sigma_R**2)))\n", "\n", " #Scale so the total area under curve is 1 and the top of the cusp is continuous\n", "\n", " dgamma = scale*_exp/(ztf.sqrt(2*np.pi))*2*(sigma_L*sigma_R)/(sigma_L+sigma_R)\n", "\n", " com = ztf.complex(dgamma, ztf.constant(0.0))\n", "\n", " return com\n", "\n", "def axiv_nonres(q):\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 = ztf.sqrt(tf.abs(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. * kabs**2. * beta**2. *tf.abs(tf.complex(C10eff, ztf.constant(0.0))*formfactor(q2, \"+\"))**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(tf.complex(C10eff, ztf.constant(0.0)) * formfactor(q2, \"0\")), 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 *ztf.sqrt(q2)\n", "\n", "def vec(q, funcs):\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 = ztf.sqrt(tf.abs(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 = kabs**2 * (1. - 1./3. * beta**2)\n", "\n", " abs_bracket = tf.abs(c9eff(q, funcs) * formfactor(q2, \"+\") + tf.complex(2.0 * C7eff * (mb + ms)/(mB + mK), ztf.constant(0.0)) * formfactor(q2, \"T\"))**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 * ztf.sqrt(q2)\n", "\n", "def c9eff(q, funcs):\n", "\n", " C9eff_nr = tf.complex(ztf.constant(pdg['C9eff']), ztf.constant(0.0))\n", "\n", " c9 = C9eff_nr\n", "\n", " c9 = c9 + 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": [ "## C_q,qbar constraint" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# r = rho_scale * rho_width/rho_mass * np.cos(rho_phase)*(1-np.tan(rho_phase)*rho_width/rho_mass)\n", "# o = omega_scale*np.cos(omega_phase)*omega_width/omega_mass\n", "# p = phi_scale*np.cos(phi_phase)*phi_width/phi_mass\n", "\n", "# # phi_s = np.linspace(-500, 5000, 100000)\n", "\n", "# # p_ = phi_s*np.cos(phi_phase)*phi_width/phi_mass\n", "\n", "# # p_y = r+o+p_\n", "\n", "# # plt.plot(phi_s, p_y)\n", "\n", "# print(r + o + p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Build pdf" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "class total_pdf(zfit.pdf.ZPDF):\n", " _N_OBS = 1 # dimension, can be omitted\n", " _PARAMS = ['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", " 'DDstar_mass', 'DDstar_scale', 'DDstar_phase',\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", " 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 resonance(q, _mass = self.params['jpsi_mass'], scale = self.params['jpsi_scale'],\n", " phase = self.params['jpsi_phase'], width = self.params['jpsi_width'])\n", "\n", " def psi2s_res(q):\n", " return resonance(q, _mass = self.params['psi2s_mass'], scale = self.params['psi2s_scale'],\n", " phase = self.params['psi2s_phase'], width = self.params['psi2s_width'])\n", " \n", " def p3770_res(q):\n", " return resonance(q, _mass = self.params['p3770_mass'], scale = self.params['p3770_scale'],\n", " phase = self.params['p3770_phase'], width = self.params['p3770_width'])\n", " \n", " def p4040_res(q):\n", " return resonance(q, _mass = self.params['p4040_mass'], scale = self.params['p4040_scale'],\n", " phase = self.params['p4040_phase'], width = self.params['p4040_width'])\n", " \n", " def p4160_res(q):\n", " return resonance(q, _mass = self.params['p4160_mass'], scale = self.params['p4160_scale'],\n", " phase = self.params['p4160_phase'], width = self.params['p4160_width'])\n", " \n", " def p4415_res(q):\n", " return resonance(q, _mass = self.params['p4415_mass'], scale = self.params['p4415_scale'],\n", " phase = self.params['p4415_phase'], 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['DDstar_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)\n", "\n", " axiv_nr = axiv_nonres(x)\n", "\n", " tot = vec_f + axiv_nr\n", "\n", " return tot" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Load data" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "x_min = 2*pdg['muon_M']\n", "x_max = (pdg[\"Bplus_M\"]-pdg[\"Ks_M\"]-0.1)\n", "\n", "obs = zfit.Space('q', limits = (x_min, x_max))\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 parameters" ] }, { "cell_type": "code", "execution_count": 8, "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": [ "#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,\n", "# upper_limit = rho_mass + rho_width)\n", "rho_w = zfit.Parameter(\"rho_w\", ztf.constant(rho_width), floating = False)\n", "rho_p = zfit.Parameter(\"rho_p\", ztf.constant(rho_phase), floating = False) #, lower_limit=-2*np.pi, upper_limit=2*np.pi)\n", "rho_s = zfit.Parameter(\"rho_s\", ztf.constant(rho_scale), floating = False) #, lower_limit=rho_scale-np.sqrt(rho_scale), upper_limit=rho_scale+np.sqrt(rho_scale))\n", "\n", "#omega\n", "\n", "omega_mass, omega_width, omega_phase, omega_scale = pdg[\"omega\"]\n", "\n", "omega_m = zfit.Parameter(\"omega_m\", ztf.constant(omega_mass), floating = False)\n", "omega_w = zfit.Parameter(\"omega_w\", ztf.constant(omega_width), floating = False)\n", "omega_p = zfit.Parameter(\"omega_p\", ztf.constant(omega_phase), floating = False) #, lower_limit=-2*np.pi, upper_limit=2*np.pi)\n", "omega_s = zfit.Parameter(\"omega_s\", ztf.constant(omega_scale), floating = False) #, lower_limit=omega_scale-np.sqrt(omega_scale), upper_limit=omega_scale+np.sqrt(omega_scale))\n", "\n", "\n", "#phi\n", "\n", "phi_mass, phi_width, phi_phase, phi_scale = pdg[\"phi\"]\n", "\n", "phi_m = zfit.Parameter(\"phi_m\", ztf.constant(phi_mass), floating = False)\n", "phi_w = zfit.Parameter(\"phi_w\", ztf.constant(phi_width), floating = False)\n", "phi_p = zfit.Parameter(\"phi_p\", ztf.constant(phi_phase), floating = False) #, lower_limit=-2*np.pi, upper_limit=2*np.pi)\n", "phi_s = zfit.Parameter(\"phi_s\", ztf.constant(phi_scale), floating = False) #, lower_limit=phi_scale-np.sqrt(phi_scale), upper_limit=phi_scale+np.sqrt(phi_scale))\n", "\n", "#jpsi\n", "\n", "jpsi_mass, jpsi_width, jpsi_phase, jpsi_scale = pdg[\"jpsi\"]\n", "# jpsi_scale *= pdg[\"factor_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), 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), 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), floating = False) #, lower_limit=-2*np.pi, upper_limit=2*np.pi)\n", "p3770_s = zfit.Parameter(\"p3770_s\", ztf.constant(p3770_scale), floating = False) #, 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), floating = False) #, lower_limit=-2*np.pi, upper_limit=2*np.pi)\n", "p4040_s = zfit.Parameter(\"p4040_s\", ztf.constant(p4040_scale), floating = False) #, 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), floating = False) #, lower_limit=-2*np.pi, upper_limit=2*np.pi)\n", "p4160_s = zfit.Parameter(\"p4160_s\", ztf.constant(p4160_scale), floating = False) #, 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), floating = False) #, lower_limit=-2*np.pi, upper_limit=2*np.pi)\n", "p4415_s = zfit.Parameter(\"p4415_s\", ztf.constant(p4415_scale), floating = False) #, 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": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "<zfit.util.temporary.TemporarilySet at 0x1795ef3a2e8>" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "_0 = jpsi_scale*np.cos(jpsi_phase)*jpsi_width/jpsi_mass**3 + psi2s_scale*np.cos(psi2s_phase)*psi2s_width/psi2s_mass**3\n", "_1 = p3770_scale*np.cos(p3770_phase)*p3770_width/p3770_mass**3 + p4040_scale*np.cos(p4040_phase)*p4040_width/p4040_mass**3\n", "_2 = p4160_scale*np.cos(p4160_phase)*p4160_width/p4160_mass**3 + p4415_scale*np.cos(p4415_phase)*p4415_width/p4415_mass**3\n", "\n", "C_pert = np.random.uniform(0.03, 0.1)\n", "# c_pert = 0.1\n", "m_c = 1300\n", "\n", "cDDstar_phase = 10\n", "\n", "\n", "DDstar_eta = 0\n", "Dbar_phase = np.random.uniform(0.0, 2*np.pi)\n", "DDstar_phase = np.random.uniform(0.0, 2*np.pi)\n", "DDstar_mass = pdg['D0_M']\n", "\n", "if Dbar_phase < np.pi:\n", " Dbar_phase = 0.0\n", "else:\n", " Dbar_phase = np.pi\n", "\n", "R = (C_pert/(m_c**2) - ((_0 + _1 + _2)))\n", "\n", "Dbar_mass = (pdg['D0_M']+pdg['Dst_M'])/2\n", "\n", "Dbar_eta = R/np.cos(Dbar_phase)*(6*Dbar_mass**2)\n", "\n", "# print(np.cos(Dbar_phase))\n", "\n", "# cDDstar_phase = R_*10*DDstar_mass**2/DDstar_eta\n", "\n", "\n", "# print(Dbar_eta)\n", "\n", "\n", "Dbar_s = zfit.Parameter(\"Dbar_s\", ztf.constant(Dbar_eta), lower_limit=-1.464, upper_limit=1.464)\n", "Dbar_m = zfit.Parameter(\"Dbar_m\", ztf.constant(Dbar_mass), floating = False)\n", "Dbar_p = zfit.Parameter(\"Dbar_p\", ztf.constant(Dbar_phase), floating = False)\n", "DDstar_s = zfit.Parameter(\"DDstar_s\", ztf.constant(DDstar_eta), floating = False)\n", "DDstar_m = zfit.Parameter(\"DDstar_m\", ztf.constant(DDstar_mass), floating = False)\n", "DDstar_p = zfit.Parameter(\"DDstar_p\", ztf.constant(DDstar_phase), floating = False)\n", "\n", "Dbar_s.set_value(0.0)\n", "DDstar_s.set_value(0.0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Tau parameters" ] }, { "cell_type": "code", "execution_count": 10, "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=-0.5, upper_limit=0.5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Setup pdf" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "total_f = total_pdf(obs=obs, 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", " DDstar_mass = DDstar_m, DDstar_scale = DDstar_s, DDstar_phase = DDstar_p,\n", " Dbar_mass = Dbar_m, Dbar_scale = Dbar_s, Dbar_phase = Dbar_p,\n", " tau_mass = tau_m, C_tt = Ctt)\n", " \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": "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, 200000)\n", "\n", "probs = total_f.pdf(test_q)\n", "\n", "calcs_test = zfit.run(probs)\n", "res_y = zfit.run(jpsi_res(test_q))\n", "f0_y = zfit.run(formfactor(test_q,\"0\"))\n", "fplus_y = zfit.run(formfactor(test_q,\"+\"))\n", "fT_y = zfit.run(formfactor(test_q,\"T\"))" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ8AAAD8CAYAAACo9anUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO29eXhd1Xnv/3nPOTqaJ0uyjSTbsrGMJ4wxxkAgQAxhCAkmBG5MQsMvkEvTkjZtevsD2iT3hoQmtGnobRNKaCAho6FAgkMgQAKBEMB4Nh4wluVJkicN1jycYd0/9j5HR9KZJOuMej/Po8f7rL32WmtvS/t73ne9611ijEFRFEVRkokj1QNQFEVRph4qPoqiKErSUfFRFEVRko6Kj6IoipJ0VHwURVGUpKPioyiKoiSduMRHRK4Rkb0i0iAi94Q5nysiT9jnN4hIXci5e+3yvSJydaw2RWSu3cY+u013HH0sE5G3RGSXiLwrInkTeRiKoihKcogpPiLiBL4HXAssBm4RkcWjqt0BdBhj5gMPAg/Y1y4G1gJLgGuAh0TEGaPNB4AHjTH1QIfddrQ+XMBPgc8bY5YAlwOecT4HRVEUJYnEY/msAhqMMY3GmCFgHbBmVJ01wOP28VPAFSIidvk6Y8ygMeYA0GC3F7ZN+5rVdhvYbd4Qo4+rgB3GmO0Axpg2Y4wv/kegKIqiJBtXHHVqgCMhn5uACyLVMcZ4RaQTqLDL3x51bY19HK7NCuCUMcYbpn6kPhYARkReBKqwxO6fR9+EiNwJ3AlQWFh43sKFC+O49anBkNdPe+8QM0vVW6mM5N3mTqYX5zKjZOzvxqDHz/snupk9rYDS/JwUjM7iRPcgx7sGWFpdiohVtquli/KCHM4ozWdnSycAZ9eUpmyM2cLmzZtbjTFVk9FWPOIjYcpG5+SJVCdSeTiLK1r9aH24gEuA84E+4PcistkY8/sRFY15BHgEYOXKlWbTpk1hmpuaXP3g63Qd72bdly5l/vTiVA9HSRP8fsO8f3iev7mynr+5csGY840ne1j9r6/x7bXLWbO8JkwLyeG7r+zj2y+9z1vfuBa3y3q1nHvfS3zsnGruuXYhi7/6IgAbv/kRRMK9RpR4EZFDk9VWPG63JmBWyOdaoCVSHXsOphRoj3JtpPJWoMxuY3Rf0fp4zRjTaozpA54HVsRxX4qN1+9P9RCUNMRn5310OcK/sF0O6/Xh9aU2P2QgPWWorjgdDrx+gydkbJ4Uj1MZSTzisxGot6PQ3FgBBOtH1VkP3GYf3wS8YqyMpeuBtXak2lygHngnUpv2Na/abWC3+WyMPl4ElolIgS1KlwG7438EiqKEIyAqTkf414TTab3tff4Ui4/9ryNEfVwOweczI8bm8emXrHQiptvNnl/5AtZL3gk8ZozZJSL3AZuMMeuBR4GfiEgDljWy1r52l4g8iSUGXuCuQDBAuDbtLu8G1onIN4CtdttE6aNDRL6DJWgGeN4Y85vTeiqKogQt4siWj9j1Uis+ftv0CR2l0yF4/QZviOCo+KQX8cz5YIx5HsudFVr21ZDjAeDmCNfeD9wfT5t2eSNWNNzo8mh9/BQr3FpRlEkiYDW4nOHFx+kIWD6pfamHc7u5nILP7x8hjIl2u3k8HpqamhgYGEhoP8kgLy+P2tpacnISF0gSl/goijL1CLy4093yCUYkhajPsOUzPLahBFs+TU1NFBcXU1dXl9GBDcYY2traaGpqYu7cuQnrR9PrKGNCFxUFhi2fiHM+jjSZ8zGG0e96l0Pw+c2IYBpvgsVnYGCAioqKjBYesES8oqIi4Racio+iKGGJZfk408XyMWPXYQSi3ZLpdgMyXngCJOM+VHwURQmLLxjtFmvOJ9VuNzMi0g1CLB+fRrulKyo+iqKEJRjtFiHgILDOJ9Xi4zeMcbsF53xGuN3UwQzwhz/8gY9+9KMADA4OcuWVV7J8+XKeeOKJpI5DAw6UsKkjFGV4zif8b0igOD3cbuEsn5HRbokOOMhEtm7disfjYdu2bUnvWy0fRQMOlLDEmvMRkeBLPpUYxgYcOB2C1zfS7ZbogIN04ODBgyxcuJDbbruNZcuWcdNNN9HX18dvf/tbFi5cyCWXXMIzzzwDwIkTJ7j11lvZtm0by5cvZ//+/Ukdq1o+iqKEJVa0m3VO0sPyGR3t5hQGPf4Rbrdkptf52q93sbula1LbXFxdwv/+2JKY9fbu3cujjz7KxRdfzO233853vvMdvv/97/PKK68wf/58PvnJTwIwffp0fvCDH/Dtb3+b5557blLHGg9q+SiKEpZYlk/gnC/lud3MGLdbMNotNOBgiuQwnDVrFhdffDEAt956K5s2bWLu3LnU19cjItx6660pHqGFWj6KooQl4E6LNOcTOJcOls/oIQai3UbkdvMmT3zisVASxegw6c7OzrQMAVfLR1GUsATcVNEsH6f9kk8lVrTbaMtH7KzWIdFuKR5nsjh8+DBvvfUWAL/4xS+48sorOXDgQHBO5xe/+EUqhxdExUdRlLDEinazzjlS/lI3mLGLTMUKhJiKWa0XLVrE448/zrJly2hvb+dv//ZveeSRR7juuuu45JJLmDNnTqqHCKjbTVGUCATnfJyRv6O6HII/1eITbp2P07Z8kpzhIB1wOBw8/PDDI8quueYa3nvvvTF1L7/8ci6//PIkjWwkavkoIaSfX1hJHb4YWypAusz5mDFut4Ao+vy6pUK6ouKjKEpYvDHS68Dw1gWpxBAlw8EUW+dTV1fHzp07Uz2MuFDxURQlLLH284H0sHz8xuCMlNttRIaDxI/TmOxw7SXjPlR8lBCy4w9HmRziXueTcvEZG+3mcjrw+EZvJpdYyycvL4+2traMF6DAfj55eXkJ7UcDDhRFCUt8GQ7SINrNmDHrfNxOB0Ne/whXW6LdbrW1tTQ1NXHy5MmE9pMMAjuZJhIVHyUEDThQhonH8nE60iCrtZ8xWyq4XQ6GfCNDrRPtdsvJyUnozp/ZhrrdFEUJS3wZDlJv+fgjWD4e38g5n6kQcJBJqPgoihKWeOd8Ur3OJ9ycT47Tgc9vGPQMh4trqHV6oeKjKEpY4stwICnfJ8cYw+hpKbfLKuj3+ADIy3FOmUWmmYKKj6IoYfEGc7tFfk24nY6Uu7Mst9toy8f63D/kxSGQ63Ko5ZNmqPgoGmGthCWwF44zyjqfHGc6rPMZG3CQG2L5uBwOXE7RbbTTDBUfRVHCEs+cT44d0pxK/GZsYtEcOx9d35APp0PIcarlk26o+CiKEhZfHFsq5KSBOytcYtHgnM+QD5dTrOi3KbKlQqag4qMoSli8cQQcBEKaU4kh3JzPsOXjcggupyR1MzklNnGJj4hcIyJ7RaRBRO4Jcz5XRJ6wz28QkbqQc/fa5XtF5OpYbYrIXLuNfXab7mh9iEidiPSLyDb7Z2QucSU2urZUCYPPb62fibYLZo4z9SHMkRaZAvR5fDgdDnKcjuAclpIexBQfEXEC3wOuBRYDt4jI4lHV7gA6jDHzgQeBB+xrFwNrgSXANcBDIuKM0eYDwIPGmHqgw247Yh82+40xy+2fz4/rCSiKEhaP3x+0ICLhSoO5FL8xY91uzoDbzWtbPo6kJBZV4icey2cV0GCMaTTGDAHrgDWj6qwBHrePnwKuEOvr0hpgnTFm0BhzAGiw2wvbpn3NarsN7DZviNGHoigJwOM1wZd4JNxpEXAQxfIJzvlIykPClZHEIz41wJGQz012Wdg6xhgv0AlURLk2UnkFcMpuY3RfkfoAmCsiW0XkNRH5YBz3pISiXwiVMHh8fnJc0V8R6RBqHW6RacBiG/DYcz6O1FtoykjiEZ9w1sXo37ZIdSarPFofR4HZxphzgS8BPxeRktEVReROEdkkIpuyIeusoiQaj88fXKwZiXQIYQ63yDTU8nE6xI7K029Z6UQ84tMEzAr5XAu0RKojIi6gFGiPcm2k8lagzG5jdF9h+7Bdem0AxpjNwH5gweibMMY8YoxZaYxZWVVVFcdtK8rUZsgXe84nx452S+UeNuFzu1mfrWg3Bzma2y3tiEd8NgL1dhSaGyuAYP2oOuuB2+zjm4BXjPXbuB5Ya0eqzQXqgXcitWlf86rdBnabz0brQ0Sq7AAGRGSe3Udj/I9AUZRweHxxzPnYFkYqrYpwWa1zQ9yFLqe1yDRZGQ6++8o+HnvjQFL6ymRi7udjjPGKyBeAFwEn8JgxZpeI3AdsMsasBx4FfiIiDVgWz1r72l0i8iSwG/ACdxljfADh2rS7vBtYJyLfALbabROpD+BS4D4R8QI+4PPGmPaJPxJFUQA83ngsH+ut7/H5g0KUbEyYgIPQcQfX+STB8vH7Dd9+6X0Abr9E9/aJRlybyRljngeeH1X21ZDjAeDmCNfeD9wfT5t2eSNWNNzo8rB9GGOeBp6OeRNKTDR2UAnFCjiI/ksRSDqaSpdW2P18QoTQ6QhkOEj8GI93DwSPO/s8lBbkJLzPTEUzHCiKEpa45nzsl3wqt1Ww1vmMCjgYYfk47AwHiXe7tZwaFp8TIUKkjEXFR1GUsHjiEB930O2W2oCD0ZZPYe6wU8ftSl6Gg64BT/C4tWco4f1lMio+SpAUBiwpaUg8AQcBcUrlAk4TJtQ61+UIupFzbfFJxmLYrv5h8WnrHUx4f5mMio+iKGGJd51PoG6qCJfhQEQodFvWT26Ow85Bl/hvV10D3uBxm1o+UVHxURQlLB6fiWudD8BQEuZTIhEutxtAgdsJWPM/rmS53UIsn86QY2UscUW7KYoy9YgnvY7bNRxqnSrCWT4wLD65LueIxbCJTAnZ1e8h1+XA6RC6B1R8oqHioyhKWDw+f8w5n3QItTZhQq0BCmy3m9vlCC46HfL5yXU5EzaWrgEPJfk5OAS6+r2xL5jCqPgoIWjEgTKMtcg0vjmfVIdaR7d8HOTlWMcDngSLT7+XkjwXDpERkW/KWHTOR1GUsAzFMecz7HZL4ZyPP/yGd/mBOR+Xg7yc4SzXiaRrwENxXg4l+TkqPjFQy0cJoqHWSijxrPNJh1DrcBkOYHj77+K8HPJty6d/KLHi0zvopSjXRY5TaOvVaLdoqOWjKEpYMiXUOlxuNxiejyrJdw273byJFZ++IR8FbifFeTkjIt+UsajlowRRw0cJZTyWTyq3qPaH2UwOoDjPDjhwOpJm+fQN+SjMdVGY6xyx5kcZi1o+iqKMwRgT1zofd3CdT3rldgO4cUUNRbkuzp1dPiLgIJH0DXkpcDspsS2fVO5zlO6o5aMoyhgCAQSxtklI1kR+NCK53T5YX8W2r34Yl9MRXHOT6HH2DlqWT0l+Dl6/od/jC4Z8KyNRy0cJol/SlACBOZxYcz6BsOXBFFs+4QIOAFy2ZRaIfOtPoPj4bLHJz7EsH9C1PtFQ8VEUZQzD4hP9FZGbBpZPpAwHoeS5Am63xI0zIGyFuU5K8i1rR8OtI6PiowQxGnKg2ATmcGKKj+2WS7XlEytjTjIsn74hy8opcLtCLB8Vn0io+CiKMoaAmAQm6iMhIuS6HAym4ZxPKHlJiHbrGxy2fErzLfHR5KKRUfFRguicjxIg4J4KBBREIy/HmWK3W+Q5nwCB+0ikhdYbavnY4qNut8io+CiKMobASzqePGh5OY6Uu91iWT5upwOHJNjysdu2Qq3tOR8NOIiIio+iMz3KGAa98Vs+ua5UWz7EnPMREfJynAmd8+kdHGv5qNstMio+ShB1uykBAosx47V8Er14Mxp+f2zLByxRCAQFJIKA5VOYa+0fVOB2asBBFFR8FEUZQ8DyyY2xyBSsOZ/BBOdMi4bPmGAS0WiU5LkSmvImYPkEtu8uzc9RyycKKj5KEA21VgIELJlY0W5gCVQqLR+fPz7xKc5z0Z1A8Qm49AL7CJXk6bYK0VDxUUjcpsJKppJRlo/f4IzD7Vacl5PQra17BwPiY1k+JfkuDTiIgoqPEkTnfJQAg+OyfJypt3xipAECy/LpSaDl0zfkRWQ4SEPdbtFR8VEUZQyB6LV4LJ/cHEfC98mJRvyWT2Ldbr2DPgrdrmCGbXW7RUfFR9GZHmUMwXU+8SwydTmDllIqiDfgoCg3sW63wHYKAUrU8olKXOIjIteIyF4RaRCRe8KczxWRJ+zzG0SkLuTcvXb5XhG5OlabIjLXbmOf3aY7Vh/2+dki0iMi/2u8D0FRlJGMf5Fpaiwfv99gDHEHHPQO+fD5E/N1q9feSC5ASX4OPYNe/AnqL9OJKT4i4gS+B1wLLAZuEZHFo6rdAXQYY+YDDwIP2NcuBtYCS4BrgIdExBmjzQeAB40x9UCH3XbEPkJ4EHgh3htXFCUyAx4fOU6J66Weyjkfnz1RGa/bDUjYvE//kDe4YypYod3GkFBXXyYTj+WzCmgwxjQaY4aAdcCaUXXWAI/bx08BV4jl+FwDrDPGDBpjDgANdnth27SvWW23gd3mDTH6QERuABqBXfHfujIaDThQAgx6/XFZPZBayydgxTjiWudjZR3oHkyMK6xn0Eth7vAzK9X8blGJR3xqgCMhn5vssrB1jDFeoBOoiHJtpPIK4JTdxui+wvYhIoXA3cDXot2EiNwpIptEZNPJkydj3LKiTG0GPL64UuuAFRHn8Rm8vuRbPwHxccUjPvYeO4mah+kZ9FI0yu2WyP4ynXh+u8L9r47+jhypzmSVR+vja1huup4w54crGvOIMWalMWZlVVVVtKpTFl1kqgQYj+UTmGTvS0F+t6DbLQ7xqSjKBaCtZyghY+nq9wYFB0IsHxWfsMSzuXgTMCvkcy3QEqFOk4i4gFKgPca14cpbgTIRcdnWTWj9SH1cANwkIv8MlAF+ERkwxnw3jntTFCUMg15/XJFuQHCSvXfQG3RtJYvAZH5c4lPoBqCtdzAhY+ke8ATnlQDKCqxn0dGn4hOOeH67NgL1dhSaGyuAYP2oOuuB2+zjm4BXjDHGLl9rR6rNBeqBdyK1aV/zqt0GdpvPRuvDGPNBY0ydMaYO+Dfgn1R4FOX0GPD44rZ8hsUn+ZaPdzzik0DLxxhD94CX4hDxrbT7a+1JjNhlOjEtH2OMV0S+ALwIOIHHjDG7ROQ+YJMxZj3wKPATEWnAskbW2tfuEpEngd2AF7jLGOMDCNem3eXdwDoR+Qaw1W6bSH0ok4cGHCgB+oa8FLrjFB+7XiCxZjIJWD7xZLUuyXOR4xRaEyA+Ax4/Xr8ZYfmUF7gRgTYVn7DE43bDGPM88Pyosq+GHA8AN0e49n7g/njatMsbsaLhRpdH7COkzv+Jdl5RlPjoHfSNeJFGI5DLrDeB2xVEwjuOgAMRoaIwNyFiEFi8Gmr5OB3CtAI3rb2JmWPKdDTDgRJEDR8lQP+QL7g1QCyKUuh2G0+oNUBFkZu2BIhBYKuGklGCXVmUGLHLBlR8FEUZQ++oVDHRKLDXtiRyo7ZI+E38lg9Y8z6JmIMJWD6jAy4qitwJcfNlAyo+ShCjkz6KTf+QLygqsQhYPj0pmPMZT8ABwIziXI53DUz6OAKWz2hXZYVaPhFR8VEUZQy9Q9643W7BdT4pcLuNJ+AAoKY8n+Ndg5OekSGwlid0nQ9AZZE7YeuKMh0VH0VRRuDzGwY8fvLjdbu5U2/5xOt2qynLB+Doqcm1ftrteaTAWqIAlUW5dA96g1tUKMOo+ChB1OmmwPB20PFaPk6HkJ/jTMmcz3gDDmrKLfFpPtU/7r5++vYhzr3vJbYe7hhzrq1nEIdAWcFI8QmIka71GYuKj6IoI+izLZh4LR+wFpr2pMLtNs6Ag9qyAgCaO8YvPo/96QAdfR7+64+NY8619g5RXuAeM/c0szQPICHzTJmOio8SROMNFLD2pQFGZGiORWGuMyWLTL3jtHxmlubhEGjq6BtXP+29QzSe7AXgj++3jkmi2tYzSEWRe8x1ATdf8yS7+bIBFR9FUUYQcJ/l58TndgPLRZfKDAfx7OcD4HY5qCnPp7G1d1z97D9p5S2+6bxauge97GrpGnG+vXeIisLcMdedYYtPywTcfNmOio8Sgpo+CvRNwPIpznOlZNO08WypEGDB9GL2HY+aBH8MAUvppvNqAXirsW3E+ZPd4S2folwXJXkuFZ8wqPgoijKCgAVTEGfAAVgZnE/1Jz+keLwBBwD1M4ppbO3BM479h5raLfFYPquMM6sKeTtEfPx+Q8upgWAww2iqy/JVfMKg4qMoyggClk+8GQ7A2rsmFZumDdkCkuOM/1W2YEYRHp/hUFv8rremjn6qinPJy3Fy0ZkVbDzQHpz3OdkzyJDPT215Qdhra8ryadE5nzGo+ChBNOBAgcgLJqNRVuDmVAr2rfH6rF9a97jEpxiA9451x31N06k+am3L5sJ5FfQO+Xi3udM6Z7vkassiWz4TCe3OdlR8FEUZQXeEJJnRKM3PYdDrT/piyoDrLMc1jjmfGcW4XQ62HzkV9zVNHf1By+bCeRXA8LxPIApudkV4y2fWtHw6+z26nfYoVHyUIGr4KABdAx4cEv8iUxjeMjrZL9iA283liP9V5nY5WFZTyuZDYxeLhsPnN7Sc6g9aPpVFudRPL+Kt/Zb47D7aRX6Ok7qKwrDXz6ssAqDx5PiCHLIdFR9FUUbQ1e+hOC9nXJP4gS2jk+1680zA7QawYk45O5u74srxdqJ7AI/PBMUHYPXC6by1v43WnkF2NHWy6IziiMlN51VZorT/5PjCu7MdFR8liM75KGBlaC7Jj9/qASjLt8KMk235eCfgdgNYMbuMIZ+fd5s6Y9ZtsrMh1ITM6dx0Xi1ev+G7rzSw9XAHl9RXRbx+1rQCXA5Ry2cUKj6Kooygq99DcW78wQYw7HY71ZfccGvPBKLdAC6aV4nTIfxh78mYdQOpeEItn/oZxVy5aAY/evMgADcsr454fY7TwZyKguDckGKh4qPoPj7KCLonYvkUpGrOx/rdHa/4lBbkcN7scl7deyJm3UCkWvWoaLZ/uWkZd146j4c+fR7zqoqitjGvqiiYJSEaO5s7+deX9o5J35ONqPgoQVSEFLACDkbvyBmLUlt8OlJm+YzP7QZw+cIqdrV0cawz+hqc5lP9TCt0j1l0W17o5h8+sohrls6M2ddZM4ppbO2NGQ34j7/ayX+80sCb+9ui1ssGVHwUJM68WMrUoKvfM641PgDFuS5yXY6kbxnt8U7M7QZwzRJLNNZvb45ar7mjn+qyvPEPLoSlNSX4/Ia9MdYWBRa+jicMPFNR8VGCqN2jgB1wME7LR0SoKs7lRJK3DvBMILdbgHlVRSyfVcbTm5ujWv3Np/pHBBtMhKU1pQDBhanh8PkNPfYaqwPjTHyaiaj4KIoSxOvz0zM4/jkfgOnFuZxM8qZpHp8ft9MxYev9E+fVsvd4N9sjRL15fX4OtfVGXMMTLzVl+ZQV5LCrJbL4tPYMBreIaJrAfkOZhoqPEkSnfJT2vvDbQceDZfkkWXy8flwTmO8JcMPyaopzXfzX62M3iAM43N6Hx2eYPz16QEEsRISza0qjWj6BwIZphe5x7zeUiaj4KBpooARp77XEZ1qYvWliMb04LyWWz0TmewIU5+Vw60VzeGHn0bCuroYTVoTa6YoPWBmx9xztpifCvkeBzNer6qZxtGuAIW92R7yp+Cg616MEabMDBsLtTROLquJcTvV54soaMFl4/Oa0xAfgsxfXkety8i8vvjfm3O6jXYhMjvhcMLcCn99ETOsTEJ/z507DmOzfelvFRwliVIamPG29p+d2A5Ia8ebx+nGfhtsNLIvtLy8/k+ffPcabDa0jzm062MHCmSUUjzMAIxwr5pThcsiIvYBCaTk1QFGui4Uzrazb2Z4JOy7xEZFrRGSviDSIyD1hzueKyBP2+Q0iUhdy7l67fK+IXB2rTRGZa7exz27THa0PEVklItvsn+0i8vGJPoypinrdlADtttusomj8brcZJdY1sdbNTCZDPj85rtP/Dv0/L53HnIoC/v6pHcEsDb2DXjYf6mBVXflptw/W5nzLakvZEFF8rJDu6imy9XbM/zURcQLfA64FFgO3iMjiUdXuADqMMfOBB4EH7GsXA2uBJcA1wEMi4ozR5gPAg8aYeqDDbjtiH8BOYKUxZrndx/dFZPyhOor63xTaeodwCJSNc50PENxyIJmT5QMeH3mu+De9i0RejpN/X3suJ7oH+MufbaFvyMv67S30e3x85OwzJmGkFhfOq2BHU2fYTBAtnf1Ul+VzRqm1pqg5yyPe4vnKsApoMMY0GmOGgHXAmlF11gCP28dPAVeIFfu4BlhnjBk0xhwAGuz2wrZpX7PabgO7zRui9WGM6TPGBGbw8tBXqKJMmLbeIcoL3OPKaB0gkPvsSHvyxKff4ydvHDuuRuOcWWU88IllvN3Yxoe/8zpff243y2eVsWrutElpH6xs2F6/4fX3x+aUazk1wBml+eTlOKksctPSqeJTAxwJ+dxkl4WtYwtBJ1AR5dpI5RXAqRAxCe0rUh+IyAUisgt4F/h8yPVBROROEdkkIptOnoydTHAqoqqttPcMTSjYACy3UmVRLkfak/fStCyfyZu6vnFFLT++/QLmVhbywfpKHvr0iknNAHLu7HIqCt38bs/xEeV9Q17ae4eCAm7tfprdAQfxuKfCPfnR76lIdSKVh/ttiVY/6jiMMRuAJSKyCHhcRF4wxoz4nzPGPAI8ArBy5Up9z4aggQZKgBPdA1ROYL4nwKxp+RxJstutvGBiYhmJS+oruaS+clLbDOB0CKsXTufFXcdGhIkH5ncCmRSqS/NpyPItGOL5ytAEzAr5XAu0RKpjz7eUAu1Rro1U3gqUhczZhPYVqY8gxpg9QC+wNI77Umw04EAJcKzTcv1MlFnlBRxOotttwOMjP2dy3G7J4sOLZ9A14OVPIZF1wT2DQiyfllP9Wb0GLx7x2QjU21FobqwAgvWj6qwHbrOPbwJeMdZTWw+stSPV5gL1wDuR2rSvedVuA7vNZ6P1YbfhAhCROcBZwMG4n4ASJIt/z5U48PkNx7sHgxPeE2H2tAKOdiZvgeSAx09eTmatGLnsrCrKCnJ4estwQtPRG9ZVl+XRN+RL+hYVySTm/5o9f/IF4EVgD/CkMWaXiNwnItfb1R4FKkSkAfgScI997ZZSTqEAACAASURBVC7gSWA38FvgLmOML1Kbdlt3A1+y26qw247YB3AJsF1EtgG/BP7SGDMyWF+JioqOAnCyexCf3zDzNMSnfkYRPr+hsTU5LqN+j4/8SQo4SBa5LidrzqnmxV3HguLScKKHAreTmSXWsw+IUDav9YkrJNkY8zzw/Kiyr4YcDwA3R7j2fuD+eNq0yxuxouFGl4ftwxjzE+AnMW9CiYnO/UxtjnYGNk2buPicZS+Q3Husm4UzSyZlXNEY8PjInYRQ62Rz03mzePytQ/z3piN87oPzeO9YF2fNLA5GGQ6v9RlgSXVpKoeaMDLLXlUUJWEEFofOLJn4nM+8yiJcDuG9GPvWTBYDHh95GTbnA3B2bSkXzJ3GI6830j3gYUdTJ2fXDItMQHyaszjBqIqPEkTdb1ObFlt8TmfOx+1yMH96UcxN0yYDn9/g8ZmMCzgI8MUr6znRPciND71J35CPD501PXiussiN2+UI/p9kIyo+iqIA1uLQQreTsoLTy2O26IwSdjZ3JjxSq3fIWs5XkGFzPgE+cGYld1wyl30nelhWW8qlC6qC50SEmrJ8nfNRsptsDudU4udgWy9zKgpPe1Hlijnl/HJrM00d/cyaVjBJoxtLt73rZ3Fe5r7GvnzdIm5ZNZva8nyco7JKVJflZXV+N7V8lCAqQVObg629zK08vR07Ac63E3FuPNgeo+bp0T1gRYpNRsbpVCEizJ9eFHbeqro0X8VHyW5UdBSPz09TRz91ladvqSyYXkxxnisJ4pP5lk80qsvyOdE9mNT9kZKJio8SRN1vU5fmjn68fsOcitO3fBwO4cJ5Fbz+fmtCf6eGLZ/sFJ+5lYUYA4fasjPiTcVHUZTgFtKT4XYD+PCiGTSf6mf30a5JaS8cw5ZP5rrdohHYPTWwlXe2oeKjBFG7Z+qy55glEgtmFE9Ke6sXTUcEXtp1PHblCdJli09JfnZaPmdWFSEC+46r+ChZinrblN0tXdSW51M6gU3kwlFZlMtF8yp4eksTfn9ifsE67R1HS7LU8sl3O6kpy97s1io+iqbVUdhztItFZ0xuOpxPXTCbpo5+Xt+XmP2zWnuGKM5zZWSGg3ipn17EvuPJyRaRbFR8lGFUg6Yk/UM+DrT2Trr4XLV4JlXFuTz82v6EBB6c7B6kqnjiew9lAvUzimls7cWXIOsxlaj4KOp2m+LsaunEb2BJ9eSKj9vl4K7Lz+TtxnZe3zf5ieZP9gxSdRob32UC9dOLGPL6OdjWm+qhTDoqPkoQdb9NTTYe7ABg5ZzySW/7lgtmM3taAV99die9g2N2tz8tWrsHqcxyyyeQ0Xpnc2eKRzL5qPgoKjlTnI0H2zmzqpCKBFgRuS4n/3zTMg639/GVX+2cNPeb329oPtVP9WkkQc0E6mcU4XY5VHyU7Ebdb1MPv9+w6WA759dNS1gfF86r4ItX1PPM1ma+9cJ7kyJAR7sGGPT6qZukdUnpSo7TwaKZxexsTtx6qVSRnQHyiqLExbvNnXQNeLlgXuLEB+CLV9TT2jPI919v5FjXAN+4YelpLQ49GFgUOwkZGdKdpTWlrN/egjHmtJO+phNq+Shq8Uxhfv/eCRwCly2YHrvyaSAifH3NUv7+6rNYv72Fqx98nfXbWyYcxbW7xbIE6idpUWw6s7SmlO4BL4fbsyvNjoqPEkRFaOrxynvHWTG7nGmF7oT3JSLc9aH5PP0XH6A4L4e//sVWrnrwNX74pwO09w6Nq60thzuYNS0/60OtgeAOpzuasmveR8VHQUMOpiZNHX3sbO5i9aLEWj2jWTG7nBe++EG++6lzKcx18bVf7+aCf/odt/9oI09sPExbz2DU6wc8Pt7Y18pF8yqSNOLUsmBGMbkuB9uOnEr1UCYVnfNRgqgETS2e2dIMwPXnVCe9b4dD+Oiyaj66rJr3jnXx1KYmXth5jFfeO4FD3uX8umlcvWQmH1o4nbqKghFzHeu3tdA96OVjKRh3KnC7HJxTW8bmQx2pHsqkouKjqLttCmKM4ZktTVw4bxq15YnbbTQeFs4s4csfXcw/XreIXS1dvLTrGC/uOs59z+3mvud2M2taPpctqGJZbRld/R4efPl9Vswu45L5lSkddzJZMaecR99oZMDjy5p0Qio+ShDdz2fq8Ob+Ng629XHXh+aneihBRISlNaUsrSnlS1edxaG2Xl5//ySvvX+SZ7Y089O3DwOwfFYZD336vKyK/IrFeXPKefg1w7vNnQkNi08mKj6KMgX5/uuNVBblprXrak5FIX92USF/dlEdQ14/xzoHcDqF6tK8KSU8YIkPwOZDHVkjPhpwoOhczxRjV0snr79/ks9eXJcxLhy3y8HsigJqyvKnnPAATCt0M6+ykE0Hs2feR8VHCaIilP0YY/jWC+9RVpDDrRfOSfVwlHGwYk45Ww53ZI17XMVHyZpfZiU2f9h7kj/ua+WvV9dP2sZxSnI4b0457b1DHGzLjsWmcYmPiFwjIntFpEFE7glzPldEnrDPbxCRupBz99rle0Xk6lhtishcu419dpvuaH2IyIdFZLOIvGv/u3qiD2OqoxqU3XQPePjyr3Yyr6pQrZ4MJDDvs+lge4pHMjnEFB8RcQLfA64FFgO3iMjiUdXuADqMMfOBB4EH7GsXA2uBJcA1wEMi4ozR5gPAg8aYeqDDbjtiH0Ar8DFjzNnAbcBPxvcIFNWcqcH9v9nD0c5+/uWmc3C71OmRacyvKqI0Pydr1vvE8xu4CmgwxjQaY4aAdcCaUXXWAI/bx08BV4g1K7gGWGeMGTTGHAAa7PbCtmlfs9puA7vNG6L1YYzZaoxpsct3AXkikv05NyaRYYtHZShbeXpzE+s2HuHOS88MfoNWMguHQzhvTjkbp4rlA9QAR0I+N9llYesYY7xAJ1AR5dpI5RXAKbuN0X1F6iOUTwBbjTFj8nOIyJ0isklENp08mZg95RUlHdl25BT3/vJdLppXwd9dtSDVw1FOg5V15ew/2TvuXHjpSDziEy6ucfRX5Eh1Jqs85jhEZAmWK+7Pw9TDGPOIMWalMWZlVVVVuCqKknW8d6yLz/7wHWaU5PK9T68gx6nutkxm5RxrjU82uN7i+U1sAmaFfK4FWiLVEREXUAq0R7k2UnkrUGa3MbqvSH0gIrXAL4HPGGP2x3FPShg04CC7aDjRza0/2IDb5eCnd1yQlMzVSmJZVluK2+lg06HMd73FIz4bgXo7Cs2NFUCwflSd9ViT/QA3Aa8YK353PbDWjlSbC9QD70Rq077mVbsN7DafjdaHiJQBvwHuNcb8aTw3r1hoqHX2sfFgO5/4z7cA4Wefu5A5U2DTtalAXo6TpTUlWbHYNKb42PMrXwBeBPYATxpjdonIfSJyvV3tUaBCRBqALwH32NfuAp4EdgO/Be4yxvgitWm3dTfwJbutCrvtiH3Y7cwHviIi2+yf5OaIz3AC2qMSlB08t6OFT/9gAxWFbp75iw8wf3pRqoekTCLn103j3aZOBjy+VA/ltIgrt5sx5nng+VFlXw05HgBujnDt/cD98bRplzdiRcONLg/bhzHmG8A3Yt6EEhG/Wj5ZwZDXzzdf2MMP/3SQ8+aU84PPrKRcXW1Zx3lzyvn+640Zn2RUE4sqTHAnYyWNONzWx1+t28r2I6e4/eK53HPtQl3Lk6UMLzbN7CSjKj5K0PJRAyjz8PsNP37rIA/8di8uh/Cfn17BtWefkephKQmkoiiXeVWFdqaDM1M9nAmj4qOo6GQoB1p7ufvpHbxzoJ3LFlTxzRvPprosP9XDUpLA+XOm8eLuY/j9BocjM7N8q/gow5aPhhxkBH1DXr73agP/9foBcnMc/MtNy7jpvNopudXAVOW8unKe2HSE/Sd7qJ9RnOrhTAgVH0UDDjIEYwy/3nGUbz6/h6OdA9x4bg13X7uQGSV5qR6akmQCcz2bDnWo+CiZSyDgQDUofdl6uINvvvAe7xxoZ0l1Cf9xy7mszODJZuX0qKsooKLQzcaD7dyyanaqhzMhVHwUJY15/3g3335xLy/tPk5lkZv7P76UtefPxpmhfn5lchCxkoxmcpodFR9FSUOaOvr4t9/t45ktTRS6Xfzdhxdw+yVzKczVP1nF4vy6aby0+zgnugeYXpx5rlf9TVaCqNct9bSc6ufh1/az7p0jIPC5D87jLy47UxeLKmM4r85a77P5YEdGhter+ChKGnC4rY//fK2BpzY3YQzcvLKWv1pdr6HTSkSWVpeS63KwUcVHyXQ0wWjy2X+yh4de3c+vtjXjFGHt+bP5/OVnUqOio8TA7XJwzqwyNmdohmsVH0VJAe8d6+J7r+7nNztacLsc3HZRHX9+2TwNm1bGxfl15Tz8WiN9Q14K3Jn1Os+s0SpKBmOM4a39bXz/9UZee/8kBW4n//PSeXzuknlUFevO78r4WTlnGj7/frYdOcUHzqxM9XDGhYqPoiQYj8/P8+8e5b/+2MjO5i4qi9z83YcXcOuFczSQQDktVswuR8RKMqrioygKAD2DXp7YeITH3jhA86l+5lUV8q0bz+aGc2vIy3GmenhKFlBakMOC6cVsysD1Pio+ShCNN5gcjncN8KM3D/Kztw/RNeBl1dxpfO36JaxeOD1jk0Aq6cvKunKe3daCz28yavGxio+iTALGGLYeOcWP/nSQ5989it8Yrl16Bp/74FzOnV2e6uEpWczKunJ+tuEwe491s7i6JNXDiRsVHyWIZrUeP4NeH7/ZcZQfvXmQHU2dFOe5uO0Dddx2UR2zKwpSPTxlCnDuLOvLzfamUyo+ipLtnOga4KcbDvPzDYdo7RnizKpCvn7DUm48t0ZT4ChJZU5FAaX5OexoOpVRSUb1ryQLMMbw/vEezpp5eqnVdc4nNlsPd/CjNw/ymx1H8RnD6rOm8/9dXMcl8yt1Px0lJYgIy2pL2XakM9VDGRcqPlnAuo1HuPeZd/n55y7gA/PHF26pWQ1i0z/k47kdLfz07UNsb+qkONfFZy6q4zMXzaGusjDVw1MUltWW8vBrjfQP+ch3Z0YkpYpPFrCjyfrG09jaO27x8fhUfCLRcKKHn284zFObj9A14LVca2uW8PEVtRSpa01JI86pLcPnN+w+2sl5czJjnyf9C8oCvD4/ADnO8bt9vH5/8FiNIBjy+nlx1zF+tuEQbze2k+MUrll6BrdeMJtVc6epa01JS86ZVQbA9iMqPkoS8dpbkbocjnFf6/Gq4gAcae/j5+8c5r83HaG1Z4hZ0/K5+5qF3LyylsoiTX2jpDczSvKYUZLLjqZTqR5K3Kj4ZAEe2/JxTcDy8YRaPpM2oszA6/Pzynsn+NmGw7y+7yQOEa5YOJ1PXziHD86v1AWhSkaxrLYs6ILPBFR8sgBP0O02AcvH549dKcs43NbHk5uO8NTmJo51DTCzJI8vXlHPJ8+fxRmlupWBkpmcU1vKy7uP09nvoTQ/J9XDiYmKTxbgtYMGJvI93RsScJDNkW8DHh+/3XmMJzYe4a3GNhwCly2o4mtrlnDFwum4JiDcipJOLKu15n3eberkkvr0TzIa11+ciFwjIntFpEFE7glzPldEnrDPbxCRupBz99rle0Xk6lhtishcu419dpvuaH2ISIWIvCoiPSLy3Yk+iExmyLZeAnM/E7k2W9nZ3MlXfrWTVff/jr95YhvNp/r5X1ct4M17ruCHn13F1UtmqvAoWcGy2lLAynSQCcS0fETECXwP+DDQBGwUkfXGmN0h1e4AOowx80VkLfAA8EkRWQysBZYA1cDvRGSBfU2kNh8AHjTGrBORh+22/zNSH8AA8BVgqf0z5QhYL/4JWC7Z6Hbr7PPwq23NPLHxCLuPdpHrcnDt0pn8j/NnceHcCp3LUbKSsgI3cyoKMiboIB632yqgwRjTCCAi64A1QKj4rAH+j338FPBdsWJS1wDrjDGDwAERabDbI1ybIrIHWA18yq7zuN3uf0bqwxjTC7whIvPHcd9ZRSBc2juBNTv9Q77gcSY73fx+w5v723hi0xFe3HWMIa+fpTUlfH3NEq5fXpMRPnBFOV2W1Zax+WBmbKsdj/jUAEdCPjcBF0SqY4zxikgnUGGXvz3q2hr7OFybFcApY4w3TP1IfbTGcQ9ZTWChqG8CbrdQ8clEDrb28szWZp7Z0kRTRz+l+Tnccv4s/sf5s1hSXZrq4SlKUllaXcKvt7fQ3jvEtDTfqDAe8Qnnoxj9lotUJ1J5OCd7tPrxjiMiInIncCfA7NmZk3wvHgKuM98E3G69oeKTIaZPZ5+H595t4ZktzWw+1IFD4OL5lfz91Wdx9ZKZulGbMmUJfOHa1dLJB+urUjya6MQjPk3ArJDPtUBLhDpNIuICSoH2GNeGK28FykTEZVs/ofUj9REXxphHgEcAVq5cmSGv2fgIuNsmEnDQN+SNXSkN8Pj8vLb3JM9sbeJ3u08w5POzYEYR9167kDXLa5hZmpfqISpKyllib6mwq6UrK8RnI1AvInOBZqwAgk+NqrMeuA14C7gJeMUYY0RkPfBzEfkOVsBBPfAOlhUzpk37mlftNtbZbT4brY+J3XZ2EbR8JhA80Ddizie9Hqcxhp3NXTy9pYlfb2+hrXeIikI3t144hxtX1LCkukTT3ShKCOWFbmrK8tnV0pXqocQkpvjY8ytfAF4EnMBjxphdInIfsMkYsx54FPiJHVDQjiUm2PWexApO8AJ3GWN8AOHatLu8G1gnIt8AttptE6kPu62DQAngFpEbgKtGReNlNf0eS0AmkiO0dzD9LJ9jnQP8aps1j/P+8R7cTgcfXjyDG1fUcOmCqgktplWUqcLi6hJ2taR/poO4FpkaY54Hnh9V9tWQ4wHg5gjX3g/cH0+bdnkjwxFxoeXR+qiLegNZTs+AJSA+//gtn1N9nskezoToG/Ly4q5jPLOlmTcaWjEGVs4p558+fjbXnX0GpQUaraYo8bCkuoTf7TlO76A3rTc2TN+RKXHTbVsvE5nzaesdDB4n24np9xvebmzj6S3NvLDzKH1DPmZNy+evV9fz8XNrdK8cRZkAS6tLMQb2HO1iZV36ZrhW8clwBjzDczb+CYhPa88Qpfk5dPYnzwJqONHDM1ua+NXWZlo6ByjOdXH9OdV84rxaVs4p13kcRTkNltQMBx2o+CgJo7Vn2HKZiOVzonuQyiI3nf2ehIYbtPcO8evtLTyzpYntTZ04HcKl9ZXc+5FFfHjxDA2PVpRJYmZJHtMK3Wk/76Pik+Ecae8PHo93kakxhgMne1g1dxr7T/ZO9tAY9Pp4Zc8Jnt7SzB/2nsDrNyw+o4QvX7eI65dXM71Yw6MVZbIREZZUl6R9xJuKT4ZzpL0veDxe8WnvHaJrwMu8qiLYc2JS5nyMMWw53MHTW5r5zY6jdPZ7mF6cyx2XzOXjK2pYOLPk9DtRFCUqS6pLefSNRoa8ftyu9IwOVfHJcHY0n6Io18WQ1z9u8QlsPBVYmDaRDAkBmjr6eHpzM89sbeJQWx/5OU6uXjKDG1fUcvH8SpyazFNRksaS6hI8PsP7x7tZWpOeaaZUfDIYYwxvNrRx7uwythzqGPecz1uNbeQ4hfPtScnxLlId8vr53Z7j/OKdw7zRYKXYu2heBX+1up5rls6kKI3DPBUlmwl8odzd0qXio0w+Gw6009jay59fNo/tR06Ny/Lx+Q2/3t7CxfMrKcqzfg3iFa8Drb384p3DPL25ibbeIapL8/jr1fXcvLKW2vKCCd2LoiiTR11FIYVupx10MCtm/VSg4pOhDHh83Pfr3cwsyeNj51TzrRfeG5f4PL2liaOdA/zvjy3BZbvEol1vjOGtxjYe/eMBXtl7AqcIVy6awSdXzeLS+ip1qylKGuFwCIurS9iZxkEHKj4ZiNfn5++e3M6eY13815+tpMDtwulwxG25HG7r45+e38O5s8u4esmMqDuhGmP4w96T/OvLe9nZ3EVFoZu/Wl3PrRfO1mg1RUljllSX8uSmI/j8Ji2/HKr4ZBgDHh9f+PlWfrfnOP/4kUVcuXgGAC6HxLXI9EBrL595bAPGwL99cjkiQo7DioYZvRnd+8e7+fIvd/LOwXZmTyvgWzeezQ3n1uiaHEXJAJZUl9A35ONAay/zpxelejhjUPHJIA609nLXz7aw51gX961ZwmcuqgueczokquVjjOHZbS189dmduJwOfnz7KuZUWOlrHA5BZDg3nDGGR984wAO/fY/ivBy+ccNSPnn+LE3oqSgZRCDQYFdLp4qPMjF8fsPP3znMt57fQ47LwQ8+s5IrFs0YUcfpkIiJRd871sU//3Yvr7x3ghWzy3jwk8uDwhPA5RA8foMxhvue280P/3SQqxbP4Js3nk1FUW7C7k1RlMQwf3oRbpeDXS1drFleE/uCJKPik+ZsOdzB1369m+1HTvGBMyv49s3nUF2WP6aeyyFjtlQ41NbLgy+/z7PbWyjKdfEPH1nIHZfMC+v/tcTL8PibB/nhnw7y2Yvr+Mp1i3Gkoa9YUZTY5DgdLJpZzM7m9Eyzo+KTpmw93MG//W4fr71/ksqiXP7v2uVcf051xKSbjhDL52T3IP/39++z7p0juJzCn196Jp+/bB5lBZH3dM9xODjS3sdP3jrEh86qUuFRlCxgSU0pz21vwRiTdgl7VXzSiCGvnxd2HuUnbx1i06EOphW6uefahfzZhXNi7svhcghen2H99hb+4Zl3GfD4WLtqFn+1up4ZJbGj0pxO4YWdx3A6hPvWLFXhUZQsYGl1KT/fcJimjn5mTUuvNXgqPmnA3mPdPLutmSc3NdHaM8icigK+fN0iblk1O+7NoJwO4Y2GVl7ec5yVc8p54BPLrJxtcRKQmg+dVZV2v6SKokyMpfb2CjubO9Pu71rFJ0XsP9nDc9uP8tyOFvad6MEhcPlZ0/mzi+ZwWX3VuC0Pp0PoG/Ixe1oBP779AvLd4wuH7rB3NP3EitpxXacoSvqyYEYxToews6WTa88+I9XDGYGKTxI53NbHr3e08NyOo+w52oUInD9nGl9fs4Rrlp5BVfHEo8oaTvQA8NmL68YtPABXLZ7BH/e1snrR9AmPQVGU9CIvx0n99CJ2NqdfpgMVnwTTcqqf3+ywLJztdhbpc2eX8ZWPLua6s89gZunkZAm4bEEVL+w8xg0TDKn8l5vPoX/IR65LF5AqSjaxtKaUP+w9kXZBByo+CaC1Z5Df7DjK+u0tbD7UAcDZNaXce+1Crlt2RkKSb963Zil3fWg+5YWRI9qiUZqfQ2l+ziSPSlGUVLO0uoSnNjdxonswruCjZKHiM0kEItWe2tzEm/vb8PkNC2cW8/dXn8V1Z59BXWVh7EZOg6ri3NNy2ymKkp0EMh3sbO5U8ckmege9PPbGAX789iFOdg8ya1o+n79sHtefU8NZM4tTPTxFUaY4i84oQQR2NneNyYySSlR8ToMX3j3KV57dRWvPIJefVcVnL57LB+dX6hoZRVHShsJcF/MqC9nZkl6ZDlR8JoAxhgd+u5eHX9vPstpSHvnMeayYXZ7qYSmKooRlaU0pGw+0p3oYI9A0xRPg+6838vBr+/nUBbN5+i8+oMKjKEpac3ZNKS2dAxzrHEj1UIKo+IyTxpM9/OtLe/nI2TO5/4alus2AoihpzwVzKwDYcKAtxSMZRt+c4+QHbxzAIcLXrl+aVjHziqIokVhcXUJxrou3G9PH9RaX+IjINSKyV0QaROSeMOdzReQJ+/wGEakLOXevXb5XRK6O1aaIzLXb2Ge36Z5oH5ONMYaXdx/nqiUzNaxZUZSMwekQVtaVZ5blIyJO4HvAtcBi4BYRWTyq2h1AhzFmPvAg8IB97WJgLbAEuAZ4SEScMdp8AHjQGFMPdNhtj7uP8T6IeDjRPcjJ7kFWztE5HkVRMouL51fSeLI3mIor1cRj+awCGowxjcaYIWAdsGZUnTXA4/bxU8AVYvmk1gDrjDGDxpgDQIPdXtg27WtW221gt3nDBPuYdFp7BgGYUaJWj6IomcX1y6txOoR///2+VA8FiC/UugY4EvK5CbggUh1jjFdEOoEKu/ztUdcGko+Fa7MCOGWM8YapP5E+gojIncCd9sceEWkDWiPedRSufWAiV6U1lUzwWWQh+iws9DkMk1XP4j+A//jUhC6tBOZM1jjiEZ9ws+omzjqRysNZXNHqT6SPkQXGPAI8EvgsIpuMMSvDXDvl0GcxjD4LC30Ow+izsLCfQ91ktReP260JmBXyuRZoiVRHRFxAKdAe5dpI5a1Amd3G6L7G24eiKIqSpsQjPhuBejsKzY01ub9+VJ31wG328U3AK8YYY5evtSPV5gL1wDuR2rSvedVuA7vNZyfYh6IoipKmxHS72fMrXwBeBJzAY8aYXSJyH7DJGLMeeBT4iYg0YFkja+1rd4nIk8BuwAvcZYzxAYRr0+7ybmCdiHwD2Gq3zUT6iMEjsatMGfRZDKPPwkKfwzD6LCwm9TmIZTwoiqIoSvLQDAeKoihK0lHxURRFUZLOlBSfWOmCsgEReUxETojIzpCyaSLysp266GURKbfLRUT+3X4eO0RkRcg1t9n194nIbeH6SmdEZJaIvCoie0Rkl4h80S6fUs9CRPJE5B0R2W4/h6/Z5WmbzirR2NlWtorIc/bnKfksROSgiLwrIttEZJNdlvi/D2PMlPrBCnDYD8wD3MB2YHGqx5WA+7wUWAHsDCn7Z+Ae+/ge4AH7+CPAC1hrpi4ENtjl04BG+99y+7g81fc2zudwBrDCPi4G3sdK6TSlnoV9P0X2cQ6wwb6/J4G1dvnDwF/Yx38JPGwfrwWesI8X238zucBc+2/Jmer7m+Az+RLwc+A5+/OUfBbAQaByVFnC/z6mouUTT7qgjMcY8zpWVGAooSmKRqcu+rGxeBtrrdUZwNXAy8aYdmNMB/AyVv68jMEYc9QYs8U+7gb2YGXAmFLPwr6fQFKvHPvHkMbprBKJiNQC1wE/sD+ndWqvFJDw+PvPqAAAAjNJREFUv4+pKD7h0gWNSceTpcwwxhwF66UMTLfLIz2TrHpWtrvkXKxv/VPuWdhupm3ACayXw37iTGcFhKazyujnYPNvwP8P+O3Pcaf2IvuehQFeEpHNYqUhgyT8fUzFbbTjSsczxTit1EWZgIgUAU8Df2OM6ZLIezFl7bMw1vq35SJSBvwSWBSumv1v1j4HEfkocMIYs1lELg8Uh6ma9c/C5mJjTIuITAdeFpH3otSdtGcxFS2fqZyO57htImP/e8IuH28apIxCRHKwhOdnxphn7OIp+SwAjDGngD9g+eynYjqri4HrReQgltt9NZYlNBWfBcaYFvvfE1hfSlaRhL+PqSg+8aQLylZCUxSNTl30GTuS5UKg0za1XwSuEpFyO9rlKrssY7B9848Ce4wx3wk5NaWehYhU2RYPIpIPXIk1/zXl0lkZY+41xtQaK0nmWqx7+zRT8FmISKGIFAeOsX6vd5KMv49UR1qk4gcrYuN9LJ/3P6Z6PAm6x18ARwEP1reSO7D81L8H9tn/TrPrCtbmfvuBd4GVIe3cjjWR2gB8NtX3NYHncAmW+b8D2Gb/fGSqPQtgGVa6qh32y+Wrdvk8rBdmA/DfQK5dnmd/brDPzwtp6x/t57MXuDbV93aaz+VyhqPdptyzsO95u/2zK/A+TMbfh6bXURRFUZLOVHS7KYqiKClGxUdRFEVJOio+iqIoStJR8VEURVGSjoqPoiiKknRUfBRFUZSko+KjKIqiJJ3/B8fWcVIpEC9fAAAAAElFTkSuQmCC\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, 6e-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=200000, mc_sampler=None)\n", "# inte = total_f.integrate(limits = (2000, x_max), norm_range=False)\n", "# inte_fl = zfit.run(inte)\n", "# print(inte_fl)\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][3]*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", "## Toys" ] }, { "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": [ "# Analysis" ] }, { "cell_type": "code", "execution_count": 39, "metadata": { "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "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 = 4.145E+04 | Ncalls=189 (189 total) |\n", "| EDM = 0.000434 (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", "Toy 0: Fitting finished\n", "Function minimum: 41446.17176214419\n", "---------------------------------------------------------------------------------------------\n", "| | Name | Value | Hesse Err | Minos Err- | Minos Err+ | Limit- | Limit+ | Fixed |\n", "---------------------------------------------------------------------------------------------\n", "| 0 | psi2s_p | -4.4 | 0.4 | | |-6.28319 | 6.28319 | |\n", "| 1 | Dbar_s | -0.7 | 1.8 | | | -1.464 | 1.464 | |\n", "| 2 | Ctt | -0.5 | 1.0 | | | -0.5 | 0.5 | |\n", "| 3 | jpsi_p | -1.9 | 0.7 | | |-6.28319 | 6.28319 | |\n", "| 4 | jpsi_s | 0.980E4 | 0.011E4 | | | 9797.52 | 9996.48 | |\n", "| 5 | psi2s_s | 1430 | 50 | | | 1358.64 | 1433.36 | |\n", "---------------------------------------------------------------------------------------------\n", "-------------------------------------------------------------\n", "| | psi2s_p Dbar_s Ctt jpsi_p jpsi_s psi2s_s |\n", "-------------------------------------------------------------\n", "| psi2s_p | 1.000 0.896 0.021 0.977 -0.145 -0.023 |\n", "| Dbar_s | 0.896 1.000 0.236 0.941 -0.223 0.004 |\n", "| Ctt | 0.021 0.236 1.000 0.013 -0.053 -0.007 |\n", "| jpsi_p | 0.977 0.941 0.013 1.000 -0.133 -0.004 |\n", "| jpsi_s | -0.145 -0.223 -0.053 -0.133 1.000 -0.001 |\n", "| psi2s_s | -0.023 0.004 -0.007 -0.004 -0.001 1.000 |\n", "-------------------------------------------------------------\n", "Hesse errors: OrderedDict([(<zfit.Parameter 'psi2s_p' floating=True>, {'error': 0.41840024605493187}), (<zfit.Parameter 'Dbar_s' floating=True>, {'error': 1.8317143176053508}), (<zfit.Parameter 'Ctt' floating=True>, {'error': 0.996289918615876}), (<zfit.Parameter 'jpsi_p' floating=True>, {'error': 0.7473453412456892}), (<zfit.Parameter 'jpsi_s' floating=True>, {'error': 113.9126570143626}), (<zfit.Parameter 'psi2s_s' floating=True>, {'error': 47.11513567257623})])\n", "Toy 1/10\n", "Time taken: 2 min, 32 s\n", "Projected time left: 22 min, 50 s\n", "Toy 1: Generating data...\n", "Toy 1: Data generation finished\n", "Toy 1: Loading data...\n", "Toy 1: Loading data finished\n", "Toy 1: Fitting pdf...\n", "------------------------------------------------------------------\n", "| FCN = 4.145E+04 | Ncalls=174 (174 total) |\n", "| EDM = 0.000246 (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", "Toy 1: Fitting finished\n", "Function minimum: 41446.21391501494\n", "---------------------------------------------------------------------------------------------\n", "| | Name | Value | Hesse Err | Minos Err- | Minos Err+ | Limit- | Limit+ | Fixed |\n", "---------------------------------------------------------------------------------------------\n", "| 0 | psi2s_p | 1.51 | 0.17 | | |-6.28319 | 6.28319 | |\n", "| 1 | Dbar_s | -0.7 | 1.8 | | | -1.464 | 1.464 | |\n", "| 2 | Ctt | -0.5 | 0.8 | | | -0.5 | 0.5 | |\n", "| 3 | jpsi_p | 2.0 | 0.5 | | |-6.28319 | 6.28319 | |\n", "| 4 | jpsi_s | 0.980E4 | 0.013E4 | | | 9797.52 | 9996.48 | |\n", "| 5 | psi2s_s | 1430 | 50 | | | 1358.64 | 1433.36 | |\n", "---------------------------------------------------------------------------------------------\n", "-------------------------------------------------------------\n", "| | psi2s_p Dbar_s Ctt jpsi_p jpsi_s psi2s_s |\n", "-------------------------------------------------------------\n", "| psi2s_p | 1.00 0.80 0.84 0.72 0.16 -0.27 |\n", "| Dbar_s | 0.80 1.00 0.98 0.88 0.17 -0.22 |\n", "| Ctt | 0.84 0.98 1.00 0.94 0.21 -0.23 |\n", "| jpsi_p | 0.72 0.88 0.94 1.00 0.17 -0.22 |\n", "| jpsi_s | 0.16 0.17 0.21 0.17 1.00 -0.05 |\n", "| psi2s_s | -0.27 -0.22 -0.23 -0.22 -0.05 1.00 |\n", "-------------------------------------------------------------\n", "Hesse errors: OrderedDict([(<zfit.Parameter 'psi2s_p' floating=True>, {'error': 0.17027941762338372}), (<zfit.Parameter 'Dbar_s' floating=True>, {'error': 1.7969490420133767}), (<zfit.Parameter 'Ctt' floating=True>, {'error': 0.8149920548026126}), (<zfit.Parameter 'jpsi_p' floating=True>, {'error': 0.4859306435087767}), (<zfit.Parameter 'jpsi_s' floating=True>, {'error': 129.95009225722333}), (<zfit.Parameter 'psi2s_s' floating=True>, {'error': 47.350217126767575})])\n", "Toy 2/10\n", "Time taken: 5 min, 16 s\n", "Projected time left: 21 min, 7 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 = 4.145E+04 | Ncalls=204 (204 total) |\n", "| EDM = 0.00128 (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 | True | True | False |\n", "------------------------------------------------------------------\n", "Toy 2: Fitting finished\n", "Function minimum: 41446.21446623224\n", "---------------------------------------------------------------------------------------------\n", "| | Name | Value | Hesse Err | Minos Err- | Minos Err+ | Limit- | Limit+ | Fixed |\n", "---------------------------------------------------------------------------------------------\n", "| 0 | psi2s_p | -4.78 | 0.18 | | |-6.28319 | 6.28319 | |\n", "| 1 | Dbar_s | -0.7 | 1.8 | | | -1.464 | 1.464 | |\n", "| 2 | Ctt | -0.5 | 0.7 | | | -0.5 | 0.5 | |\n", "| 3 | jpsi_p | -4.3 | 0.5 | | |-6.28319 | 6.28319 | |\n", "| 4 | jpsi_s | 0.980E4 | 0.013E4 | | | 9797.52 | 9996.48 | |\n", "| 5 | psi2s_s | 1430 | 40 | | | 1358.64 | 1433.36 | |\n", "---------------------------------------------------------------------------------------------\n", "-------------------------------------------------------------\n", "| | psi2s_p Dbar_s Ctt jpsi_p jpsi_s psi2s_s |\n", "-------------------------------------------------------------\n", "| psi2s_p | 1.00 0.82 0.86 0.76 0.20 0.32 |\n", "| Dbar_s | 0.82 1.00 0.98 0.89 0.21 0.26 |\n", "| Ctt | 0.86 0.98 1.00 0.95 0.24 0.27 |\n", "| jpsi_p | 0.76 0.89 0.95 1.00 0.21 0.26 |\n", "| jpsi_s | 0.20 0.21 0.24 0.21 1.00 0.07 |\n", "| psi2s_s | 0.32 0.26 0.27 0.26 0.07 1.00 |\n", "-------------------------------------------------------------\n", "Hesse errors: OrderedDict([(<zfit.Parameter 'psi2s_p' floating=True>, {'error': 0.17884003769401335}), (<zfit.Parameter 'Dbar_s' floating=True>, {'error': 1.7953075848497315}), (<zfit.Parameter 'Ctt' floating=True>, {'error': 0.6507663138579277}), (<zfit.Parameter 'jpsi_p' floating=True>, {'error': 0.49871222470060794}), (<zfit.Parameter 'jpsi_s' floating=True>, {'error': 127.32075889710268}), (<zfit.Parameter 'psi2s_s' floating=True>, {'error': 38.958773985027506})])\n", "Toy 3/10\n", "Time taken: 7 min, 58 s\n", "Projected time left: 18 min, 37 s\n", "Toy 3: Generating data...\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Toy 3: Data generation finished\n", "Toy 3: Loading data...\n", "Toy 3: Loading data finished\n", "Toy 3: Fitting pdf...\n", "------------------------------------------------------------------\n", "| FCN = 4.145E+04 | Ncalls=233 (233 total) |\n", "| EDM = 0.000301 (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", "Toy 3: Fitting finished\n", "Function minimum: 41446.213143894536\n", "---------------------------------------------------------------------------------------------\n", "| | Name | Value | Hesse Err | Minos Err- | Minos Err+ | Limit- | Limit+ | Fixed |\n", "---------------------------------------------------------------------------------------------\n", "| 0 | psi2s_p | -4.77 | 0.20 | | |-6.28319 | 6.28319 | |\n", "| 1 | Dbar_s | -0.7 | 1.8 | | | -1.464 | 1.464 | |\n", "| 2 | Ctt | -0.5 | 1.0 | | | -0.5 | 0.5 | |\n", "| 3 | jpsi_p | -4.3 | 0.5 | | |-6.28319 | 6.28319 | |\n", "| 4 | jpsi_s | 0.980E4 | 0.012E4 | | | 9797.52 | 9996.48 | |\n", "| 5 | psi2s_s | 1430 | 50 | | | 1358.64 | 1433.36 | |\n", "---------------------------------------------------------------------------------------------\n", "-------------------------------------------------------------\n", "| | psi2s_p Dbar_s Ctt jpsi_p jpsi_s psi2s_s |\n", "-------------------------------------------------------------\n", "| psi2s_p | 1.00 0.89 0.92 0.85 0.41 -0.32 |\n", "| Dbar_s | 0.89 1.00 0.98 0.93 0.43 -0.26 |\n", "| Ctt | 0.92 0.98 1.00 0.97 0.46 -0.27 |\n", "| jpsi_p | 0.85 0.93 0.97 1.00 0.43 -0.26 |\n", "| jpsi_s | 0.41 0.43 0.46 0.43 1.00 -0.12 |\n", "| psi2s_s | -0.32 -0.26 -0.27 -0.26 -0.12 1.00 |\n", "-------------------------------------------------------------\n", "Hesse errors: OrderedDict([(<zfit.Parameter 'psi2s_p' floating=True>, {'error': 0.2028644492314955}), (<zfit.Parameter 'Dbar_s' floating=True>, {'error': 1.8154514362409575}), (<zfit.Parameter 'Ctt' floating=True>, {'error': 0.9807202234423524}), (<zfit.Parameter 'jpsi_p' floating=True>, {'error': 0.5159537418531039}), (<zfit.Parameter 'jpsi_s' floating=True>, {'error': 117.77974601317328}), (<zfit.Parameter 'psi2s_s' floating=True>, {'error': 45.73863717851896})])\n", "Toy 4/10\n", "Time taken: 10 min, 46 s\n", "Projected time left: 16 min, 9 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 = 4.145E+04 | Ncalls=180 (180 total) |\n", "| EDM = 0.000259 (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", "Toy 4: Fitting finished\n", "Function minimum: 41446.213230540496\n", "---------------------------------------------------------------------------------------------\n", "| | Name | Value | Hesse Err | Minos Err- | Minos Err+ | Limit- | Limit+ | Fixed |\n", "---------------------------------------------------------------------------------------------\n", "| 0 | psi2s_p | -4.77 | 0.21 | | |-6.28319 | 6.28319 | |\n", "| 1 | Dbar_s | -0.7 | 1.8 | | | -1.464 | 1.464 | |\n", "| 2 | Ctt | -0.5 | 0.7 | | | -0.5 | 0.5 | |\n", "| 3 | jpsi_p | -4.3 | 0.5 | | |-6.28319 | 6.28319 | |\n", "| 4 | jpsi_s | 0.980E4 | 0.012E4 | | | 9797.52 | 9996.48 | |\n", "| 5 | psi2s_s | 1430 | 60 | | | 1358.64 | 1433.36 | |\n", "---------------------------------------------------------------------------------------------\n", "-------------------------------------------------------------\n", "| | psi2s_p Dbar_s Ctt jpsi_p jpsi_s psi2s_s |\n", "-------------------------------------------------------------\n", "| psi2s_p | 1.00 0.89 0.92 0.87 0.45 -0.45 |\n", "| Dbar_s | 0.89 1.00 0.98 0.93 0.47 -0.37 |\n", "| Ctt | 0.92 0.98 1.00 0.97 0.50 -0.38 |\n", "| jpsi_p | 0.87 0.93 0.97 1.00 0.47 -0.38 |\n", "| jpsi_s | 0.45 0.47 0.50 0.47 1.00 -0.19 |\n", "| psi2s_s | -0.45 -0.37 -0.38 -0.38 -0.19 1.00 |\n", "-------------------------------------------------------------\n", "Hesse errors: OrderedDict([(<zfit.Parameter 'psi2s_p' floating=True>, {'error': 0.21090063108815427}), (<zfit.Parameter 'Dbar_s' floating=True>, {'error': 1.8277281246469832}), (<zfit.Parameter 'Ctt' floating=True>, {'error': 0.7144740345016457}), (<zfit.Parameter 'jpsi_p' floating=True>, {'error': 0.5130723061292268}), (<zfit.Parameter 'jpsi_s' floating=True>, {'error': 120.0500270283901}), (<zfit.Parameter 'psi2s_s' floating=True>, {'error': 59.20149657515435})])\n", "Toy 5/10\n", "Time taken: 13 min, 31 s\n", "Projected time left: 13 min, 31 s\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 = 4.145E+04 | Ncalls=304 (304 total) |\n", "| EDM = 0.00152 (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", "Toy 5: Fitting finished\n", "Function minimum: 41446.26376718663\n", "---------------------------------------------------------------------------------------------\n", "| | Name | Value | Hesse Err | Minos Err- | Minos Err+ | Limit- | Limit+ | Fixed |\n", "---------------------------------------------------------------------------------------------\n", "| 0 | psi2s_p | 4.67 | 0.32 | | |-6.28319 | 6.28319 | |\n", "| 1 | Dbar_s | -0.8 | 1.8 | | | -1.464 | 1.464 | |\n", "| 2 | Ctt | -0.5 | 0.5 | | | -0.5 | 0.5 | |\n", "| 3 | jpsi_p | -2.1 | 0.5 | | |-6.28319 | 6.28319 | |\n", "| 4 | jpsi_s | 0.980E4 | 0.015E4 | | | 9797.52 | 9996.48 | |\n", "| 5 | psi2s_s | 1430 | 50 | | | 1358.64 | 1433.36 | |\n", "---------------------------------------------------------------------------------------------\n", "-------------------------------------------------------------\n", "| | psi2s_p Dbar_s Ctt jpsi_p jpsi_s psi2s_s |\n", "-------------------------------------------------------------\n", "| psi2s_p | 1.000 -0.915 -0.954 0.850 -0.143 0.107 |\n", "| Dbar_s | -0.915 1.000 0.967 -0.848 0.117 -0.004 |\n", "| Ctt | -0.954 0.967 1.000 -0.937 0.157 -0.019 |\n", "| jpsi_p | 0.850 -0.848 -0.937 1.000 -0.126 0.030 |\n", "| jpsi_s | -0.143 0.117 0.157 -0.126 1.000 -0.005 |\n", "| psi2s_s | 0.107 -0.004 -0.019 0.030 -0.005 1.000 |\n", "-------------------------------------------------------------\n", "Hesse errors: OrderedDict([(<zfit.Parameter 'psi2s_p' floating=True>, {'error': 0.31586736018571226}), (<zfit.Parameter 'Dbar_s' floating=True>, {'error': 1.7931181136417407}), (<zfit.Parameter 'Ctt' floating=True>, {'error': 0.505275767707557}), (<zfit.Parameter 'jpsi_p' floating=True>, {'error': 0.5046726165950495}), (<zfit.Parameter 'jpsi_s' floating=True>, {'error': 147.58315697207308}), (<zfit.Parameter 'psi2s_s' floating=True>, {'error': 45.95130250092461})])\n", "Toy 6/10\n", "Time taken: 16 min, 29 s\n", "Projected time left: 10 min, 59 s\n", "Toy 6: Generating data...\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Toy 6: Data generation finished\n", "Toy 6: Loading data...\n", "Toy 6: Loading data finished\n", "Toy 6: Fitting pdf...\n", "------------------------------------------------------------------\n", "| FCN = 4.145E+04 | Ncalls=216 (216 total) |\n", "| EDM = 0.000883 (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", "Toy 6: Fitting finished\n", "Function minimum: 41446.21393607266\n", "---------------------------------------------------------------------------------------------\n", "| | Name | Value | Hesse Err | Minos Err- | Minos Err+ | Limit- | Limit+ | Fixed |\n", "---------------------------------------------------------------------------------------------\n", "| 0 | psi2s_p | -4.77 | 0.19 | | |-6.28319 | 6.28319 | |\n", "| 1 | Dbar_s | -0.7 | 1.8 | | | -1.464 | 1.464 | |\n", "| 2 | Ctt | -0.5 | 0.5 | | | -0.5 | 0.5 | |\n", "| 3 | jpsi_p | 2.0 | 0.5 | | |-6.28319 | 6.28319 | |\n", "| 4 | jpsi_s | 0.980E4 | 0.010E4 | | | 9797.52 | 9996.48 | |\n", "| 5 | psi2s_s | 1430 | 40 | | | 1358.64 | 1433.36 | |\n", "---------------------------------------------------------------------------------------------\n", "-------------------------------------------------------------\n", "| | psi2s_p Dbar_s Ctt jpsi_p jpsi_s psi2s_s |\n", "-------------------------------------------------------------\n", "| psi2s_p | 1.00 0.87 0.90 0.83 0.49 0.16 |\n", "| Dbar_s | 0.87 1.00 0.98 0.92 0.51 0.13 |\n", "| Ctt | 0.90 0.98 1.00 0.96 0.56 0.13 |\n", "| jpsi_p | 0.83 0.92 0.96 1.00 0.51 0.13 |\n", "| jpsi_s | 0.49 0.51 0.56 0.51 1.00 0.07 |\n", "| psi2s_s | 0.16 0.13 0.13 0.13 0.07 1.00 |\n", "-------------------------------------------------------------\n", "Hesse errors: OrderedDict([(<zfit.Parameter 'psi2s_p' floating=True>, {'error': 0.1917856359807728}), (<zfit.Parameter 'Dbar_s' floating=True>, {'error': 1.8156013659952255}), (<zfit.Parameter 'Ctt' floating=True>, {'error': 0.5117900295417959}), (<zfit.Parameter 'jpsi_p' floating=True>, {'error': 0.49818141138885563}), (<zfit.Parameter 'jpsi_s' floating=True>, {'error': 100.16204116151494}), (<zfit.Parameter 'psi2s_s' floating=True>, {'error': 38.14017204277275})])\n", "Toy 7/10\n", "Time taken: 19 min, 23 s\n", "Projected time left: 8 min, 18 s\n", "Toy 7: Generating data...\n", "Toy 7: Data generation finished\n", "Toy 7: Loading data...\n", "Toy 7: Loading data finished\n", "Toy 7: Fitting pdf...\n", "------------------------------------------------------------------\n", "| FCN = 4.145E+04 | Ncalls=202 (202 total) |\n", "| EDM = 0.000126 (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", "Toy 7: Fitting finished\n", "Function minimum: 41446.171485007675\n", "---------------------------------------------------------------------------------------------\n", "| | Name | Value | Hesse Err | Minos Err- | Minos Err+ | Limit- | Limit+ | Fixed |\n", "---------------------------------------------------------------------------------------------\n", "| 0 | psi2s_p | 1.9 | 0.4 | | |-6.28319 | 6.28319 | |\n", "| 1 | Dbar_s | -0.7 | 1.8 | | | -1.464 | 1.464 | |\n", "| 2 | Ctt | -0.5 | 1.0 | | | -0.5 | 0.5 | |\n", "| 3 | jpsi_p | -1.9 | 0.7 | | |-6.28319 | 6.28319 | |\n", "| 4 | jpsi_s | 0.980E4 | 0.012E4 | | | 9797.52 | 9996.48 | |\n", "| 5 | psi2s_s | 1430 | 60 | | | 1358.64 | 1433.36 | |\n", "---------------------------------------------------------------------------------------------\n", "-------------------------------------------------------------\n", "| | psi2s_p Dbar_s Ctt jpsi_p jpsi_s psi2s_s |\n", "-------------------------------------------------------------\n", "| psi2s_p | 1.000 0.861 0.035 0.970 -0.112 -0.115 |\n", "| Dbar_s | 0.861 1.000 0.322 0.916 -0.184 0.022 |\n", "| Ctt | 0.035 0.322 1.000 0.018 -0.058 -0.048 |\n", "| jpsi_p | 0.970 0.916 0.018 1.000 -0.102 -0.013 |\n", "| jpsi_s | -0.112 -0.184 -0.058 -0.102 1.000 -0.004 |\n", "| psi2s_s | -0.115 0.022 -0.048 -0.013 -0.004 1.000 |\n", "-------------------------------------------------------------\n", "Hesse errors: OrderedDict([(<zfit.Parameter 'psi2s_p' floating=True>, {'error': 0.39988882358765476}), (<zfit.Parameter 'Dbar_s' floating=True>, {'error': 1.8338950485749796}), (<zfit.Parameter 'Ctt' floating=True>, {'error': 0.9813742089966919}), (<zfit.Parameter 'jpsi_p' floating=True>, {'error': 0.709065375300536}), (<zfit.Parameter 'jpsi_s' floating=True>, {'error': 116.26239022795016}), (<zfit.Parameter 'psi2s_s' floating=True>, {'error': 64.86664680863748})])\n", "Toy 8/10\n", "Time taken: 22 min, 12 s\n", "Projected time left: 5 min, 33 s\n", "Toy 8: Generating data...\n", "Toy 8: Data generation finished\n", "Toy 8: Loading data...\n", "Toy 8: Loading data finished\n", "Toy 8: Fitting pdf...\n", "------------------------------------------------------------------\n", "| FCN = 4.145E+04 | Ncalls=185 (185 total) |\n", "| EDM = 7.47E-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", "Toy 8: Fitting finished\n", "Function minimum: 41446.171426607005\n", "---------------------------------------------------------------------------------------------\n", "| | Name | Value | Hesse Err | Minos Err- | Minos Err+ | Limit- | Limit+ | Fixed |\n", "---------------------------------------------------------------------------------------------\n", "| 0 | psi2s_p | -4.4 | 0.4 | | |-6.28319 | 6.28319 | |\n", "| 1 | Dbar_s | -0.7 | 1.8 | | | -1.464 | 1.464 | |\n", "| 2 | Ctt | -0.5 | 1.0 | | | -0.5 | 0.5 | |\n", "| 3 | jpsi_p | 4.4 | 0.7 | | |-6.28319 | 6.28319 | |\n", "| 4 | jpsi_s | 0.980E4 | 0.012E4 | | | 9797.52 | 9996.48 | |\n", "| 5 | psi2s_s | 1430 | 60 | | | 1358.64 | 1433.36 | |\n", "---------------------------------------------------------------------------------------------\n", "-------------------------------------------------------------\n", "| | psi2s_p Dbar_s Ctt jpsi_p jpsi_s psi2s_s |\n", "-------------------------------------------------------------\n", "| psi2s_p | 1.000 0.850 0.027 0.966 -0.108 -0.122 |\n", "| Dbar_s | 0.850 1.000 0.309 0.913 -0.184 0.021 |\n", "| Ctt | 0.027 0.309 1.000 0.012 -0.053 -0.042 |\n", "| jpsi_p | 0.966 0.913 0.012 1.000 -0.097 -0.019 |\n", "| jpsi_s | -0.108 -0.184 -0.053 -0.097 1.000 -0.004 |\n", "| psi2s_s | -0.122 0.021 -0.042 -0.019 -0.004 1.000 |\n", "-------------------------------------------------------------\n", "Hesse errors: OrderedDict([(<zfit.Parameter 'psi2s_p' floating=True>, {'error': 0.37581893654929166}), (<zfit.Parameter 'Dbar_s' floating=True>, {'error': 1.8374920961421453}), (<zfit.Parameter 'Ctt' floating=True>, {'error': 0.995543793163386}), (<zfit.Parameter 'jpsi_p' floating=True>, {'error': 0.6592821643495412}), (<zfit.Parameter 'jpsi_s' floating=True>, {'error': 118.1086035526505}), (<zfit.Parameter 'psi2s_s' floating=True>, {'error': 62.73818368176751})])\n", "Toy 9/10\n", "Time taken: 25 min, 2 s\n", "Projected time left: 2 min, 46 s\n", "Toy 9: Generating data...\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Toy 9: Data generation finished\n", "Toy 9: Loading data...\n", "Toy 9: Loading data finished\n", "Toy 9: Fitting pdf...\n", "------------------------------------------------------------------\n", "| FCN = 4.145E+04 | Ncalls=181 (181 total) |\n", "| EDM = 0.000129 (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", "Toy 9: Fitting finished\n", "Function minimum: 41446.27791906745\n", "---------------------------------------------------------------------------------------------\n", "| | Name | Value | Hesse Err | Minos Err- | Minos Err+ | Limit- | Limit+ | Fixed |\n", "---------------------------------------------------------------------------------------------\n", "| 0 | psi2s_p | -1.94 | 0.26 | | |-6.28319 | 6.28319 | |\n", "| 1 | Dbar_s | -0.7 | 0.8 | | | -1.464 | 1.464 | |\n", "| 2 | Ctt | -0.5 | 0.8 | | | -0.5 | 0.5 | |\n", "| 3 | jpsi_p | 1.8 | 0.4 | | |-6.28319 | 6.28319 | |\n", "| 4 | jpsi_s | 0.980E4 | 0.014E4 | | | 9797.52 | 9996.48 | |\n", "| 5 | psi2s_s | 1430 | 40 | | | 1358.64 | 1433.36 | |\n", "---------------------------------------------------------------------------------------------\n", "-------------------------------------------------------------\n", "| | psi2s_p Dbar_s Ctt jpsi_p jpsi_s psi2s_s |\n", "-------------------------------------------------------------\n", "| psi2s_p | 1.000 -0.766 -0.254 0.902 0.074 0.130 |\n", "| Dbar_s | -0.766 1.000 0.391 -0.817 -0.169 0.028 |\n", "| Ctt | -0.254 0.391 1.000 -0.005 -0.040 -0.028 |\n", "| jpsi_p | 0.902 -0.817 -0.005 1.000 0.041 0.019 |\n", "| jpsi_s | 0.074 -0.169 -0.040 0.041 1.000 -0.004 |\n", "| psi2s_s | 0.130 0.028 -0.028 0.019 -0.004 1.000 |\n", "-------------------------------------------------------------\n", "Hesse errors: OrderedDict([(<zfit.Parameter 'psi2s_p' floating=True>, {'error': 0.2574778505665698}), (<zfit.Parameter 'Dbar_s' floating=True>, {'error': 0.8041361051173593}), (<zfit.Parameter 'Ctt' floating=True>, {'error': 0.800071724878314}), (<zfit.Parameter 'jpsi_p' floating=True>, {'error': 0.40940340284378474}), (<zfit.Parameter 'jpsi_s' floating=True>, {'error': 136.70046081025066}), (<zfit.Parameter 'psi2s_s' floating=True>, {'error': 41.38771918974419})])\n", "Toy 10/10\n", "Time taken: 27 min, 51 s\n", "Projected time left: \n" ] } ], "source": [ "# zfit.run.numeric_checks = False \n", "\n", "Ctt_list = []\n", "Ctt_error_list = []\n", "\n", "nr_of_toys = 10\n", "nevents = int(pdg[\"number_of_decays\"])\n", "nevents = pdg[\"number_of_decays\"]\n", "event_stack = 1000000\n", "# zfit.settings.set_verbosity(10)\n", "calls = int(nevents/event_stack + 1)\n", "\n", "total_samp = []\n", "\n", "start = time.time()\n", "\n", "sampler = total_f.create_sampler(n=event_stack)\n", "\n", "for toy in range(nr_of_toys):\n", " \n", " ### Generate data\n", " \n", " 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", " 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", " with open(\"data/zfit_toys/toy_{0}/{1}.pkl\".format(toy, call), \"wb\") as f:\n", " pkl.dump(sam, f, pkl.HIGHEST_PROTOCOL)\n", " \n", " print(\"Toy {}: Data generation finished\".format(toy))\n", " \n", " ### Load data\n", " \n", " print(\"Toy {}: Loading data...\".format(toy))\n", "\n", " for call in range(calls):\n", " with open(r\"data/zfit_toys/toy_0/{}.pkl\".format(call), \"rb\") as input_file:\n", " sam = pkl.load(input_file)\n", " total_samp = np.append(total_samp, sam)\n", "\n", " total_samp = total_samp.astype('float64')\n", "\n", " 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 = (jpsi_mass+50.0, psi2s_mass-50.0))\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", " print(\"Toy {0}/{1}\".format(toy+1, nr_of_toys))\n", " print(\"Time taken: {}\".format(display_time(int(time.time() - start))))\n", " print(\"Projected time left: {}\".format(display_time(int((time.time() - start)/(c+calls*(toy))*((nr_of_toys-toy)*calls-c)))))\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Mean Ctt value = -0.4955320725102208\n", "Mean Ctt error = 0.7951298069508159\n" ] } ], "source": [ "print('Mean Ctt value = {}'.format(np.mean(Ctt_list)))\n", "print('Mean Ctt error = {}'.format(np.mean(Ctt_error_list)))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.3" } }, "nbformat": 4, "nbformat_minor": 2 }