{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Import"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\util\\execution.py:57: UserWarning: Not running on Linux. Determining available cpus for thread can failand be overestimated. Workaround (only if too many cpus are used):`zfit.run.set_n_cpu(your_cpu_number)`\n",
" warnings.warn(\"Not running on Linux. Determining available cpus for thread can fail\"\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"WARNING: The TensorFlow contrib module will not be included in TensorFlow 2.0.\n",
"For more information, please see:\n",
" * https://github.com/tensorflow/community/blob/master/rfcs/20180907-contrib-sunset.md\n",
" * https://github.com/tensorflow/addons\n",
"If you depend on functionality not listed there, please file an issue.\n",
"\n"
]
}
],
"source": [
"import os\n",
"\n",
"# os.environ[\"CUDA_VISIBLE_DEVICES\"] = \"-1\"\n",
"import random\n",
"import numpy as np\n",
"from pdg_const1 import pdg\n",
"import matplotlib\n",
"import matplotlib.pyplot as plt\n",
"import pickle as pkl\n",
"import sys\n",
"import time\n",
"from helperfunctions import display_time, prepare_plot\n",
"import cmath as c\n",
"import scipy.integrate as integrate\n",
"from scipy.optimize import fminbound\n",
"from array import array as arr\n",
"import collections\n",
"from itertools import compress\n",
"import tensorflow as tf\n",
"import zfit\n",
"from zfit import ztf\n",
"# from IPython.display import clear_output\n",
"import os\n",
"import tensorflow_probability as tfp\n",
"tfd = tfp.distributions\n",
"\n",
"from matplotlib.pyplot import figure"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# chunksize = 10000\n",
"# zfit.run.chunking.active = True\n",
"# zfit.run.chunking.max_n_points = chunksize"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Build model and graphs\n",
"## Create graphs"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"def formfactor(q2, subscript, b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2): #returns real value\n",
" #check if subscript is viable\n",
"\n",
" if subscript != \"0\" and subscript != \"+\" and subscript != \"T\":\n",
" raise ValueError('Wrong subscript entered, choose either 0, + or T')\n",
"\n",
" #get constants\n",
"\n",
" mK = ztf.constant(pdg['Ks_M'])\n",
" mbstar0 = ztf.constant(pdg[\"mbstar0\"])\n",
" mbstar = ztf.constant(pdg[\"mbstar\"])\n",
"\n",
"\n",
" mmu = ztf.constant(pdg['muon_M'])\n",
" mb = ztf.constant(pdg['bquark_M'])\n",
" ms = ztf.constant(pdg['squark_M'])\n",
" mB = ztf.constant(pdg['Bplus_M'])\n",
"\n",
" #N comes from derivation in paper\n",
"\n",
" N = 3\n",
"\n",
" #some helperfunctions\n",
"\n",
" tpos = (mB - mK)**2\n",
" tzero = (mB + mK)*(ztf.sqrt(mB)-ztf.sqrt(mK))**2\n",
"\n",
" z_oben = ztf.sqrt(tpos - q2) - ztf.sqrt(tpos - tzero)\n",
" z_unten = ztf.sqrt(tpos - q2) + ztf.sqrt(tpos - tzero)\n",
" z = tf.divide(z_oben, z_unten)\n",
"\n",
" #calculate f0\n",
"\n",
" if subscript == \"0\":\n",
" prefactor = 1/(1 - q2/(mbstar0**2))\n",
" _sum = 0\n",
" b0 = [b0_0, b0_1, b0_2]\n",
"\n",
" for i in range(N):\n",
" _sum += b0[i]*(tf.pow(z,i))\n",
"\n",
" return ztf.to_complex(prefactor * _sum)\n",
"\n",
" #calculate f+ or fT\n",
"\n",
" else:\n",
" prefactor = 1/(1 - q2/(mbstar**2))\n",
" _sum = 0\n",
"\n",
" if subscript == \"T\":\n",
" bT = [bT_0, bT_1, bT_2]\n",
" for i in range(N):\n",
" _sum += bT[i] * (tf.pow(z, i) - ((-1)**(i-N)) * (i/N) * tf.pow(z, N))\n",
" else:\n",
" bplus = [bplus_0, bplus_1, bplus_2]\n",
" for i in range(N):\n",
" _sum += bplus[i] * (tf.pow(z, i) - ((-1)**(i-N)) * (i/N) * tf.pow(z, N))\n",
"\n",
" return ztf.to_complex(prefactor * _sum)\n",
"\n",
"def resonance(q, _mass, width, phase, scale):\n",
"\n",
" q2 = tf.pow(q, 2)\n",
"\n",
" mmu = ztf.constant(pdg['muon_M'])\n",
"\n",
" p = 0.5 * ztf.sqrt(q2 - 4*(mmu**2))\n",
"\n",
" p0 = 0.5 * ztf.sqrt(_mass**2 - 4*mmu**2)\n",
"\n",
" gamma_j = tf.divide(p, q) * _mass * width / p0\n",
"\n",
" #Calculate the resonance\n",
"\n",
" _top = tf.complex(_mass * width, ztf.constant(0.0))\n",
"\n",
" _bottom = tf.complex(_mass**2 - q2, -_mass*gamma_j)\n",
"\n",
" com = _top/_bottom\n",
"\n",
" #Rotate by the phase\n",
"\n",
" r = ztf.to_complex(scale*tf.abs(com))\n",
"\n",
" _phase = tf.angle(com)\n",
"\n",
" _phase += phase\n",
"\n",
" com = r * tf.exp(tf.complex(ztf.constant(0.0), _phase))\n",
"\n",
" return com\n",
"\n",
"\n",
"def axiv_nonres(q, b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2):\n",
"\n",
" GF = ztf.constant(pdg['GF'])\n",
" alpha_ew = ztf.constant(pdg['alpha_ew'])\n",
" Vtb = ztf.constant(pdg['Vtb'])\n",
" Vts = ztf.constant(pdg['Vts'])\n",
" C10eff = ztf.constant(pdg['C10eff'])\n",
"\n",
" mmu = ztf.constant(pdg['muon_M'])\n",
" mb = ztf.constant(pdg['bquark_M'])\n",
" ms = ztf.constant(pdg['squark_M'])\n",
" mK = ztf.constant(pdg['Ks_M'])\n",
" mB = ztf.constant(pdg['Bplus_M'])\n",
"\n",
" q2 = tf.pow(q, 2)\n",
"\n",
" #Some helperfunctions\n",
"\n",
" beta = 1. - 4. * mmu**2. / q2\n",
"\n",
" kabs = ztf.sqrt(mB**2. + tf.pow(q2, 2)/mB**2. + mK**4./mB**2. - 2. * (mB**2. * mK**2. + mK**2. * q2 + mB**2. * q2) / mB**2.)\n",
"\n",
" #prefactor in front of whole bracket\n",
"\n",
" prefactor1 = GF**2. *alpha_ew**2. * (tf.abs(Vtb*Vts))**2. * kabs * beta / (128. * np.pi**5.)\n",
"\n",
" #left term in bracket\n",
"\n",
" bracket_left = 2./3. * tf.pow(kabs,2) * tf.pow(beta,2) * tf.pow(tf.abs(ztf.to_complex(C10eff)*formfactor(q2, \"+\", b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2)),2)\n",
"\n",
" #middle term in bracket\n",
"\n",
" _top = 4. * mmu**2. * (mB**2. - mK**2.) * (mB**2. - mK**2.)\n",
"\n",
" _under = q2 * mB**2.\n",
"\n",
" bracket_middle = _top/_under *tf.pow(tf.abs(ztf.to_complex(C10eff) * formfactor(q2, \"0\", b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2)), 2)\n",
" \n",
" #Note sqrt(q2) comes from derivation as we use q2 and plot q\n",
"\n",
" return prefactor1 * (bracket_left + bracket_middle) * 2 * q\n",
"\n",
"def vec(q, funcs, b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2):\n",
" \n",
" q2 = tf.pow(q, 2)\n",
"\n",
" GF = ztf.constant(pdg['GF'])\n",
" alpha_ew = ztf.constant(pdg['alpha_ew'])\n",
" Vtb = ztf.constant(pdg['Vtb'])\n",
" Vts = ztf.constant(pdg['Vts'])\n",
" C7eff = ztf.constant(pdg['C7eff'])\n",
"\n",
" mmu = ztf.constant(pdg['muon_M'])\n",
" mb = ztf.constant(pdg['bquark_M'])\n",
" ms = ztf.constant(pdg['squark_M'])\n",
" mK = ztf.constant(pdg['Ks_M'])\n",
" mB = ztf.constant(pdg['Bplus_M'])\n",
"\n",
" #Some helperfunctions\n",
"\n",
" beta = 1. - 4. * mmu**2. / q2\n",
"\n",
" kabs = ztf.sqrt(mB**2. + tf.pow(q2, 2)/mB**2. + mK**4./mB**2. - 2 * (mB**2 * mK**2 + mK**2 * q2 + mB**2 * q2) / mB**2)\n",
" \n",
" #prefactor in front of whole bracket\n",
"\n",
" prefactor1 = GF**2. *alpha_ew**2. * (tf.abs(Vtb*Vts))**2 * kabs * beta / (128. * np.pi**5.)\n",
"\n",
" #right term in bracket\n",
"\n",
" prefactor2 = tf.pow(kabs,2) * (1. - 1./3. * beta)\n",
"\n",
" abs_bracket = tf.pow(tf.abs(c9eff(q, funcs) * formfactor(q2, \"+\", b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2) + ztf.to_complex(2.0 * C7eff * (mb + ms)/(mB + mK)) * formfactor(q2, \"T\", b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2)),2)\n",
"\n",
" bracket_right = prefactor2 * abs_bracket\n",
"\n",
" #Note sqrt(q2) comes from derivation as we use q2 and plot q\n",
"\n",
" return prefactor1 * bracket_right * 2 * q\n",
"\n",
"def c9eff(q, funcs):\n",
"\n",
" C9eff_nr = ztf.to_complex(ztf.constant(pdg['C9eff']))\n",
"\n",
" c9 = C9eff_nr + funcs\n",
"\n",
" return c9"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def G(y):\n",
" \n",
" def inner_rect_bracket(q):\n",
" return tf.log(ztf.to_complex((1+tf.sqrt(q))/(1-tf.sqrt(q)))-tf.complex(ztf.constant(0), -1*ztf.constant(np.pi))) \n",
" \n",
" def inner_right(q):\n",
" return ztf.to_complex(2 * tf.atan(1/tf.sqrt(tf.math.real(-q))))\n",
" \n",
" big_bracket = tf.where(tf.math.real(y) > ztf.constant(0.0), inner_rect_bracket(y), inner_right(y))\n",
" \n",
" return ztf.to_complex(tf.sqrt(tf.abs(y))) * big_bracket\n",
"\n",
"def h_S(m, q):\n",
" \n",
" return ztf.to_complex(2) - G(ztf.to_complex(1) - ztf.to_complex(4*tf.pow(m, 2)) / ztf.to_complex(tf.pow(q, 2)))\n",
"\n",
"def h_P(m, q):\n",
" \n",
" return ztf.to_complex(2/3) + (ztf.to_complex(1) - ztf.to_complex(4*tf.pow(m, 2)) / ztf.to_complex(tf.pow(q, 2))) * h_S(m,q)\n",
"\n",
"def two_p_ccbar(mD, m_D_bar, m_D_star, q):\n",
" \n",
" \n",
" #Load constants\n",
" nu_D_bar = ztf.to_complex(pdg[\"nu_D_bar\"])\n",
" nu_D = ztf.to_complex(pdg[\"nu_D\"])\n",
" nu_D_star = ztf.to_complex(pdg[\"nu_D_star\"])\n",
" \n",
" phase_D_bar = ztf.to_complex(pdg[\"phase_D_bar\"])\n",
" phase_D = ztf.to_complex(pdg[\"phase_D\"])\n",
" phase_D_star = ztf.to_complex(pdg[\"phase_D_star\"])\n",
" \n",
" #Calculation\n",
" left_part = nu_D_bar * tf.exp(tf.complex(ztf.constant(0.0), phase_D_bar)) * h_S(m_D_bar, q) \n",
" \n",
" right_part_D = nu_D * tf.exp(tf.complex(ztf.constant(0.0), phase_D)) * h_P(m_D, q) \n",
" \n",
" right_part_D_star = nu_D_star * tf.exp(tf.complex(ztf.constant(0.0), phase_D_star)) * h_P(m_D_star, q) \n",
"\n",
" return left_part + right_part_D + right_part_D_star"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Build pdf"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"class total_pdf_cut(zfit.pdf.ZPDF):\n",
" _N_OBS = 1 # dimension, can be omitted\n",
" _PARAMS = ['b0_0', 'b0_1', 'b0_2', \n",
" 'bplus_0', 'bplus_1', 'bplus_2', \n",
" 'bT_0', 'bT_1', 'bT_2', \n",
" 'rho_mass', 'rho_scale', 'rho_phase', 'rho_width',\n",
" 'jpsi_mass', 'jpsi_scale', 'jpsi_phase', 'jpsi_width',\n",
" 'psi2s_mass', 'psi2s_scale', 'psi2s_phase', 'psi2s_width',\n",
" 'p3770_mass', 'p3770_scale', 'p3770_phase', 'p3770_width',\n",
" 'p4040_mass', 'p4040_scale', 'p4040_phase', 'p4040_width',\n",
" 'p4160_mass', 'p4160_scale', 'p4160_phase', 'p4160_width',\n",
" 'p4415_mass', 'p4415_scale', 'p4415_phase', 'p4415_width',\n",
" 'omega_mass', 'omega_scale', 'omega_phase', 'omega_width',\n",
" 'phi_mass', 'phi_scale', 'phi_phase', 'phi_width',\n",
" 'Dbar_mass', 'Dbar_scale', 'Dbar_phase',\n",
" 'Dstar_mass', 'DDstar_scale', 'DDstar_phase', 'D_mass',\n",
" 'tau_mass', 'C_tt']\n",
"# the name of the parameters\n",
"\n",
" def _unnormalized_pdf(self, x):\n",
" \n",
" x = x.unstack_x()\n",
" \n",
" b0 = [self.params['b0_0'], self.params['b0_1'], self.params['b0_2']]\n",
" bplus = [self.params['bplus_0'], self.params['bplus_1'], self.params['bplus_2']]\n",
" bT = [self.params['bT_0'], self.params['bT_1'], self.params['bT_2']]\n",
" \n",
" def rho_res(q):\n",
" return resonance(q, _mass = self.params['rho_mass'], scale = self.params['rho_scale'],\n",
" phase = self.params['rho_phase'], width = self.params['rho_width'])\n",
" \n",
" def omega_res(q):\n",
" return resonance(q, _mass = self.params['omega_mass'], scale = self.params['omega_scale'],\n",
" phase = self.params['omega_phase'], width = self.params['omega_width'])\n",
" \n",
" def phi_res(q):\n",
" return resonance(q, _mass = self.params['phi_mass'], scale = self.params['phi_scale'],\n",
" phase = self.params['phi_phase'], width = self.params['phi_width'])\n",
"\n",
" def jpsi_res(q):\n",
" return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['jpsi_mass'], 2)) * resonance(q, _mass = self.params['jpsi_mass'], \n",
" scale = self.params['jpsi_scale'],\n",
" phase = self.params['jpsi_phase'], \n",
" width = self.params['jpsi_width'])\n",
" def psi2s_res(q):\n",
" return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['psi2s_mass'], 2)) * resonance(q, _mass = self.params['psi2s_mass'], \n",
" scale = self.params['psi2s_scale'],\n",
" phase = self.params['psi2s_phase'], \n",
" width = self.params['psi2s_width'])\n",
" def p3770_res(q):\n",
" return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['p3770_mass'], 2)) * resonance(q, _mass = self.params['p3770_mass'], \n",
" scale = self.params['p3770_scale'],\n",
" phase = self.params['p3770_phase'], \n",
" width = self.params['p3770_width'])\n",
" \n",
" def p4040_res(q):\n",
" return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['p4040_mass'], 2)) * resonance(q, _mass = self.params['p4040_mass'], \n",
" scale = self.params['p4040_scale'],\n",
" phase = self.params['p4040_phase'], \n",
" width = self.params['p4040_width'])\n",
" \n",
" def p4160_res(q):\n",
" return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['p4160_mass'], 2)) * resonance(q, _mass = self.params['p4160_mass'], \n",
" scale = self.params['p4160_scale'],\n",
" phase = self.params['p4160_phase'], \n",
" width = self.params['p4160_width'])\n",
" \n",
" def p4415_res(q):\n",
" return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['p4415_mass'], 2)) * resonance(q, _mass = self.params['p4415_mass'], \n",
" scale = self.params['p4415_scale'],\n",
" phase = self.params['p4415_phase'], \n",
" width = self.params['p4415_width'])\n",
" \n",
" def P2_D(q):\n",
" Dbar_contrib = ztf.to_complex(self.params['Dbar_scale'])*tf.exp(tf.complex(ztf.constant(0.0), self.params['Dbar_phase']))*ztf.to_complex(h_S(self.params['Dbar_mass'], q))\n",
" DDstar_contrib = ztf.to_complex(self.params['DDstar_scale'])*tf.exp(tf.complex(ztf.constant(0.0), self.params['DDstar_phase']))*(ztf.to_complex(h_P(self.params['Dstar_mass'], q)) + ztf.to_complex(h_P(self.params['D_mass'], q)))\n",
" return Dbar_contrib + DDstar_contrib\n",
" \n",
" def ttau_cusp(q):\n",
" return ztf.to_complex(self.params['C_tt'])*(ztf.to_complex((h_S(self.params['tau_mass'], q))) - ztf.to_complex(h_P(self.params['tau_mass'], q)))\n",
" \n",
"\n",
" funcs = rho_res(x) + omega_res(x) + phi_res(x) + jpsi_res(x) + psi2s_res(x) + p3770_res(x) + p4040_res(x)+ p4160_res(x) + p4415_res(x) + P2_D(x) + ttau_cusp(x)\n",
"\n",
" vec_f = vec(x, funcs, b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2)\n",
"\n",
" axiv_nr = axiv_nonres(x, b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2)\n",
"\n",
" tot = vec_f + axiv_nr\n",
" \n",
" #Cut out jpsi and psi2s\n",
" \n",
" tot = tf.where(tf.math.logical_or(x < ztf.constant(jpsi_mass-60.), x > ztf.constant(jpsi_mass+70.)), tot, 0.0*tot)\n",
" \n",
" tot = tf.where(tf.math.logical_or(x < ztf.constant(psi2s_mass-50.), x > ztf.constant(psi2s_mass+50.)), tot, 0.0*tot)\n",
" \n",
" return tot\n",
" \n",
"class total_pdf_full(zfit.pdf.ZPDF):\n",
" _N_OBS = 1 # dimension, can be omitted\n",
" _PARAMS = ['b0_0', 'b0_1', 'b0_2', \n",
" 'bplus_0', 'bplus_1', 'bplus_2', \n",
" 'bT_0', 'bT_1', 'bT_2', \n",
" 'rho_mass', 'rho_scale', 'rho_phase', 'rho_width',\n",
" 'jpsi_mass', 'jpsi_scale', 'jpsi_phase', 'jpsi_width',\n",
" 'psi2s_mass', 'psi2s_scale', 'psi2s_phase', 'psi2s_width',\n",
" 'p3770_mass', 'p3770_scale', 'p3770_phase', 'p3770_width',\n",
" 'p4040_mass', 'p4040_scale', 'p4040_phase', 'p4040_width',\n",
" 'p4160_mass', 'p4160_scale', 'p4160_phase', 'p4160_width',\n",
" 'p4415_mass', 'p4415_scale', 'p4415_phase', 'p4415_width',\n",
" 'omega_mass', 'omega_scale', 'omega_phase', 'omega_width',\n",
" 'phi_mass', 'phi_scale', 'phi_phase', 'phi_width',\n",
" 'Dbar_mass', 'Dbar_scale', 'Dbar_phase',\n",
" 'Dstar_mass', 'DDstar_scale', 'DDstar_phase', 'D_mass',\n",
" 'tau_mass', 'C_tt']\n",
"# the name of the parameters\n",
"\n",
" def _unnormalized_pdf(self, x):\n",
" \n",
" x = x.unstack_x()\n",
" \n",
" b0 = [self.params['b0_0'], self.params['b0_1'], self.params['b0_2']]\n",
" bplus = [self.params['bplus_0'], self.params['bplus_1'], self.params['bplus_2']]\n",
" bT = [self.params['bT_0'], self.params['bT_1'], self.params['bT_2']]\n",
" \n",
" def rho_res(q):\n",
" return resonance(q, _mass = self.params['rho_mass'], scale = self.params['rho_scale'],\n",
" phase = self.params['rho_phase'], width = self.params['rho_width'])\n",
" \n",
" def omega_res(q):\n",
" return resonance(q, _mass = self.params['omega_mass'], scale = self.params['omega_scale'],\n",
" phase = self.params['omega_phase'], width = self.params['omega_width'])\n",
" \n",
" def phi_res(q):\n",
" return resonance(q, _mass = self.params['phi_mass'], scale = self.params['phi_scale'],\n",
" phase = self.params['phi_phase'], width = self.params['phi_width'])\n",
"\n",
" def jpsi_res(q):\n",
" return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['jpsi_mass'], 2)) * resonance(q, _mass = self.params['jpsi_mass'], \n",
" scale = self.params['jpsi_scale'],\n",
" phase = self.params['jpsi_phase'], \n",
" width = self.params['jpsi_width'])\n",
" def psi2s_res(q):\n",
" return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['psi2s_mass'], 2)) * resonance(q, _mass = self.params['psi2s_mass'], \n",
" scale = self.params['psi2s_scale'],\n",
" phase = self.params['psi2s_phase'], \n",
" width = self.params['psi2s_width'])\n",
" def p3770_res(q):\n",
" return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['p3770_mass'], 2)) * resonance(q, _mass = self.params['p3770_mass'], \n",
" scale = self.params['p3770_scale'],\n",
" phase = self.params['p3770_phase'], \n",
" width = self.params['p3770_width'])\n",
" \n",
" def p4040_res(q):\n",
" return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['p4040_mass'], 2)) * resonance(q, _mass = self.params['p4040_mass'], \n",
" scale = self.params['p4040_scale'],\n",
" phase = self.params['p4040_phase'], \n",
" width = self.params['p4040_width'])\n",
" \n",
" def p4160_res(q):\n",
" return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['p4160_mass'], 2)) * resonance(q, _mass = self.params['p4160_mass'], \n",
" scale = self.params['p4160_scale'],\n",
" phase = self.params['p4160_phase'], \n",
" width = self.params['p4160_width'])\n",
" \n",
" def p4415_res(q):\n",
" return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['p4415_mass'], 2)) * resonance(q, _mass = self.params['p4415_mass'], \n",
" scale = self.params['p4415_scale'],\n",
" phase = self.params['p4415_phase'], \n",
" width = self.params['p4415_width'])\n",
" \n",
" def P2_D(q):\n",
" Dbar_contrib = ztf.to_complex(self.params['Dbar_scale'])*tf.exp(tf.complex(ztf.constant(0.0), self.params['Dbar_phase']))*ztf.to_complex(h_S(self.params['Dbar_mass'], q))\n",
" DDstar_contrib = ztf.to_complex(self.params['DDstar_scale'])*tf.exp(tf.complex(ztf.constant(0.0), self.params['DDstar_phase']))*(ztf.to_complex(h_P(self.params['Dstar_mass'], q)) + ztf.to_complex(h_P(self.params['D_mass'], q)))\n",
" return Dbar_contrib + DDstar_contrib\n",
" \n",
" def ttau_cusp(q):\n",
" return ztf.to_complex(self.params['C_tt'])*(ztf.to_complex((h_S(self.params['tau_mass'], q))) - ztf.to_complex(h_P(self.params['tau_mass'], q)))\n",
" \n",
"\n",
" funcs = rho_res(x) + omega_res(x) + phi_res(x) + jpsi_res(x) + psi2s_res(x) + p3770_res(x) + p4040_res(x)+ p4160_res(x) + p4415_res(x) + P2_D(x) + ttau_cusp(x)\n",
"\n",
" vec_f = vec(x, funcs, b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2)\n",
"\n",
" axiv_nr = axiv_nonres(x, b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2)\n",
"\n",
" tot = vec_f + axiv_nr\n",
" \n",
" #Cut out jpsi and psi2s\n",
" \n",
"# tot = tf.where(tf.math.logical_or(x < ztf.constant(jpsi_mass-60.), x > ztf.constant(jpsi_mass+70.)), tot, 0.0*tot)\n",
" \n",
"# tot = tf.where(tf.math.logical_or(x < ztf.constant(psi2s_mass-50.), x > ztf.constant(psi2s_mass+50.)), tot, 0.0*tot)\n",
" \n",
" return tot"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup parameters"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"WARNING:tensorflow:From C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\ops\\resource_variable_ops.py:435: colocate_with (from tensorflow.python.framework.ops) is deprecated and will be removed in a future version.\n",
"Instructions for updating:\n",
"Colocations handled automatically by placer.\n"
]
}
],
"source": [
"# formfactors\n",
"\n",
"b0_0 = zfit.Parameter(\"b0_0\", ztf.constant(0.292), floating = False) #, lower_limit = -2.0, upper_limit= 2.0)\n",
"b0_1 = zfit.Parameter(\"b0_1\", ztf.constant(0.281), floating = False) #, lower_limit = -2.0, upper_limit= 2.0)\n",
"b0_2 = zfit.Parameter(\"b0_2\", ztf.constant(0.150), floating = False) #, lower_limit = -2.0, upper_limit= 2.0)\n",
"\n",
"bplus_0 = zfit.Parameter(\"bplus_0\", ztf.constant(0.466), lower_limit = -2.0, upper_limit= 2.0)\n",
"bplus_1 = zfit.Parameter(\"bplus_1\", ztf.constant(-0.885), lower_limit = -2.0, upper_limit= 2.0)\n",
"bplus_2 = zfit.Parameter(\"bplus_2\", ztf.constant(-0.213), lower_limit = -2.0, upper_limit= 2.0)\n",
"\n",
"bT_0 = zfit.Parameter(\"bT_0\", ztf.constant(0.460), floating = False) #, lower_limit = -2.0, upper_limit= 2.0)\n",
"bT_1 = zfit.Parameter(\"bT_1\", ztf.constant(-1.089), floating = False) #, lower_limit = -2.0, upper_limit= 2.0)\n",
"bT_2 = zfit.Parameter(\"bT_2\", ztf.constant(-1.114), floating = False) #, lower_limit = -2.0, upper_limit= 2.0)\n",
"\n",
"\n",
"#rho\n",
"\n",
"rho_mass, rho_width, rho_phase, rho_scale = pdg[\"rho\"]\n",
"\n",
"rho_m = zfit.Parameter(\"rho_m\", ztf.constant(rho_mass), floating = False) #lower_limit = rho_mass - rho_width, upper_limit = rho_mass + rho_width)\n",
"rho_w = zfit.Parameter(\"rho_w\", ztf.constant(rho_width), floating = False)\n",
"rho_p = zfit.Parameter(\"rho_p\", ztf.constant(rho_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"rho_s = zfit.Parameter(\"rho_s\", ztf.constant(rho_scale), lower_limit=rho_scale-np.sqrt(rho_scale), upper_limit=rho_scale+np.sqrt(rho_scale))\n",
"\n",
"#omega\n",
"\n",
"omega_mass, omega_width, omega_phase, omega_scale = pdg[\"omega\"]\n",
"\n",
"omega_m = zfit.Parameter(\"omega_m\", ztf.constant(omega_mass), floating = False)\n",
"omega_w = zfit.Parameter(\"omega_w\", ztf.constant(omega_width), floating = False)\n",
"omega_p = zfit.Parameter(\"omega_p\", ztf.constant(omega_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"omega_s = zfit.Parameter(\"omega_s\", ztf.constant(omega_scale), lower_limit=omega_scale-np.sqrt(omega_scale), upper_limit=omega_scale+np.sqrt(omega_scale))\n",
"\n",
"\n",
"#phi\n",
"\n",
"phi_mass, phi_width, phi_phase, phi_scale = pdg[\"phi\"]\n",
"\n",
"phi_m = zfit.Parameter(\"phi_m\", ztf.constant(phi_mass), floating = False)\n",
"phi_w = zfit.Parameter(\"phi_w\", ztf.constant(phi_width), floating = False)\n",
"phi_p = zfit.Parameter(\"phi_p\", ztf.constant(phi_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"phi_s = zfit.Parameter(\"phi_s\", ztf.constant(phi_scale), lower_limit=phi_scale-np.sqrt(phi_scale), upper_limit=phi_scale+np.sqrt(phi_scale))\n",
"\n",
"#jpsi\n",
"\n",
"jpsi_mass, jpsi_width, jpsi_phase, jpsi_scale = pdg[\"jpsi\"]\n",
"\n",
"jpsi_m = zfit.Parameter(\"jpsi_m\", ztf.constant(jpsi_mass), floating = False)\n",
"jpsi_w = zfit.Parameter(\"jpsi_w\", ztf.constant(jpsi_width), floating = False)\n",
"jpsi_p = zfit.Parameter(\"jpsi_p\", ztf.constant(jpsi_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"jpsi_s = zfit.Parameter(\"jpsi_s\", ztf.constant(jpsi_scale), floating = False) #, lower_limit=jpsi_scale-np.sqrt(jpsi_scale), upper_limit=jpsi_scale+np.sqrt(jpsi_scale))\n",
"\n",
"#psi2s\n",
"\n",
"psi2s_mass, psi2s_width, psi2s_phase, psi2s_scale = pdg[\"psi2s\"]\n",
"\n",
"psi2s_m = zfit.Parameter(\"psi2s_m\", ztf.constant(psi2s_mass), floating = False)\n",
"psi2s_w = zfit.Parameter(\"psi2s_w\", ztf.constant(psi2s_width), floating = False)\n",
"psi2s_p = zfit.Parameter(\"psi2s_p\", ztf.constant(psi2s_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"psi2s_s = zfit.Parameter(\"psi2s_s\", ztf.constant(psi2s_scale), floating = False) #, lower_limit=psi2s_scale-np.sqrt(psi2s_scale), upper_limit=psi2s_scale+np.sqrt(psi2s_scale))\n",
"\n",
"#psi(3770)\n",
"\n",
"p3770_mass, p3770_width, p3770_phase, p3770_scale = pdg[\"p3770\"]\n",
"\n",
"p3770_m = zfit.Parameter(\"p3770_m\", ztf.constant(p3770_mass), floating = False)\n",
"p3770_w = zfit.Parameter(\"p3770_w\", ztf.constant(p3770_width), floating = False)\n",
"p3770_p = zfit.Parameter(\"p3770_p\", ztf.constant(p3770_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"p3770_s = zfit.Parameter(\"p3770_s\", ztf.constant(p3770_scale), lower_limit=p3770_scale-np.sqrt(p3770_scale), upper_limit=p3770_scale+np.sqrt(p3770_scale))\n",
"\n",
"#psi(4040)\n",
"\n",
"p4040_mass, p4040_width, p4040_phase, p4040_scale = pdg[\"p4040\"]\n",
"\n",
"p4040_m = zfit.Parameter(\"p4040_m\", ztf.constant(p4040_mass), floating = False)\n",
"p4040_w = zfit.Parameter(\"p4040_w\", ztf.constant(p4040_width), floating = False)\n",
"p4040_p = zfit.Parameter(\"p4040_p\", ztf.constant(p4040_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"p4040_s = zfit.Parameter(\"p4040_s\", ztf.constant(p4040_scale), lower_limit=p4040_scale-np.sqrt(p4040_scale), upper_limit=p4040_scale+np.sqrt(p4040_scale))\n",
"\n",
"#psi(4160)\n",
"\n",
"p4160_mass, p4160_width, p4160_phase, p4160_scale = pdg[\"p4160\"]\n",
"\n",
"p4160_m = zfit.Parameter(\"p4160_m\", ztf.constant(p4160_mass), floating = False)\n",
"p4160_w = zfit.Parameter(\"p4160_w\", ztf.constant(p4160_width), floating = False)\n",
"p4160_p = zfit.Parameter(\"p4160_p\", ztf.constant(p4160_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"p4160_s = zfit.Parameter(\"p4160_s\", ztf.constant(p4160_scale), lower_limit=p4160_scale-np.sqrt(p4160_scale), upper_limit=p4160_scale+np.sqrt(p4160_scale))\n",
"\n",
"#psi(4415)\n",
"\n",
"p4415_mass, p4415_width, p4415_phase, p4415_scale = pdg[\"p4415\"]\n",
"\n",
"p4415_m = zfit.Parameter(\"p4415_m\", ztf.constant(p4415_mass), floating = False)\n",
"p4415_w = zfit.Parameter(\"p4415_w\", ztf.constant(p4415_width), floating = False)\n",
"p4415_p = zfit.Parameter(\"p4415_p\", ztf.constant(p4415_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)\n",
"p4415_s = zfit.Parameter(\"p4415_s\", ztf.constant(p4415_scale), lower_limit=p4415_scale-np.sqrt(p4415_scale), upper_limit=p4415_scale+np.sqrt(p4415_scale))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Dynamic generation of 2 particle contribution"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"m_c = 1300\n",
"\n",
"Dbar_phase = 0.0\n",
"DDstar_phase = 0.0\n",
"Dstar_mass = pdg['Dst_M']\n",
"Dbar_mass = pdg['D0_M']\n",
"D_mass = pdg['D0_M']\n",
"\n",
"Dbar_s = zfit.Parameter(\"Dbar_s\", ztf.constant(0.0), lower_limit=-0.3, upper_limit=0.3)\n",
"Dbar_m = zfit.Parameter(\"Dbar_m\", ztf.constant(Dbar_mass), floating = False)\n",
"Dbar_p = zfit.Parameter(\"Dbar_p\", ztf.constant(Dbar_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)#, floating = False)\n",
"DDstar_s = zfit.Parameter(\"DDstar_s\", ztf.constant(0.0), lower_limit=-0.3, upper_limit=0.3)#, floating = False)\n",
"Dstar_m = zfit.Parameter(\"Dstar_m\", ztf.constant(Dstar_mass), floating = False)\n",
"D_m = zfit.Parameter(\"D_m\", ztf.constant(D_mass), floating = False)\n",
"DDstar_p = zfit.Parameter(\"DDstar_p\", ztf.constant(DDstar_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)#, floating = False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Tau parameters"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"tau_m = zfit.Parameter(\"tau_m\", ztf.constant(pdg['tau_M']), floating = False)\n",
"Ctt = zfit.Parameter(\"Ctt\", ztf.constant(0.0), lower_limit=-2.5, upper_limit=2.5)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Load data"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"x_min = 2*pdg['muon_M']\n",
"x_max = (pdg[\"Bplus_M\"]-pdg[\"Ks_M\"]-0.1)\n",
"\n",
"# # Full spectrum\n",
"\n",
"obs_toy = zfit.Space('q', limits = (x_min, x_max))\n",
"\n",
"# Jpsi and Psi2s cut out\n",
"\n",
"obs1 = zfit.Space('q', limits = (x_min, jpsi_mass - 60.))\n",
"obs2 = zfit.Space('q', limits = (jpsi_mass + 70., psi2s_mass - 50.))\n",
"obs3 = zfit.Space('q', limits = (psi2s_mass + 50., x_max))\n",
"\n",
"obs_fit = obs1 + obs2 + obs3\n",
"\n",
"# with open(r\"./data/slim_points/slim_points_toy_0_range({0}-{1}).pkl\".format(int(x_min), int(x_max)), \"rb\") as input_file:\n",
"# part_set = pkl.load(input_file)\n",
"\n",
"# x_part = part_set['x_part']\n",
"\n",
"# x_part = x_part.astype('float64')\n",
"\n",
"# data = zfit.data.Data.from_numpy(array=x_part, obs=obs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup pdf"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"total_f = total_pdf_cut(obs=obs_toy, jpsi_mass = jpsi_m, jpsi_scale = jpsi_s, jpsi_phase = jpsi_p, jpsi_width = jpsi_w,\n",
" psi2s_mass = psi2s_m, psi2s_scale = psi2s_s, psi2s_phase = psi2s_p, psi2s_width = psi2s_w,\n",
" p3770_mass = p3770_m, p3770_scale = p3770_s, p3770_phase = p3770_p, p3770_width = p3770_w,\n",
" p4040_mass = p4040_m, p4040_scale = p4040_s, p4040_phase = p4040_p, p4040_width = p4040_w,\n",
" p4160_mass = p4160_m, p4160_scale = p4160_s, p4160_phase = p4160_p, p4160_width = p4160_w,\n",
" p4415_mass = p4415_m, p4415_scale = p4415_s, p4415_phase = p4415_p, p4415_width = p4415_w,\n",
" rho_mass = rho_m, rho_scale = rho_s, rho_phase = rho_p, rho_width = rho_w,\n",
" omega_mass = omega_m, omega_scale = omega_s, omega_phase = omega_p, omega_width = omega_w,\n",
" phi_mass = phi_m, phi_scale = phi_s, phi_phase = phi_p, phi_width = phi_w,\n",
" Dstar_mass = Dstar_m, DDstar_scale = DDstar_s, DDstar_phase = DDstar_p, D_mass = D_m,\n",
" Dbar_mass = Dbar_m, Dbar_scale = Dbar_s, Dbar_phase = Dbar_p,\n",
" tau_mass = tau_m, C_tt = Ctt, b0_0 = b0_0, b0_1 = b0_1, b0_2 = b0_2,\n",
" bplus_0 = bplus_0, bplus_1 = bplus_1, bplus_2 = bplus_2,\n",
" bT_0 = bT_0, bT_1 = bT_1, bT_2 = bT_2)\n",
"\n",
"total_f_fit = total_pdf_full(obs=obs_fit, jpsi_mass = jpsi_m, jpsi_scale = jpsi_s, jpsi_phase = jpsi_p, jpsi_width = jpsi_w,\n",
" psi2s_mass = psi2s_m, psi2s_scale = psi2s_s, psi2s_phase = psi2s_p, psi2s_width = psi2s_w,\n",
" p3770_mass = p3770_m, p3770_scale = p3770_s, p3770_phase = p3770_p, p3770_width = p3770_w,\n",
" p4040_mass = p4040_m, p4040_scale = p4040_s, p4040_phase = p4040_p, p4040_width = p4040_w,\n",
" p4160_mass = p4160_m, p4160_scale = p4160_s, p4160_phase = p4160_p, p4160_width = p4160_w,\n",
" p4415_mass = p4415_m, p4415_scale = p4415_s, p4415_phase = p4415_p, p4415_width = p4415_w,\n",
" rho_mass = rho_m, rho_scale = rho_s, rho_phase = rho_p, rho_width = rho_w,\n",
" omega_mass = omega_m, omega_scale = omega_s, omega_phase = omega_p, omega_width = omega_w,\n",
" phi_mass = phi_m, phi_scale = phi_s, phi_phase = phi_p, phi_width = phi_w,\n",
" Dstar_mass = Dstar_m, DDstar_scale = DDstar_s, DDstar_phase = DDstar_p, D_mass = D_m,\n",
" Dbar_mass = Dbar_m, Dbar_scale = Dbar_s, Dbar_phase = Dbar_p,\n",
" tau_mass = tau_m, C_tt = Ctt, b0_0 = b0_0, b0_1 = b0_1, b0_2 = b0_2,\n",
" bplus_0 = bplus_0, bplus_1 = bplus_1, bplus_2 = bplus_2,\n",
" bT_0 = bT_0, bT_1 = bT_1, bT_2 = bT_2)\n",
" \n",
"# print(total_pdf.obs)\n",
"\n",
"# print(calcs_test)\n",
"\n",
"# for param in total_f.get_dependents():\n",
"# print(zfit.run(param))"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"# total_f_fit.normalization(obs_fit)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Test if graphs actually work and compute values"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"# def total_test_tf(xq):\n",
"\n",
"# def jpsi_res(q):\n",
"# return resonance(q, jpsi_m, jpsi_s, jpsi_p, jpsi_w)\n",
"\n",
"# def psi2s_res(q):\n",
"# return resonance(q, psi2s_m, psi2s_s, psi2s_p, psi2s_w)\n",
"\n",
"# def cusp(q):\n",
"# return bifur_gauss(q, cusp_m, sig_L, sig_R, cusp_s)\n",
"\n",
"# funcs = jpsi_res(xq) + psi2s_res(xq) + cusp(xq)\n",
"\n",
"# vec_f = vec(xq, funcs)\n",
"\n",
"# axiv_nr = axiv_nonres(xq)\n",
"\n",
"# tot = vec_f + axiv_nr\n",
" \n",
"# return tot\n",
"\n",
"# def jpsi_res(q):\n",
"# return resonance(q, jpsi_m, jpsi_s, jpsi_p, jpsi_w)\n",
"\n",
"# calcs = zfit.run(total_test_tf(x_part))\n",
"\n",
"\n",
"\n",
"test_q = np.linspace(x_min, x_max, int(2e6))\n",
"\n",
"probs = total_f_fit.pdf(test_q, norm_range=False)\n",
"\n",
"calcs_test = zfit.run(probs)\n",
"\n",
"Ctt.set_value(0.5)\n",
"\n",
"calcs_test1 = zfit.run(probs)\n",
"\n",
"Ctt.set_value(0.0)\n",
"\n",
"Dbar_s.set_value(0.3)\n",
"\n",
"DDstar_s.set_value(0.3)\n",
"\n",
"calcs_test2 = zfit.run(probs)\n",
"# res_y = zfit.run(jpsi_res(test_q))\n",
"# b0 = [b0_0, b0_1, b0_2]\n",
"# bplus = [bplus_0, bplus_1, bplus_2]\n",
"# bT = [bT_0, bT_1, bT_2]\n",
"# f0_y = zfit.run(tf.math.real(formfactor(test_q,\"0\", b0, bplus, bT)))\n",
"# fplus_y = zfit.run(tf.math.real(formfactor(test_q,\"+\", b0, bplus, bT)))\n",
"# fT_y = zfit.run(tf.math.real(formfactor(test_q,\"T\", b0, bplus, bT)))"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEMCAYAAADNtWEcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAciUlEQVR4nO3dfZAdV3km8OexZBsnTGxA4xSy5MhsJIwS2zG6ZVzF1iIHArIZSVTKm1gigSSGiQkmH1aFmCJlFOefxcT5YNeETIACUiuMIWyQp0R5XcRaUgEb3wnGWBYSY2Hj2fGiwTbOJCE4Ru/+0T0z7ftxTuvqnHvOTD+/KtX0Pfedft/b0tV7u2/3aZoZRESkeU5LXYCIiKShBiAi0lBqACIiDaUGICLSUGoAIiINpQYgItJQq1MlXrNmjW3YsCFVehGRZWlqaup7ZjYaYl3JGsCGDRvQbrdTpRdZPo4eBTZtOvWYWFy5U9a1QpF8LNS6dAhIJHdjY2FiYnHlTlmXeKkBiOTu6NEwMbG4cqesS7zUAERyt3dvmJhYXLlT1iVeagAiIg3FVJPBtVot05fAIiInh+SUmbVCrEt7ACK5q3MWTcozbVy5dQZQ1tQARHI3Ofm8h4ef+Gd07bl3xAyVK3fKusRLDUAkd/Pzi4v3fPM4rvyLf8Bn2jN9Y4bOlTtlXeKlBiCSu/HxxcVH5v4FAPDN/zffN2boXLlT1iVeagAiuZuaChMTiyt3yrrESw1AJHd79oSJicWVO2Vd4qUGIJK7tWvDxMTiyp2yLvFSAxDJnfYAJBI1AJHcaQ9AIlEDEMldnSvmU15V78qtq/2zpgYgkrseM2oazBszNJoNdNlSAxDJXWVGTZLemKHTbKDLlhqASO4OHgwTE4srd8q6xEsNQCR3da6m1ZXAMgBvAyD5MZLHST7U5/k3k3yw/PNlkpeEL1OkwVo1Zv6tExOLK3fKusSrzh7AxwFsczz/bQCvMbOLAfwxgIkAdYnIAu0BSCTeBmBmXwLwlOP5L5vZ0+XDewGsC1SbiADAyEiYmFhcuVPWJV6hvwO4FsAX+j1Jcpxkm2R7bm4ucGqRFWp2tmuo60Z+PWKGxpU7ZV3iFawBkLwCRQP4g34xZjZhZi0za42OjoZKLbKyVc6k6XMSqM4CkoEEaQAkLwbwEQA7zezJEOsUkdJEja/V6sTE4sqdsi7xOuUGQPJ8AJ8D8Ktmpsv+REK7884wMbG4cqesS7zqnAb6KQBfAfBykjMkryV5HcnrypCbALwEwIdIPkBSk3+IhLR7d5iYWFy5U9YlXqt9AWa2y/P82wC8LVhFIvJ8Y2NhYmJx5U5Zl3jpSmCR3GkPQCJRAxDJXb8J4E42JhZX7pR1iZcagEjuuk76HzAmFlfulHWJlxqASO727Vtc7PuBuhIzdK7cKesSLzUAkdxNToaJicWVO2Vd4qUGIJK7Op+itQcgA1ADEMnd9u1hYmJx5U5Zl3ipAYjkTtNBSyRqACK527o1TEwsrtwp6xIvNQCR3K1d2zVknadX9ogZGlfulHWJlxqASO7m5xcX+15WVYkZOlfulHWJlxqASO40HbREogYgkrt2jQl268TE4sqdsi7xUgMQyZ32ACQSNQCR3OksIIlEDUAkd3v3homJxZU7ZV3ipQYgkrtNm8LExOLKnbIu8VIDEMldq9U11DXJco+YoXHlTlmXeKkBiORudnZxkf3mg67EDJ0rd8q6xEsNQCR3t94aJiYWV+6UdYmXGoBI7up8itYegAzA2wBIfozkcZIP9XmeJD9IcprkgyRfGb5MkQbTHoBEUmcP4OMAtjmevxLAxvLPOIC/PPWyRGTRli1hYmJx5U5Zl3h5G4CZfQnAU46QnQA+aYV7AZxD8qWhChRpPF0JLJGE+A7gPACPVx7PlGMiEsLISJiYWFy5U9YlXiEaQK/z0rpOUwYAkuMk2yTbc3NzAVKLNMDYWNdQ5+0AesUMjSt3yrrEK0QDmAGwvvJ4HYCeX/2b2YSZtcysNTo6GiC1SAMcPbq42O8ygGrM0Llyp6xLvEI0gP0A3lKeDXQ5gGfM7IkA6xURQHMBSTSrfQEkPwVgK4A1JGcAvA/A6QBgZh8GcADAVQCmAfwbgF+PVayIiITjbQBmtsvzvAF4Z7CKROT5tAcgkehKYJHcaTZQiUQNQCR3k5NhYmJx5U5Zl3ipAYjkbn6+a8g6z7TuETM0rtwp6xIvNQCR3I2PLy72Owu0GjN0rtwp6xIvNQCR3E1NhYmJxZU7ZV3ipQYgkrs9e8LExOLKnbIu8VIDEMnd2rVhYmJx5U5Zl3ipAYjkTnsAEokagEjutAcgkagBiOSu3Q4TE4srd8q6xEsNQCR3PWbU7JoOWrOBygDUAERyV51Pp9980JoLSAagBiCSu4MHw8TE4sqdsi7xUgMQyV2dq2l1JbAMQA1AJHetVpiYWFy5U9YlXmoAIrnTHoBEogYgkruRkTAxsbhyp6xLvNQARHI3O9s11HkWaK+YoXHlTlmXeKkBiOSuciZN3+mgdRaQDEANQCR3ExNhYmJx5U5Zl3ipAYjk7s47w8TE4sqdsi7xUgMQyd3u3WFiYnHlTlmXeNVqACS3kTxCcprkjT2eP5/kPSS/RvJBkleFL1WkocbGwsTE4sqdsi7x8jYAkqsA3AbgSgCbAewiubkj7A8B3GFmlwK4BsCHQhcq0ljaA5BI6uwBXAZg2syOmdmzAG4HsLMjxgD8RLl8NgCd+yUSSr8J4E42JhZX7pR1iVedBnAegMcrj2fKsaq9AH6F5AyAAwDe1WtFJMdJtkm25+bmBihXpIG65n7uMdQjZmhcuVPWJV51GkCvFt75t7oLwMfNbB2AqwD8DcmudZvZhJm1zKw1Ojp68tWKNNG+fYuLfT9QV2KGzpU7ZV3iVacBzABYX3m8Dt2HeK4FcAcAmNlXALwAwJoQBYo03uRkmJhYXLlT1iVedRrA/QA2kryA5BkovuTd3xHzHQCvBQCSr0DRAHSMRySEOp+itQcgA/A2ADN7DsD1AO4CcBjF2T6HSN5MckcZtgfA20l+HcCnAPyamQ7+iQSxfXuYmFhcuVPWJV6r6wSZ2QEUX+5Wx26qLD8M4NVhSxMRAJoOWqLRlcAiudu6NUxMLK7cKesSLzUAkdytXdtj0GrEDIkrd8q6xEsNQCR38/OLi+w3IXQlZuhcuVPWJV5qACK5q0ypbN23gumKGTpNB71sqQGI5K7d7jHIGjFD4sqdsi7xUgMQyZ1uCCORqAGI5E5nAUkkagAiudu7N0xMLK7cKesSLzUAkdxt2hQmJhZX7pR1iZcagEjuWq0eg1YjZkhcuVPWJV5qACK5m12afLfvdQCzCe/B5Mqdsi7xUgMQyd2tt4aJicWVO2Vd4qUGIJK7Op+itQcgA1ADEMmd9gAkEjUAkdxt2RImJhZX7pR1iZcagEjudCWwRKIGIJK7kZEwMbG4cqesS7zUAERyNzbWNdR1w9UeMUPjyp2yLvFSAxDJ3dGji4vscxlANWboXLlT1iVeagAiudNcQBKJGoCISEPVagAkt5E8QnKa5I19Yn6J5MMkD5HcF7ZMkQbTHoBE4m0AJFcBuA3AlQA2A9hFcnNHzEYA7wHwajP7GQC/G6FWkWbSbKASSZ09gMsATJvZMTN7FsDtAHZ2xLwdwG1m9jQAmNnxsGWKNNjkZJiYWFy5U9YlXnUawHkAHq88ninHqjYB2ETyH0neS3JbqAJFGm9+vmuo6zTQHjFD48qdsi7xqtMAep141vnPbzWAjQC2AtgF4CMkz+laETlOsk2yPTc3d7K1ijTT+PjiYr+zQKsxQ+fKnbIu8arTAGYArK88Xgegc4q/GQCfN7P/MLNvAziCoiE8j5lNmFnLzFqjo6OD1izSLFNTYWJiceVOWZd41WkA9wPYSPICkmcAuAbA/o6YvwNwBQCQXIPikNCxkIWKNNaePWFiYnHlTlmXeHkbgJk9B+B6AHcBOAzgDjM7RPJmkjvKsLsAPEnyYQD3APh9M3syVtEijbJ2bZiYWFy5U9YlXrSub5OGo9VqWbvdTpJbZLm6/avfwY2f+wZ+ubUe77/64tTlSAIkp8wsyM2WdSWwSO60ByCRqAGI5K7OnnLKvWlXbu3lZ00NQCR3PWbUtM4zsTUbqAxADUAkd5X5dPpOB625gGQAagAiuTt4MExMLK7cKesSLzUAkdzVuZpWVwLLANQARHLXqnHGX52YWFy5U9YlXmoAIrnTHoBEogYgkruRkTAxsbhyp6xLvNQARHI32zn3Yo/poHvEDI0rd8q6xEsNQCR3lTNp2G9CaJ0FJANQAxDJ3cREmJhYXLlT1iVeagAiubvzzjAxsbhyp6xLvNQARHK3e3eYmFhcuVPWJV5qACK5GxsLExOLK3fKusRLDUAkd9oDkEjUAERy13cGuJOMicWVO2Vd4qUGIJK7Hnft6xpJdGc/b+6UdYmXGoBI7vbtW1ru94G6GjNsrtwp6xIvNQCR3E1OhomJxZU7ZV3ipQYgkrs6n6K1ByADUAMQyd327WFiYnHlTlmXeNVqACS3kTxCcprkjY64q0kaSU0CLhKKpoOWSLwNgOQqALcBuBLAZgC7SG7uETcC4LcB3Be6SJFG27o1TEwsrtwp6xKvOnsAlwGYNrNjZvYsgNsB7OwR98cAbgHw7wHrE5G1a8PExOLKnbIu8arTAM4D8Hjl8Uw5tojkpQDWm5nzK3+S4yTbJNtzc3MnXaxII83Pdw11nV7fI2ZoXLlT1iVedRpArzOPF//5kTwNwJ8B2ONbkZlNmFnLzFqjo6P1qxRpssqUyn2vq9V00DKAOg1gBsD6yuN1AKq3+RkB8LMADpJ8FMDlAPbri2CRQNrtMDGxuHKnrEu86jSA+wFsJHkByTMAXANg/8KTZvaMma0xsw1mtgHAvQB2mJn+5kVC0A1hJBJvAzCz5wBcD+AuAIcB3GFmh0jeTHJH7AJFGk9nAUkkq+sEmdkBAAc6xm7qE7v11MsSkUV794aJicWVO2Vd4qUrgUVyt2lTmJhYXLlT1iVeagAiuWt1n09hnRNC94gZGlfulHWJlxqASO5ml066Y78brMzO9h4fBlfulHWJlxqASO5uvTVMTCyu3CnrEi81AJHc1fkUrT0AGYAagEjutAcgkagBiORuy5YwMbG4cqesS7zUAERypyuBJRI1AJHcjYyEiYnFlTtlXeKlBiCSu7Gx7rHO6aB7xQyLK3fKusRLDUAkd0ePLi72nQ66EjN0rtwp6xIvNQCR3GkuIIlEDUBEpKHUAERypz0AiUQNQCR3mg1UIlEDEMnd5GSYmFhcuVPWJV5qACK5m5/vGuo8C7RXzNC4cqesS7zUAERyNz6+uNhvNuhqzNC5cqesS7zUAERyNzUVJiYWV+6UdYmXGoBI7vbsCRMTiyt3yrrESw1AJHdr14aJicWVO2Vd4lWrAZDcRvIIyWmSN/Z4/gaSD5N8kOQXSf5U+FJFGkp7ABKJtwGQXAXgNgBXAtgMYBfJzR1hXwPQMrOLAXwWwC2hCxVpLO0BSCR19gAuAzBtZsfM7FkAtwPYWQ0ws3vM7N/Kh/cCWBe2TJEGa7fDxMTiyp2yLvGq0wDOA/B45fFMOdbPtQC+cCpFiUhFjxk1zcwbMzSaDXTZqtMAep153HUdCgCQ/BUALQAf6PP8OMk2yfbc3Fz9KkWarDKfTt/rADQXkAygTgOYAbC+8ngdgNnOIJKvA/BeADvM7Ie9VmRmE2bWMrPW6OjoIPWKNM/Bg2FiYnHlTlmXeNVpAPcD2EjyApJnALgGwP5qAMlLAfwViv/8j4cvU6TBKlfTdh756RUzdLoSeNnyNgAzew7A9QDuAnAYwB1mdojkzSR3lGEfAPBCAJ8h+QDJ/X1WJyInq9XqGmLnsaAeMUPjyp2yLvFaXSfIzA4AONAxdlNl+XWB6xKRBXU+RWsPQAagK4FFcjcyEiYmFlfulHWJlxqASO5mu865GCwmFlfulHWJlxqASO50FpBEogYgkruJiTAxsbhyp6xLvNQARHJ3551hYmJx5U5Zl3ipAYjkbvfuMDGxuHKnrEu81ABEcjc2FiYmFlfulHWJlxqASO4qn6L7XgmsPQAZgBqASO56zADXNdJ3lrghcOVOWZd4qQGI5K7Hx/6ukb67BkPgyp2yLvFSAxDJ3b59i4t9P1BXYobOlTtlXeKlBiCSu8nJMDGxuHKnrEu81ABEclf5FN33iIr2AGQAagAiudu+vWuo60hQj5ihceVOWZd4qQGI5E7TQUskagAiudu6NUxMLK7cKesSLzUAkdytXRsmJhZX7pR1iZcagEju5ufDxMTiyp2yLvFSAxDJnaaDlkjUAERy126HiYmlkvuLh7+Lx578157PSX5q3RReRBJaRnsA136i+A//0f/2xq7nJD/aAxDJnc4CkkhqNQCS20geITlN8sYez59J8tPl8/eR3BC6UJHG2rt3cbHv1GqVmKFz5U5Zl3h5GwDJVQBuA3AlgM0AdpHc3BF2LYCnzeynAfwZgPeHLlSksTZtWlz80YkTAIBVp7FvzNC5cqesS7zqfAdwGYBpMzsGACRvB7ATwMOVmJ0A9pbLnwXwP0jSLPxcsA88/n184suP9nyuXzpXEc6ZbPv+Tv9fcr7gPk+a47cGmWl34PUN8DuuV9y/PsfaMvg7dAm+zWus70Pv3oHfumU/AODeY08BAD4zNYPnTizF7H3HL2DvX97dP9EAzAwnDDhhBit/nijHzAw/OlEs//kNb8S7/mQSJyov5vc+/QAI4H3v+AX8UeC6VoIrLjwX2y9Jf40EfW8EklcD2GZmbysf/yqAV5nZ9ZWYh8qYmfLxI2XM9zrWNQ5gHADOP//8LY899thJF3zPN4/jffsPOertM+5YJx03rej7jOseGAPkcv+OK1ef9Q14H47g9fX9+3Bs8wH+Dl1F9HvGvV1dqYbzd9hr+MQJQ/uxpzE6cibOOn2VI+OpMRhOI3EaCRLlMsrHS8unETjtNC4u3//o0wCA9S8+K1ptK8GbX/VTuO41/2mg3yU5ZWatEHXU2QPo9a+zs2vUiYGZTQCYAIBWqzXQR7ArLjwXV1x47iC/KrI83XorsGfPqcfE4sqdsi7xqvMl8AyA9ZXH6wDM9oshuRrA2QCeClGgSOPNdr7dBoyJxZU7ZV3iVecQ0GoARwG8FsD/BXA/gN1mdqgS804AF5nZdSSvAfCLZvZLrvW2Wi1r6yIREZGTEvIQkHcPwMyeA3A9gLsAHAZwh5kdInkzyR1l2EcBvITkNIAbAHSdKioiA9qyJUxMLK7cKesSr1pXApvZAQAHOsZuqiz/O4D/GrY0EQGwrK4EPqnnJDldCSySu5GRMDGxuHKnrEu81ABEcjc2FiYmFlfulHWJlxqASO6OHg0TE4srd8q6xEsNQCR3debT0VxAMgDvaaDREpNzAP4VwPd8sQ2xBtoWC7QtCtoOS7QtCmsA/LiZjYZYWbIGAAAk26HOZ13utC2WaFsUtB2WaFsUQm8HHQISEWkoNQARkYZK3QB0lcgSbYsl2hYFbYcl2haFoNsh6XcAIiKSTuo9ABERSUQNQESkoZI1AN+N5pc7kh8jeby8W9rC2ItJ3k3yW+XPF5XjJPnBcls8SPKVld95axn/LZJvTfFaThXJ9STvIXmY5CGSv1OON2p7kHwBya+S/Hq5Hf6oHL+A5H3la/o0yTPK8TPLx9Pl8xsq63pPOX6E5BvSvKJTR3IVya+RnCwfN3JbkHyU5DdIPkCyXY7Ff3+Y2dD/AFgF4BEALwNwBoCvA9icopaIr/G/AHglgIcqY7cAuLFcvhHA+8vlqwB8AcWd1S4HcF85/mIAx8qfLyqXX5T6tQ2wLV4K4JXl8giK+0tsbtr2KF/PC8vl0wHcV76+OwBcU45/GMA7yuXfAvDhcvkaAJ8ulzeX75kzAVxQvpdWpX59A26TGwDsAzBZPm7ktgDwKIA1HWPR3x+p9gAWbzRvZs8CWLjR/IphZl9C913RdgL4RLn8CQBvqox/0gr3AjiH5EsBvAHA3Wb2lJk9DeBuANviVx+WmT1hZv9ULs+juK/EeWjY9ihfz7+UD08v/xiAnwfw2XK8czssbJ/PAngtixsS7wRwu5n90My+DWAaxXtqWSG5DsAbAXykfEw0dFv0Ef39kaoBnAfg8crjmXJspftJM3sCKP5TBLBwc+N+22PFbady1/1SFJ9+G7c9ykMeDwA4juIN+giA71tx4yXg+a9p8fWWzz8D4CVYAduh9OcA3g3gRPn4JWjutjAA/5vkFMnxciz6+6PWDWEiqHUT+Qbptz1W1HYi+UIAfwvgd83sn4sPcL1De4ytiO1hZj8C8HMkzwHwvwC8oldY+XPFbgeSYwCOm9kUya0Lwz1CV/y2KL3azGZJngvgbpLfdMQG2xap9gDq3Gh+JfpuuauG8ufxcrzf9lgx24nk6Sj+8/+fZva5crix28PMvg/gIIpjuOewuPc28PzXtPh6y+fPRnFYcSVsh1cD2EHyURSHgH8exR5BE7cFzGy2/HkcxQeDyzCE90eqBnA/gI3lN/5noPhSZ3+iWoZpP4CFb+bfCuDzlfG3lN/uXw7gmXKX7y4Aryf5ovIMgNeXY8tKeaz2owAOm9mfVp5q1PYgOVp+8gfJswC8DsX3IfcAuLoM69wOC9vnagB/b8W3ffsBXFOeGXMBgI0AvjqcVxGGmb3HzNaZ2QYU7/+/N7M3o4HbguSPkxxZWEbx7/ohDOP9kfBb76tQnA3yCID3pqoj4uv7FIAnAPwHis58LYpjll8E8K3y54vLWAK4rdwW3wDQqqznN1B8sTUN4NdTv64Bt8V/RrEr+iCAB8o/VzVtewC4GMDXyu3wEICbyvGXofhPaxrAZwCcWY6/oHw8XT7/ssq63ltunyMArkz92k5xu2zF0llAjdsW5Wv+evnn0ML/h8N4f2gqCBGRhtKVwCIiDaUGICLSUGoAIiINpQYgItJQagAiIg2lBiAi0lBqALLikHwTyb8m+XmSr29KbpGTpQYgK46Z/Z2ZvR3ArwH4ZaCYhI7kD8qJ2FCO/SZJI/maytj15djr+q2f5E+T/EbH2Jkkvw3gaI/cZ5XzvD9Lck3AlypyStQAZCX7QxRXTC54xMx+rvL4YhRX5b4CAEj+GIortudQXGHZzzEA60lW3z/jAP6PmT3cmdvMflDmXXZz1MjKpgYgyxrJjSQPkmyTvKW8SxJJvh/AF6y8D0EfF6GYsuPC8vFvo5hu4ISZfbdc/wXl4Zw2i7t5vdzMTgD4DoANZcxZAPYA2HsSuUWSUwOQZYvkKgCfBHCDmbUAnIViLpV3oZho7WqS1zlW8QoUd6C6kOTZKA7ZfBnFPD0LM5h+pLL+vSjuzAQUk7gtNI53AthvZo+eRG6R5FLdD0AkhDcBeLjySfswihuKfBDAB12/SHI9gCfN7Fg5B/u7Afx3AJtQHBZaWP/PAPjb8t4FqwH8QyXXy0l+CUUDuBwA6uQWyYUagCxnl6KYWXTBJSjuslXHxVg6zj+P4tZ5l6GYk36hoVyCYmbGj/b4/cMo5rD/HRT3OPjuyZUukp4OAcly9iTKwzAkXwXgLVj69O5zEZYawAcAXG/F3bouqqzjCQBvWPiyl+RFXLqN2WEUDeM3yt8XWXbUAGQ5+xsArfKUzF9E0RCma/7uRSiP9ZvZpJl9pRzfDGDhTJ6PoXiPHC5PH/0DW5o//Ui5jgkze+aUX4lIAjoEJMuWmX0PwKuAxWP6W8szdOr87pv7jJ9bWf4Blu5O1Rn3Q+j9I8uc9gBkpbgE7sM/PwJwdvVCsGFZuBAMwOkAajUokWHQHcFERBpKewAiIg2lBiAi0lBqACIiDaUGICLSUGoAIiINpQYgItJQagAiIg2lBiAi0lBqACIiDfX/AcGzTh/TaOqtAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.clf()\n",
"# plt.plot(x_part, calcs, '.')\n",
"plt.plot(test_q, calcs_test)#, label = 'pdf (Ctt = 0.0)')\n",
"# plt.plot(test_q, calcs_test1, label = 'pdf (Ctt = 0.5)')\n",
"# plt.plot(test_q, calcs_test2, label = 'pdf (D-contribs = 0.3)')\n",
"# plt.plot(test_q, f0_y, label = '0')\n",
"# plt.plot(test_q, fT_y, label = 'T')\n",
"# plt.plot(test_q, fplus_y, label = '+')\n",
"# plt.plot(test_q, res_y, label = 'res')\n",
"plt.axvline(x=jpsi_mass -70,color='red', linewidth=1.0, linestyle = 'dotted')\n",
"plt.axvline(x=jpsi_mass +70,color='red', linewidth=1.0, linestyle = 'dotted')\n",
"plt.axvline(x=psi2s_mass -50,color='red', linewidth=1.0, linestyle = 'dotted')\n",
"plt.axvline(x=psi2s_mass +50,color='red', linewidth=1.0, linestyle = 'dotted')\n",
"# plt.legend()\n",
"# plt.ylim(0.0, 1.5e-6)\n",
"plt.xlabel(r'$q^2 [MeV^2]$')\n",
"# plt.yscale('log')\n",
"# plt.xlim(770, 785)\n",
"plt.savefig('test.png')\n",
"# print(jpsi_width)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"\n",
"\n",
"# probs = mixture.prob(test_q)\n",
"# probs_np = zfit.run(probs)\n",
"# probs_np *= np.max(calcs_test) / np.max(probs_np)\n",
"# plt.figure()\n",
"# plt.semilogy(test_q, probs_np,label=\"importance sampling\")\n",
"# plt.semilogy(test_q, calcs_test, label = 'pdf')\n"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"# 0.213/(0.00133+0.213+0.015)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Adjust scaling of different parts"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"total_f.update_integration_options(draws_per_dim=2000000, mc_sampler=None)\n",
"# inte = total_f.integrate(limits = (950., 1050.), norm_range=False)\n",
"# inte_fl = zfit.run(inte)\n",
"# print(inte_fl/4500)\n",
"# print(pdg[\"jpsi_BR\"]/pdg[\"NR_BR\"], inte_fl*pdg[\"psi2s_auc\"]/pdg[\"NR_auc\"])"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"# # print(\"jpsi:\", inte_fl)\n",
"# # print(\"Increase am by factor:\", np.sqrt(pdg[\"jpsi_BR\"]/pdg[\"NR_BR\"]*pdg[\"NR_auc\"]/inte_fl))\n",
"# # print(\"New amp:\", pdg[\"jpsi\"][3]*np.sqrt(pdg[\"jpsi_BR\"]/pdg[\"NR_BR\"]*pdg[\"NR_auc\"]/inte_fl))\n",
"\n",
"# # print(\"psi2s:\", inte_fl)\n",
"# # print(\"Increase am by factor:\", np.sqrt(pdg[\"psi2s_BR\"]/pdg[\"NR_BR\"]*pdg[\"NR_auc\"]/inte_fl))\n",
"# # print(\"New amp:\", pdg[\"psi2s\"][3]*np.sqrt(pdg[\"psi2s_BR\"]/pdg[\"NR_BR\"]*pdg[\"NR_auc\"]/inte_fl))\n",
"\n",
"# name = \"phi\"\n",
"\n",
"# print(name+\":\", inte_fl)\n",
"# print(\"Increase am by factor:\", np.sqrt(pdg[name+\"_BR\"]/pdg[\"NR_BR\"]*pdg[\"NR_auc\"]/inte_fl))\n",
"# print(\"New amp:\", pdg[name][0]*np.sqrt(pdg[name+\"_BR\"]/pdg[\"NR_BR\"]*pdg[\"NR_auc\"]/inte_fl))\n",
"\n",
"\n",
"# print(x_min)\n",
"# print(x_max)\n",
"# # total_f.update_integration_options(draws_per_dim=2000000, mc_sampler=None)\n",
"# total_f.update_integration_options(mc_sampler=lambda dim, num_results,\n",
"# dtype: tf.random_uniform(maxval=1., shape=(num_results, dim), dtype=dtype),\n",
"# draws_per_dim=1000000)\n",
"# # _ = []\n",
"\n",
"# # for i in range(10):\n",
"\n",
"# # inte = total_f.integrate(limits = (x_min, x_max))\n",
"# # inte_fl = zfit.run(inte)\n",
"# # print(inte_fl)\n",
"# # _.append(inte_fl)\n",
"\n",
"# # print(\"mean:\", np.mean(_))\n",
"\n",
"# _ = time.time()\n",
"\n",
"# inte = total_f.integrate(limits = (x_min, x_max))\n",
"# inte_fl = zfit.run(inte)\n",
"# print(inte_fl)\n",
"# print(\"Time taken: {}\".format(display_time(int(time.time() - _))))\n",
"\n",
"# print(pdg['NR_BR']/pdg['NR_auc']*inte_fl)\n",
"# print(0.25**2*4.2/1000)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Sampling\n",
"## Mixture distribution for sampling"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"\n",
" \n",
"# print(list_of_borders[:9])\n",
"# print(list_of_borders[-9:])\n",
"\n",
"\n",
"class UniformSampleAndWeights(zfit.util.execution.SessionHolderMixin):\n",
" def __call__(self, limits, dtype, n_to_produce):\n",
" # n_to_produce = tf.cast(n_to_produce, dtype=tf.int32)\n",
" low, high = limits.limit1d\n",
" low = tf.cast(low, dtype=dtype)\n",
" high = tf.cast(high, dtype=dtype)\n",
"# uniform = tfd.Uniform(low=low, high=high)\n",
"# uniformjpsi = tfd.Uniform(low=tf.constant(3080, dtype=dtype), high=tf.constant(3112, dtype=dtype))\n",
"# uniformpsi2s = tfd.Uniform(low=tf.constant(3670, dtype=dtype), high=tf.constant(3702, dtype=dtype))\n",
"\n",
"# list_of_borders = []\n",
"# _p = []\n",
"# splits = 10\n",
"\n",
"# _ = np.linspace(x_min, x_max, splits)\n",
"\n",
"# for i in range(splits):\n",
"# list_of_borders.append(tf.constant(_[i], dtype=dtype))\n",
"# _p.append(tf.constant(1/splits, dtype=dtype))\n",
" \n",
"# mixture = tfd.MixtureSameFamily(mixture_distribution=tfd.Categorical(probs=_p[:(splits-1)]),\n",
"# components_distribution=tfd.Uniform(low=list_of_borders[:(splits-1)], \n",
"# high=list_of_borders[-(splits-1):]))\n",
" mixture = tfd.MixtureSameFamily(mixture_distribution=tfd.Categorical(probs=[tf.constant(0.05, dtype=dtype),\n",
" tf.constant(0.93, dtype=dtype),\n",
" tf.constant(0.05, dtype=dtype),\n",
" tf.constant(0.065, dtype=dtype),\n",
" tf.constant(0.04, dtype=dtype),\n",
" tf.constant(0.05, dtype=dtype)]),\n",
" components_distribution=tfd.Uniform(low=[tf.constant(x_min, dtype=dtype), \n",
" tf.constant(3090, dtype=dtype),\n",
" tf.constant(3681, dtype=dtype), \n",
" tf.constant(3070, dtype=dtype),\n",
" tf.constant(1000, dtype=dtype),\n",
" tf.constant(3660, dtype=dtype)], \n",
" high=[tf.constant(x_max, dtype=dtype),\n",
" tf.constant(3102, dtype=dtype), \n",
" tf.constant(3691, dtype=dtype),\n",
" tf.constant(3110, dtype=dtype),\n",
" tf.constant(1040, dtype=dtype),\n",
" tf.constant(3710, dtype=dtype)]))\n",
"# dtype = tf.float64\n",
"# mixture = tfd.MixtureSameFamily(mixture_distribution=tfd.Categorical(probs=[tf.constant(0.04, dtype=dtype),\n",
"# tf.constant(0.90, dtype=dtype),\n",
"# tf.constant(0.02, dtype=dtype),\n",
"# tf.constant(0.07, dtype=dtype),\n",
"# tf.constant(0.02, dtype=dtype)]),\n",
"# components_distribution=tfd.Uniform(low=[tf.constant(x_min, dtype=dtype), \n",
"# tf.constant(3089, dtype=dtype),\n",
"# tf.constant(3103, dtype=dtype), \n",
"# tf.constant(3681, dtype=dtype),\n",
"# tf.constant(3691, dtype=dtype)], \n",
"# high=[tf.constant(3089, dtype=dtype),\n",
"# tf.constant(3103, dtype=dtype), \n",
"# tf.constant(3681, dtype=dtype),\n",
"# tf.constant(3691, dtype=dtype), \n",
"# tf.constant(x_max, dtype=dtype)]))\n",
"# mixture = tfd.Uniform(tf.constant(x_min, dtype=dtype), tf.constant(x_max, dtype=dtype))\n",
"# sample = tf.random.uniform((n_to_produce, 1), dtype=dtype)\n",
" sample = mixture.sample((n_to_produce, 1))\n",
"# sample = tf.random.uniform((n_to_produce, 1), dtype=dtype)\n",
" weights = mixture.prob(sample)[:,0]\n",
"# weights = tf.broadcast_to(tf.constant(1., dtype=dtype), shape=(n_to_produce,))\n",
" # sample = tf.expand_dims(sample, axis=-1)\n",
"# print(sample, weights)\n",
" \n",
"# weights = tf.ones(shape=(n_to_produce,), dtype=dtype)\n",
" weights_max = None\n",
" thresholds = tf.random_uniform(shape=(n_to_produce,), dtype=dtype)\n",
" return sample, thresholds, weights, weights_max, n_to_produce"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"# total_f._sample_and_weights = UniformSampleAndWeights"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Constraints"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [],
"source": [
"# 1. Constraint - Real part of sum of Psi contrib and D contribs\n",
"\n",
"sum_list = []\n",
"\n",
"sum_list.append(ztf.to_complex(jpsi_s) * tf.exp(tf.complex(ztf.constant(0.0), jpsi_p)) * ztf.to_complex(jpsi_w / (tf.pow(jpsi_m,3))))\n",
"sum_list.append(ztf.to_complex(psi2s_s) * tf.exp(tf.complex(ztf.constant(0.0), psi2s_p)) * ztf.to_complex(psi2s_w / (tf.pow(psi2s_m,3))))\n",
"sum_list.append(ztf.to_complex(p3770_s) * tf.exp(tf.complex(ztf.constant(0.0), p3770_p)) * ztf.to_complex(p3770_w / (tf.pow(p3770_m,3))))\n",
"sum_list.append(ztf.to_complex(p4040_s) * tf.exp(tf.complex(ztf.constant(0.0), p4040_p)) * ztf.to_complex(p4040_w / (tf.pow(p4040_m,3))))\n",
"sum_list.append(ztf.to_complex(p4160_s) * tf.exp(tf.complex(ztf.constant(0.0), p4160_p)) * ztf.to_complex(p4160_w / (tf.pow(p4160_m,3))))\n",
"sum_list.append(ztf.to_complex(p4415_s) * tf.exp(tf.complex(ztf.constant(0.0), p4415_p)) * ztf.to_complex(p4415_w / (tf.pow(p4415_m,3))))\n",
"sum_list.append(ztf.to_complex(DDstar_s) * tf.exp(tf.complex(ztf.constant(0.0), DDstar_p)) * (ztf.to_complex(1.0 / (10.0*tf.pow(Dstar_m,2)) + 1.0 / (10.0*tf.pow(D_m,2)))))\n",
"sum_list.append(ztf.to_complex(Dbar_s) * tf.exp(tf.complex(ztf.constant(0.0), Dbar_p)) * ztf.to_complex(1.0 / (6.0*tf.pow(Dbar_m,2))))\n",
"\n",
"sum_ru_1 = ztf.to_complex(ztf.constant(0.0))\n",
"\n",
"for part in sum_list:\n",
" sum_ru_1 += part\n",
"\n",
"sum_1 = tf.math.real(sum_ru_1)\n",
"# constraint1 = zfit.constraint.GaussianConstraint(params = sum_1, mu = ztf.constant(1.7*10**-8), \n",
"# sigma = ztf.constant(2.2*10**-8))\n",
"\n",
"constraint1 = tf.pow((sum_1-ztf.constant(1.7*10**-8))/ztf.constant(2.2*10**-8),2)/ztf.constant(2.)\n",
"\n",
"# 2. Constraint - Abs. of sum of Psi contribs and D contribs\n",
"\n",
"sum_2 = tf.abs(sum_ru_1)\n",
"constraint2 = tf.cond(tf.greater_equal(sum_2, 5.0e-8), lambda: 100000., lambda: 0.)\n",
"\n",
"# 3. Constraint - Maximum eta of D contribs\n",
"\n",
"constraint3_0 = tf.cond(tf.greater_equal(tf.abs(Dbar_s), 0.2), lambda: 100000., lambda: 0.)\n",
"\n",
"constraint3_1 = tf.cond(tf.greater_equal(tf.abs(DDstar_s), 0.2), lambda: 100000., lambda: 0.)\n",
"\n",
"# 4. Constraint - Formfactor multivariant gaussian covariance fplus\n",
"\n",
"Cov_matrix = [[ztf.constant( 1.), ztf.constant( 0.45), ztf.constant( 0.19), ztf.constant(0.857), ztf.constant(0.598), ztf.constant(0.531), ztf.constant(0.752), ztf.constant(0.229), ztf.constant(0,117)],\n",
" [ztf.constant( 0.45), ztf.constant( 1.), ztf.constant(0.677), ztf.constant(0.708), ztf.constant(0.958), ztf.constant(0.927), ztf.constant(0.227), ztf.constant(0.443), ztf.constant(0.287)],\n",
" [ztf.constant( 0.19), ztf.constant(0.677), ztf.constant( 1.), ztf.constant(0.595), ztf.constant(0.770), ztf.constant(0.819),ztf.constant(-0.023), ztf.constant( 0.07), ztf.constant(0.196)],\n",
" [ztf.constant(0.857), ztf.constant(0.708), ztf.constant(0.595), ztf.constant( 1.), ztf.constant( 0.83), ztf.constant(0.766), ztf.constant(0.582), ztf.constant(0.237), ztf.constant(0.192)],\n",
" [ztf.constant(0.598), ztf.constant(0.958), ztf.constant(0.770), ztf.constant( 0.83), ztf.constant( 1.), ztf.constant(0.973), ztf.constant(0.324), ztf.constant(0.372), ztf.constant(0.272)],\n",
" [ztf.constant(0.531), ztf.constant(0.927), ztf.constant(0.819), ztf.constant(0.766), ztf.constant(0.973), ztf.constant( 1.), ztf.constant(0.268), ztf.constant(0.332), ztf.constant(0.269)],\n",
" [ztf.constant(0.752), ztf.constant(0.227),ztf.constant(-0.023), ztf.constant(0.582), ztf.constant(0.324), ztf.constant(0.268), ztf.constant( 1.), ztf.constant( 0.59), ztf.constant(0.515)],\n",
" [ztf.constant(0.229), ztf.constant(0.443), ztf.constant( 0.07), ztf.constant(0.237), ztf.constant(0.372), ztf.constant(0.332), ztf.constant( 0.59), ztf.constant( 1.), ztf.constant(0.897)],\n",
" [ztf.constant(0.117), ztf.constant(0.287), ztf.constant(0.196), ztf.constant(0.192), ztf.constant(0.272), ztf.constant(0.269), ztf.constant(0.515), ztf.constant(0.897), ztf.constant( 1.)]]\n",
"\n",
"def triGauss(val1,val2,val3,m = Cov_matrix):\n",
"\n",
" mean1 = ztf.constant(0.466)\n",
" mean2 = ztf.constant(-0.885)\n",
" mean3 = ztf.constant(-0.213)\n",
" sigma1 = ztf.constant(0.014/3.)\n",
" sigma2 = ztf.constant(0.128/3.)\n",
" sigma3 = ztf.constant(0.548/3.)\n",
" x1 = (val1-mean1)/sigma1\n",
" x2 = (val2-mean2)/sigma2\n",
" x3 = (val3-mean3)/sigma3\n",
" rho12 = m[0][1]\n",
" rho13 = m[0][2]\n",
" rho23 = m[1][2]\n",
" w = x1*x1*(rho23*rho23-1) + x2*x2*(rho13*rho13-1)+x3*x3*(rho12*rho12-1)+2*(x1*x2*(rho12-rho13*rho23)+x1*x3*(rho13-rho12*rho23)+x2*x3*(rho23-rho12*rho13))\n",
" d = 2*(rho12*rho12+rho13*rho13+rho23*rho23-2*rho12*rho13*rho23-1)\n",
" \n",
" fcn = -w/d\n",
" chisq = -2*fcn\n",
" return chisq\n",
"\n",
"constraint4 = triGauss(bplus_0, bplus_1, bplus_2)\n",
"\n",
"# mean1 = ztf.constant(0.466)\n",
"# mean2 = ztf.constant(-0.885)\n",
"# mean3 = ztf.constant(-0.213)\n",
"# sigma1 = ztf.constant(0.014)\n",
"# sigma2 = ztf.constant(0.128)\n",
"# sigma3 = ztf.constant(0.548)\n",
"# constraint4_0 = tf.pow((bplus_0-mean1)/sigma1,2)/ztf.constant(2.)\n",
"# constraint4_1 = tf.pow((bplus_1-mean2)/sigma2,2)/ztf.constant(2.)\n",
"# constraint4_2 = tf.pow((bplus_2-mean3)/sigma3,2)/ztf.constant(2.)\n",
"\n",
"# 5. Constraint - Abs. of sum of light contribs\n",
"\n",
"sum_list_5 = []\n",
"\n",
"sum_list_5.append(rho_s*rho_w/rho_m)\n",
"sum_list_5.append(omega_s*omega_w/omega_m)\n",
"sum_list_5.append(phi_s*phi_w/phi_m)\n",
"\n",
"\n",
"sum_ru_5 = ztf.constant(0.0)\n",
"\n",
"for part in sum_list_5:\n",
" sum_ru_5 += part\n",
"\n",
"constraint5 = tf.cond(tf.greater_equal(tf.abs(sum_ru_5), ztf.constant(0.02)), lambda: 100000., lambda: 0.)\n",
"\n",
"# 6. Constraint on phases of Jpsi and Psi2s for cut out fit\n",
"\n",
"\n",
"# constraint6_0 = zfit.constraint.GaussianConstraint(params = jpsi_p, mu = ztf.constant(pdg[\"jpsi_phase_unc\"]),\n",
"# sigma = ztf.constant(jpsi_phase))\n",
"# constraint6_1 = zfit.constraint.GaussianConstraint(params = psi2s_p, mu = ztf.constant(pdg[\"psi2s_phase_unc\"]),\n",
"# sigma = ztf.constant(psi2s_phase))\n",
"\n",
"constraint6_0 = tf.pow((jpsi_p-ztf.constant(jpsi_phase))/ztf.constant(pdg[\"jpsi_phase_unc\"]),2)/ztf.constant(2.)\n",
"constraint6_1 = tf.pow((psi2s_p-ztf.constant(psi2s_phase))/ztf.constant(pdg[\"psi2s_phase_unc\"]),2)/ztf.constant(2.)\n",
"\n",
"# 7. Constraint on Ctt with higher limits\n",
"\n",
"constraint7 = tf.cond(tf.greater_equal(Ctt*Ctt, 0.25), lambda: 100000., lambda: 0.)\n",
"\n",
"constraint7dtype = tf.float64\n",
"\n",
"# zfit.run(constraint6_0)\n",
"\n",
"# ztf.convert_to_tensor(constraint6_0)\n",
"\n",
"#List of all constraints\n",
"\n",
"constraints = [constraint1, constraint2, constraint3_0, constraint3_1,# constraint4, #constraint4_0, constraint4_1, constraint4_2,\n",
" constraint6_0, constraint6_1]#, constraint7]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Reset params"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"param_values_dic = {\n",
" 'jpsi_m': jpsi_mass,\n",
" 'jpsi_s': jpsi_scale,\n",
" 'jpsi_p': jpsi_phase,\n",
" 'jpsi_w': jpsi_width,\n",
" 'psi2s_m': psi2s_mass,\n",
" 'psi2s_s': psi2s_scale,\n",
" 'psi2s_p': psi2s_phase,\n",
" 'psi2s_w': psi2s_width,\n",
" 'p3770_m': p3770_mass,\n",
" 'p3770_s': p3770_scale,\n",
" 'p3770_p': p3770_phase,\n",
" 'p3770_w': p3770_width,\n",
" 'p4040_m': p4040_mass,\n",
" 'p4040_s': p4040_scale,\n",
" 'p4040_p': p4040_phase,\n",
" 'p4040_w': p4040_width,\n",
" 'p4160_m': p4160_mass,\n",
" 'p4160_s': p4160_scale,\n",
" 'p4160_p': p4160_phase,\n",
" 'p4160_w': p4160_width,\n",
" 'p4415_m': p4415_mass,\n",
" 'p4415_s': p4415_scale,\n",
" 'p4415_p': p4415_phase,\n",
" 'p4415_w': p4415_width,\n",
" 'rho_m': rho_mass,\n",
" 'rho_s': rho_scale,\n",
" 'rho_p': rho_phase,\n",
" 'rho_w': rho_width,\n",
" 'omega_m': omega_mass,\n",
" 'omega_s': omega_scale,\n",
" 'omega_p': omega_phase,\n",
" 'omega_w': omega_width,\n",
" 'phi_m': phi_mass,\n",
" 'phi_s': phi_scale,\n",
" 'phi_p': phi_phase,\n",
" 'phi_w': phi_width,\n",
" 'Dstar_m': Dstar_mass,\n",
" 'DDstar_s': 0.0,\n",
" 'DDstar_p': 0.0,\n",
" 'D_m': D_mass,\n",
" 'Dbar_m': Dbar_mass,\n",
" 'Dbar_s': 0.0,\n",
" 'Dbar_p': 0.0,\n",
" 'tau_m': pdg['tau_M'],\n",
" 'Ctt': 0.0,\n",
" 'b0_0': 0.292,\n",
" 'b0_1': 0.281,\n",
" 'b0_2': 0.150,\n",
" 'bplus_0': 0.466,\n",
" 'bplus_1': -0.885,\n",
" 'bplus_2': -0.213,\n",
" 'bT_0': 0.460,\n",
" 'bT_1': -1.089,\n",
" 'bT_2': -1.114}\n",
"\n",
"\n",
"def reset_param_values(variation = 0.05):\n",
" for param in total_f_fit.get_dependents():\n",
" if param.floating:\n",
" param.set_value(param_values_dic[param.name] + random.uniform(-1, 1) * param_values_dic[param.name]* variation)\n",
"# print(param.name)\n",
"# for param in totalf.get_dependents():\n",
"# param.set_value()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Pull dictionary"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"# BR_steps = np.linspace(0.0, 1e-3, 11)\n",
"pull_dic = {}\n",
"\n",
"mi = 2e-4\n",
"ma = 6e-4\n",
"ste = 5\n",
"\n",
"# mi = 1e-4\n",
"# ma = 3e-3\n",
"# ste = 20\n",
"\n",
"for param in total_f_fit.get_dependents():\n",
" if param.floating:\n",
" pull_dic[param.name] = []\n",
" for step in range(2*ste):\n",
" pull_dic[param.name].append([])\n",
" \n",
"\n",
"\n",
"def save_pulls(step):\n",
" for param in total_f_fit.get_dependents():\n",
" if param.floating:\n",
" pull_dic[param.name][step].append((params[param]['value'] - param_values_dic[param.name])/params[param]['minuit_hesse']['error'])\n",
"\n",
"\n",
"\n",
"# for key in pull_dic.keys():\n",
"# print(np.shape(pull_dic[key]))\n",
"# save_pulls(New_step=True)\n",
"# params[Ctt]['value']"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"36668\n",
"5404696\n"
]
}
],
"source": [
"# for param in total_f_fit.get_dependents():\n",
"# if param.floating:\n",
"# print(param.name)\n",
"\n",
"# print(params[Ctt])\n",
"\n",
"total_BR = 1.7e-10 + 4.9e-10 + 2.5e-9 + 6.02e-5 + 4.97e-6 + 1.38e-9 + 4.2e-10 + 2.6e-9 + 6.1e-10 + 4.37e-7\n",
"cut_BR = 1.0 - (6.02e-5 + 4.97e-6)/total_BR\n",
"\n",
"nevents = int(pdg[\"number_of_decays\"]*cut_BR)\n",
"\n",
"print(nevents)\n",
"print(int(pdg[\"number_of_decays\"]))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# CLS Code"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[0. 0.0003 0.0004 0.0005 0.0006]\n",
"[0. 0.26726124 0.3086067 0.34503278 0.37796447]\n",
"WARNING:tensorflow:From C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\core\\sample.py:163: to_int32 (from tensorflow.python.ops.math_ops) is deprecated and will be removed in a future version.\n",
"Instructions for updating:\n",
"Use tf.cast instead.\n"
]
}
],
"source": [
"# zfit.run.numeric_checks = False \n",
"\n",
"load = True\n",
"\n",
"bo = True\n",
"\n",
"D_contribs = True\n",
"\n",
"if not D_contribs:\n",
" Dbar_s.floating = False\n",
" Dbar_p.floating = False\n",
" DDstar_s.floating = False\n",
" DDstar_p.floating = False\n",
"\n",
"bo_set = 1\n",
"\n",
"fitting_range = 'cut'\n",
"total_BR = 1.7e-10 + 4.9e-10 + 2.5e-9 + 6.02e-5 + 4.97e-6 + 1.38e-9 + 4.2e-10 + 2.6e-9 + 6.1e-10 + 4.37e-7\n",
"cut_BR = 1.0 - (6.02e-5 + 4.97e-6)/total_BR\n",
"\n",
"Ctt_list = []\n",
"Ctt_error_list = []\n",
"\n",
"nr_of_toys = 1\n",
"nevents = int(pdg[\"number_of_decays\"]*cut_BR)\n",
"# nevents = pdg[\"number_of_decays\"]\n",
"event_stack = 1000000\n",
"# nevents *= 41\n",
"# zfit.settings.set_verbosity(10)\n",
"\n",
"# mi = 1e-4\n",
"# ma = 3e-3\n",
"# ste = 13\n",
"\n",
"BR_steps = np.linspace(mi, ma, ste)\n",
"\n",
"BR_steps[0] = 0.0\n",
"\n",
"print(BR_steps)\n",
"\n",
"Ctt_steps = np.sqrt(BR_steps/4.2*1000)\n",
"\n",
"print(Ctt_steps)\n",
"\n",
"# total_samp = []\n",
"\n",
"start = time.time()\n",
"\n",
"Nll_list = []\n",
"\n",
"sampler = total_f.create_sampler(n=nevents, fixed_params = False)\n",
"sampler.set_data_range(obs_fit)\n",
"\n",
"__ = -1\n",
"\n",
"#-----------------------------------------------------\n",
"\n",
"if not load:\n",
" for Ctt_step in Ctt_steps:\n",
" \n",
" __ += 1\n",
" \n",
" for i in range(2):\n",
" Ctt_list.append([])\n",
" Ctt_error_list.append([])\n",
" Nll_list.append([])\n",
"\n",
" for param in total_f_fit.get_dependents():\n",
" if param.floating:\n",
" pull_dic[param.name].append([])\n",
" \n",
" for toy in range(nr_of_toys): \n",
" \n",
" newset = True\n",
" \n",
" while newset:\n",
" \n",
" for floaty in [True, False]:\n",
" Ctt.floating = floaty\n",
" \n",
" for bo_step in range(bo_set):\n",
"\n",
" print('Step: {0}/{1}'.format(int(__), ste))\n",
" print('Current Ctt: {0}'.format(Ctt_step))\n",
" print('Ctt floating: {0}'.format(floaty))\n",
" \n",
" reset_param_values(variation = 0.0)\n",
"\n",
" if floaty:\n",
" print('Toy {0}/{1} - Fit {2}/{3}'.format(toy, nr_of_toys, bo_step, bo_set))\n",
" Ctt.set_value(Ctt_step)\n",
"\n",
" else:\n",
" Ctt.set_value(0.0)\n",
" print('Toy {0}/{1} - Fit {2}/{3}'.format(toy, nr_of_toys, bo_step, bo_set))\n",
"\n",
" if newset:\n",
" sampler.resample(n=nevents)\n",
" data = sampler\n",
" newset = False\n",
"\n",
" ### Fit data\n",
" \n",
" if floaty:\n",
" plt.clf()\n",
" plt.title('Ctt value: {:.2f}'.format(Ctt_step))\n",
" plt.hist(zfit.run(data), bins = int((x_max-x_min)/7), range = (x_min, x_max))\n",
" plt.savefig('data/CLs/plots/set_histo{}.png'.format(__))\n",
" _step = 2*__\n",
" else:\n",
" _step = 2*__+1\n",
"\n",
" nll = zfit.loss.UnbinnedNLL(model=total_f_fit, data=data, constraints = constraints)\n",
"\n",
" minimizer = zfit.minimize.MinuitMinimizer(verbosity = 5)\n",
" # minimizer._use_tfgrad = False\n",
" result = minimizer.minimize(nll)\n",
"\n",
" print(\"Function minimum:\", result.fmin)\n",
" print(\"Hesse errors:\", result.hesse())\n",
"\n",
" params = result.params\n",
"\n",
" if result.converged:\n",
" \n",
" save_pulls(step = _step)\n",
"\n",
" if floaty:\n",
" Nll_list[-2].append(result.fmin)\n",
" Ctt_list[-2].append(params[Ctt]['value'])\n",
" Ctt_error_list[-2].append(params[Ctt]['minuit_hesse']['error'])\n",
"\n",
" else:\n",
" Nll_list[-1].append(result.fmin)\n",
" Ctt_list[-1].append(0.0)\n",
" Ctt_error_list[-1].append(0.0)\n",
" \n",
"\n",
" else:\n",
" for _ in [1,2]:\n",
" del Nll_list[-_][toy*bo_set:]\n",
"# print(np.shape(Nll_list[-_]))\n",
" del Ctt_list[-_][toy*bo_set:]\n",
" del Ctt_error_list[-_][toy*bo_set:]\n",
" for param in total_f_fit.get_dependents():\n",
" if param.floating:\n",
" del pull_dic[param.name][_step+1-_][toy*bo_set:]\n",
" newset = True\n",
" break\n",
" \n",
" if not result.converged:\n",
" break\n",
" \n",
" print()\n",
" print('Time taken: {}'.format(display_time(int(time.time()-start))))\n",
" print('Estimated time left: {}'.format(display_time(int((time.time()-start)/(__+(toy+1)/nr_of_toys)*(ste-__-(nr_of_toys-toy-1)/nr_of_toys)))))"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Number of jobs: 300\n",
"2317442\n",
"2317446\n",
"2317447\n",
"2317455\n",
"2317458\n",
"2317459\n",
"2317461\n",
"2317462\n",
"2317464\n",
"2317465\n",
"2317931\n",
"2317932\n",
"2317933\n",
"2317934\n",
"2317935\n",
"2317936\n",
"2317937\n",
"2317938\n",
"2317939\n",
"2317940\n",
"2317941\n",
"2317942\n",
"2317943\n",
"2317944\n",
"2317945\n",
"2317946\n",
"2317947\n",
"2317948\n",
"2317949\n",
"2317950\n",
"2317951\n",
"2317952\n",
"2317953\n",
"2317954\n",
"2317955\n",
"2317956\n",
"2317957\n",
"2317958\n",
"2317959\n",
"2317960\n",
"2317961\n",
"2317962\n",
"2317963\n",
"2317964\n",
"2317965\n",
"2317966\n",
"2317967\n",
"2317968\n",
"2317969\n",
"2317970\n",
"2317971\n",
"2317972\n",
"2317973\n",
"2317974\n",
"2317975\n",
"2317976\n",
"2317977\n",
"2317978\n",
"2317979\n",
"2317980\n",
"2317981\n",
"2317982\n",
"2317983\n",
"2317984\n",
"2317985\n",
"2317986\n",
"2317987\n",
"2317988\n",
"2317989\n",
"2317990\n",
"2317991\n",
"2317992\n",
"2317993\n",
"2317994\n",
"2317995\n",
"2317996\n",
"2317997\n",
"2317998\n",
"2317999\n",
"2318000\n",
"2318001\n",
"2318002\n",
"2318003\n",
"2318004\n",
"2318005\n",
"2318006\n",
"2318007\n",
"2318008\n",
"2318009\n",
"2318010\n",
"2318011\n",
"2318012\n",
"2318013\n",
"2318014\n",
"2318015\n",
"2318016\n",
"2318017\n",
"2318018\n",
"2318019\n",
"2318020\n",
"2318021\n",
"2318022\n",
"2318023\n",
"2318024\n",
"2318025\n",
"2318026\n",
"2318027\n",
"2318028\n",
"2318029\n",
"2318030\n",
"2318031\n",
"2318032\n",
"2318033\n",
"2318034\n",
"2318035\n",
"2318036\n",
"2318037\n",
"2318038\n",
"2318039\n",
"2318040\n",
"2318041\n",
"2318042\n",
"2318043\n",
"2318044\n",
"2318045\n",
"2318046\n",
"2318047\n",
"2318048\n",
"2318049\n",
"2318050\n",
"2318051\n",
"2318052\n",
"2318053\n",
"2318054\n",
"2318055\n",
"2318056\n",
"2318057\n",
"2318058\n",
"2318059\n",
"2318060\n",
"2318061\n",
"2318062\n",
"2318063\n",
"2318064\n",
"2318065\n",
"2318066\n",
"2318067\n",
"2318068\n",
"2318069\n",
"2318070\n",
"2318071\n",
"2318072\n",
"2318073\n",
"2318074\n",
"2318075\n",
"2318076\n",
"2318077\n",
"2318078\n",
"2318079\n",
"2318080\n",
"2318081\n",
"2318082\n",
"2318083\n",
"2318084\n",
"2318085\n",
"2318086\n",
"2318087\n",
"2318088\n",
"2318089\n",
"2318090\n",
"2318091\n",
"2318092\n",
"2318093\n",
"2318094\n",
"2318095\n",
"2318096\n",
"2318097\n",
"2318098\n",
"2318099\n",
"2318100\n",
"2318101\n",
"2318102\n",
"2318103\n",
"2318104\n",
"2318105\n",
"2318106\n",
"2318107\n",
"2318108\n",
"2318109\n",
"2318110\n",
"2318111\n",
"2318112\n",
"2318113\n",
"2318114\n",
"2318115\n",
"2318116\n",
"2318117\n",
"2318118\n",
"2318119\n",
"2318120\n",
"2318121\n",
"2318122\n",
"2318123\n",
"2318124\n",
"2318125\n",
"2318126\n",
"2318127\n",
"2318128\n",
"2318129\n",
"2318130\n",
"2318131\n",
"2318132\n",
"2318133\n",
"2318134\n",
"2318135\n",
"2318136\n",
"2318137\n",
"2318138\n",
"2318139\n",
"2318140\n",
"2318141\n",
"2318142\n",
"2318143\n",
"2318144\n",
"2318145\n",
"2318146\n"
]
}
],
"source": [
"if load:\n",
" \n",
" phase_combi = '-+'\n",
" \n",
" if D_contribs:\n",
" D_dir = 'D-True'\n",
" else:\n",
" D_dir = 'D-False'\n",
"\n",
" _dir = 'data/CLs/finished/f3d41/{}/{}/'.format(phase_combi, D_dir)\n",
" \n",
" jobs = os.listdir(_dir)\n",
" \n",
" First = True\n",
" \n",
" print('Number of jobs: {}'.format(len(jobs)))\n",
" \n",
" for job in jobs:\n",
" \n",
" dirName = _dir + str(job) + '/data/CLs'\n",
" \n",
" if not os.path.exists(\"{}/{}-{}_{}s{}b{}t--CLs_Nll_list.pkl\".format(dirName, mi,ma,ste,bo_set,nr_of_toys)):\n",
" print(job)\n",
" continue\n",
" \n",
" with open(r\"{}/variab.pkl\".format(dirName), \"rb\") as input_file:\n",
" variab = pkl.load(input_file)\n",
"# print(variab)\n",
" \n",
" ### sanity check:\n",
" if variab['mi'] != mi or variab['ma'] != ma or variab['ste'] != ste or bo_set != bo_set:\n",
" print('Fitting parameters of data dont equal the ones given -- Job {} skipped!'.format(job))\n",
" continue\n",
" \n",
" with open(r\"{}/{}-{}_{}s{}b{}t--CLs_Nll_list.pkl\".format(dirName, mi,ma,ste,bo_set,nr_of_toys), \"rb\") as input_file:\n",
" _Nll_list = pkl.load(input_file)\n",
" \n",
" with open(r\"{}/{}-{}_{}s{}b{}t--Ctt_list.pkl\".format(dirName, mi,ma,ste,bo_set,nr_of_toys), \"rb\") as input_file:\n",
" _Ctt_list = pkl.load(input_file)\n",
" \n",
" with open(r\"{}/{}-{}_{}s{}b{}t--Ctt_error_list.pkl\".format(dirName, mi,ma,ste,bo_set,nr_of_toys), \"rb\") as input_file:\n",
" _Ctt_error_list = pkl.load(input_file)\n",
" \n",
" with open(r\"{}/{}-{}_{}s{}b{}t--pull_dic.pkl\".format(dirName, mi,ma,ste,bo_set,nr_of_toys), \"rb\") as input_file:\n",
" _pull_dic = pkl.load(input_file)\n",
" \n",
" with open(r\"{}/{}-{}_{}s--CLs_list.pkl\".format(dirName, mi,ma,ste), \"rb\") as input_file:\n",
" _CLs_list = pkl.load(input_file)\n",
" \n",
" \n",
" if First:\n",
" Nll_list = _Nll_list\n",
" Ctt_list = _Ctt_list\n",
" Ctt_error_list = _Ctt_error_list\n",
" pull_dic = _pull_dic\n",
"# print(_pull_dic)\n",
" CLs_list = _CLs_list\n",
" First = False\n",
" else:\n",
" for step in range(2*ste):\n",
"# print(Nll_list[step], step)\n",
" Nll_list[step].extend(_Nll_list[step])\n",
" Ctt_list[step].extend(_Ctt_list[step])\n",
" Ctt_error_list[step].extend(_Ctt_error_list[step])\n",
" for key in pull_dic.keys():\n",
"# print(key, np.shape(pull_dic[key]))\n",
" pull_dic[key][step].extend(_pull_dic[key][step])\n",
" for step in range(ste):\n",
" CLs_list[step].extend(_CLs_list[step])\n",
"\n",
"# print('----------------------')"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [],
"source": [
"dirName = 'data/CLs'\n",
"\n",
"# if bo and not load:\n",
"# for s in range(2*ste):\n",
"# Nll_list[s] = [np.min(Nll_list[s])]\n",
"\n",
"\n",
"if not load:\n",
" \n",
" if not os.path.exists(dirName):\n",
" os.mkdir(dirName)\n",
" print(\"Directory \" , dirName , \" Created \")\n",
"\n",
" with open(\"{}/{}-{}_{}s{}b{}t--CLs_Nll_list.pkl\".format(dirName, mi,ma,ste,bo_set,nr_of_toys), \"wb\") as f:\n",
" pkl.dump(Nll_list, f, pkl.HIGHEST_PROTOCOL)\n",
" \n",
" with open(\"{}/{}-{}_{}s{}b{}t--Ctt_list.pkl\".format(dirName, mi,ma,ste,bo_set,nr_of_toys), \"wb\") as f:\n",
" pkl.dump(Ctt_list, f, pkl.HIGHEST_PROTOCOL)\n",
" \n",
" with open(\"{}/{}-{}_{}s{}b{}t--Ctt_error_list.pkl\".format(dirName, mi,ma,ste,bo_set,nr_of_toys), \"wb\") as f:\n",
" pkl.dump(Ctt_error_list, f, pkl.HIGHEST_PROTOCOL)\n",
" \n",
" with open(\"{}/{}-{}_{}s{}b{}t--pull_dic.pkl\".format(dirName, mi,ma,ste,bo_set,nr_of_toys), \"wb\") as f:\n",
" pkl.dump(pull_dic, f, pkl.HIGHEST_PROTOCOL)\n",
" \n",
" variab = {'mi': mi,\n",
" 'ma': ma,\n",
" 'ste': ste,\n",
" 'bo_set': bo_set,\n",
" 'nr_of_toys': nr_of_toys}\n",
" \n",
" with open(\"{}/variab.pkl\".format(dirName), \"wb\") as f:\n",
" pkl.dump(variab, f, pkl.HIGHEST_PROTOCOL)\n",
" \n",
" CLs_values = []\n",
" \n",
" toy_size = bo_set\n",
"\n",
" print(np.shape(Nll_list))\n",
" print(Nll_list[0:1])\n",
" \n",
" for step in range(ste):\n",
" CLs_values.append([])\n",
" for toy in range(nr_of_toys):\n",
" float_min = np.min(Nll_list[2*step][toy*bo_set:(toy+1)*bo_set])\n",
" fix_min = np.min(Nll_list[2*step+1][toy*bo_set:(toy+1)*bo_set])\n",
" CLs_values[step].append(float_min-fix_min)\n",
" \n",
" \n",
" print(np.shape(CLs_values))\n",
" \n",
" with open(\"{}/{}-{}_{}s--CLs_list.pkl\".format(dirName, mi,ma,ste), \"wb\") as f:\n",
" pkl.dump(CLs_values, f, pkl.HIGHEST_PROTOCOL)"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [],
"source": [
"# print(variab['mi'] != mi)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Plot"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(5, 74)\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAANcAAAFhCAYAAADjmT9pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAMTQAADE0B0s6tTgAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deVxUZfs/8M/sC6g8Qfh1AxJFVFDAjdzRVCRJc8slFVGxNM3HJf3WQ+aSXy3jZ26P9k3U7CthaqWmhrlhpj76KC6p5A4UbhSlOPvcvz/IeUQWWeb2HuZc79drXjoz95xzzcA193UO55xLxhhjIIQ4nVx0AIS4K0ouQjih5CKEE0ouQjih5CKEE0ouQjih5CKEE0ouJ5LJZPj+++9Lff69995Dx44dn2JEVVPd4nU5rJrp0qULA8AAMLlczurVq8cmT57MjEaj6NAYALZnz55Sn7937x7Ly8t7ihGVX4cOHdjs2bOLPMYr3o0bN7LQ0FCmVqtZnTp12AcffFDm+KNHj7JWrVoxjUbDmjRpwrZv317k+Ye/D4/fcnNzS12mxWJhf//735m3tzfz9PRkI0eOZPfu3StxbFJSEgPA3nnnnQq9z2o5c02ZMgW5ubnIysrCunXrsHXrVsybN6/U8SaT6SlGVzpPT08888wzT3WdVXnvPOLdsGED3njjDUybNg3nz5/Hjh070Lp161LH5+XloXfv3oiMjMTJkycxatQoDBgwABcvXiwybtOmTcjNzS1y8/X1LXW58+bNw8aNG5Gamoq9e/fixIkTmDBhQrFxFy5cwLJlyxAaGlrxN1uhVHQBXbp0KfYN8uabb7KYmBjH/VGjRrFhw4axmTNnMm9vb9anTx/HuOeee47pdDrWrFkz9sUXXxRZjr+/P1u8eDEbOHAg0+v1LDg4mO3du7fImLS0NNa2bVum0WiYr68ve/311x3PAWBr165l3bt3ZzqdjkVERLDTp087np89ezbr0KFDkfcyY8YMlpCQwDw9PZm/vz9LSUkpsr7/+7//Yw0aNGB6vZ6NHDmSTZs2jXXp0qXUz+fhOpKSkljdunVZSEgIY4yxBQsWsODgYKbT6VijRo3Yxx9/XOTzwiPf+P7+/iXGe//+fTZmzBjm5eXFPDw8WP/+/dnNmzdLjeVxZrOZ1a5dm61bt67cr/n4449ZvXr1mN1udzzWqVMn9uabbzru4wkVw+NsNhvz8fFhn3zyieOxvXv3MoVCwe7evet4zGKxsFatWrFt27aV+Hv3JNVy5npUdnY2vv/+e7Rp06bI49988w0MBgMOHz6MpKQkAIC3tze++OILnDt3DpMmTcKIESNw9uzZIq/74IMPEBsbi4yMDHTq1AmvvvoqzGYzAOD8+fN48cUX8cILL+DUqVPYtWsXmjZtWuT1c+fOxaRJk5CRkYG6deti9OjRZca/evVqBAcH49SpU4iLi8Po0aNx+/ZtAEBmZiZGjhyJ119/HSdPnkRQUBA++eSTJ34mGRkZOHbsGNLS0rBp0yYAgEajwf/+7//ip59+wvvvv4+3334bO3fuBAB8/PHHaNu2LaZNm4bc3FwcP368xOX+/e9/x8GDB/HNN98gPT0dv/zyC0aMGOF4ft26dZDJZKXG9e9//xu3bt2CxWJBSEgIGjRogFGjRiEvL6/U1/zrX/9CVFRUkeV2794dx44de+LnUJqrV6/i7t276Natm+OxLl26AABOnDjheGzu3Llo3rw5YmNjK7eiCqWiC+jSpQtTqVTMw8ODabVaBoD16NGDmc1mx5hRo0axhg0bMpvNVuayevXqxebMmeO47+/vX2Qm+vXXXxkAdvbsWcYYYyNHjmQvvvhiqcsDwBYtWuS4/+OPPzIAjlq+pJmrd+/ejvsWi4Xp9XrHNsX06dOLjGeMseeff/6JM5enp2ep2w8PjR8/no0ePdpxv6Rtrkfj/fPPP5lSqWTffvut4/kLFy4wAOzcuXOMMca2bt3KmjRpUuo6U1JSGAAWGBjIvv32W3b48GHWtm1b1qNHj1Jf06NHDzZ16tQij61YsYI1bNjQcR8A02q1zMPDw3ELCgoqdZmHDx9mANiff/5Z5PFnn32WffbZZ4wxxo4dO8b8/f3Z77//zhgruWJ6kmo5c40bNw4ZGRk4ffo0du/ejZycHEyZMqXImJYtW0IuL/r21q9fj9atW8PHxweenp7Yu3cvsrOzi4x5tLb+r//6LwBwzCTnzp1D165dy4ytrNc/abxSqYSPj49j/KVLl9CqVasi48vaPnmocePG8PT0LPLYt99+i44dO6J27drw9PREcnJysfdelqtXr8JqtSIyMtLxWHBwMLy8vJCZmQkAePnll4ttCz3KbrcDABITExETE4P27dvjk08+wZ49e0qNhZXzpI3Vq1cjIyPDcfvuu+8AAIcOHYKnp6fjdujQoScu02QyYeTIkfjnP/8JLy+vcq2/JMpKv1Kgv/3tb2jUqBEAICgoCLNnz8bw4cOxePFi6HQ6AIBery/ymkOHDmHcuHH48MMP0blzZ9SoUQOTJk2CxWIpMk6lUjn+/7AUefhLUZ4fdFmvf9L4h695dH1llVmlefy9X716Ff3798fMmTOxZMkS1KpVC4sWLcLly5fLvczy/pKXpXbt2gCAJk2aOB57+P/s7Gw0aNCgxNc8/uV0586dYjsr6tat6/ideFTr1q2RkZHhuF+vXj388ssvAAq/9GrUqAEAsNls+O233+Dr64vc3FxkZmYWKQdtNhvS09Oxbt065OTklOv9VsuZ63FKpRI2m82xbVSSY8eOoVmzZnjzzTcRHh6Ohg0b4sqVKxVaT2hoKA4cOFDFaMsvKCgI//73v4s89vj98jh58iR0Oh3mzp2L1q1bo3Hjxrh27VqRMSqVCjabrdRlBAYGQqlU4ujRo47HLl68iPz8fAQHB5crjlatWkGlUhVJ6of/9/PzK/E1bdu2xYEDB4ok9759+9CuXbtyrVOn06FRo0aOm06nQ8OGDeHj44P9+/c7xqWnpwMoTMZ69erh7NmzRWbC1q1bY9y4cdi7d2+51gtU0+QqKCjAzZs3kZubi8OHD2P+/Pno2LEjatWqVeprAgMDkZmZiR07diAzMxOTJk3CzZs3K7TemTNnIi0tDe+88w4uXryI06dPY/ny5VV9O6UaO3Ysjhw5gkWLFuHnn3/GwoULcfbs2QrPZoGBgfjzzz+xbt06XL58GfPnzy+208Lf3x9Hjx7FL7/8gt9//73YMmrUqIH4+HhMmTIFhw4dwsmTJxEXF4cePXqgWbNmAICvvvqqzETz8vJCXFwc3n33XRw6dAhnzpzBhAkT0Lt3b9SvXx8AsHz5cnTv3t3xmuHDh6OgoABvvvkmLly4gEWLFuHo0aN47bXXiiz7999/x82bN4vcHq9KHpLL5Xj99dfx7rvvYt++ffjXv/6FyZMnY9iwYfD29oZKpUJISEiRm4eHB5599tkis+6TVMvkWrJkCerUqYN69eph4MCBaN68OVJTU8t8Tb9+/TBu3DiMGDEC7du3R40aNSq8F6hZs2bYvn07du/ejZYtW6JXr16O7Q0emjRpgvXr12P58uUIDw/H+fPnMWLECGg0mgotJzw8HO+//z7eeustRERE4Pr16xg/fnyRMdOnT0deXh4aNmyI8PDwEpfz0UcfoVOnToiNjUXnzp1Rr149bNiwwfH8H3/88cTP4+OPP0bPnj3x0ksvoXv37vD398fnn3/ueP7u3btFKgpvb2/s3LkTP/74I8LCwrB27Vps3ry5WBIPHjwYderUKXIra5Z/9913MWTIEAwaNAjdu3dHeHg4Vq5cWWbsFSVjziimyVPzwgsvoEmTJlixYoXoUMgTVMsdGlKyYsUKtG/fHp6enti0aRP27duHuXPnig6LlAMll4s7d+4c5s6di3v37iEoKAhbtmxB+/btRYdFyoHKQkI4qZY7NAipDii5COHE5be5NBoNnn32WdFhEFKiO3fulHpaj8sn17PPPlvuw01cyp9/AvXrAzk5QM2aoqMhnDz843dJqCzkxcMDOHKk8F8iSS4/c1VbCgXQvLnoKIhANHPx8uefgExW+C+RJJq5ePH0BLKzC/91Ervd7pRTP0j5yWSyYucFlhclFy8yWeGOjEqcj/U4s9mMrKysUo/yJnypVCr4+flBrVZX6HWUXLzcuwfUqgX88UeV9xZmZWWhRo0a8Pb2rtTJk6TyGGPIy8tDVlZWiSdjloWSi5caNQoT668zXSvLbrfDYrHA29sbSiX9uETw9vbGb7/9BrvdXqESkXZo8MJY4c6MKm4jPdzGohlLnIeffUW3d+mrkJf794EGDZxSFj4udtkPTl3eo7ZPostXOwvNXLzUrFk4a7nh0RlWqxVz585FcHAwmjdvjuDgYCQkJCA/Px/Xr18vdm3FJUuWlHkFrMrYsWMHgoOD0ahRIwwYMAD379936vKdgZKLF5sN+Omnwn/dzJgxY3D8+HEcOXIEP/30E86fP48ePXrgt99+eyrJdf/+fYwZMwZff/01Ll++jDp16uD999932vKdhcpCXgoKgOefd7tjCy9fvowvv/wSWVlZ+Nvf/gag8IIvgwYNAgDExMQgKysLYWFh8PPzQ+vWrfHrr79i4MCB0Gq1WLduHcLCwqoUw65du9C6dWvHdTQmTJiAmJgY/M///E/V3pyTUXLxUrOmWx6dcfLkSTRu3Bg+Pj4lPr9q1SpMnz69yGWhk5OTsXnzZoSEhJT4moEDB5Z6DcXt27cXu55hVlYW/P39HfcDAgLwyy+/VHhvHm+UXJz0XXIAjbMu4pJfMOwKJe0oKMPmzZsr/JrqsPeUkosTjcWMWcnvYsI7n8OgcJ+POSIiApcuXUJeXh68vb2dssyKzlx+fn7Yt2+f4/7169dRr149l5q1AEoubgxaPUbP2yo6DKd7uHduzJgxWLduHby8vMAYw4YNG9ChQwfUrFkTf/zxR5HXlPTYoyo6c0VHR2PixIm4ePEigoODsXLlSgwZMqRS74cnSi5O5DYrWv78b5wOagW7k2cu0SVmcnIy5s+fj3bt2kGpVIIxhs6dO+Oll16Cv78/mjRpgpCQEDRs2BDbtm3D5MmTMXr0aOj1eqfs0KhRowY+/fRT9OvXD1arFaGhoVi/fr2T3p3zuPzVn+rXr18tz0QetDgNH300HtOmrYZRo690QthsNvz8888ICgqCQqFwcpSkPMr6GZT1+0kzFydGjR4T397w5IHEbbnWFqAbUdis6HBqPxQ2q+hQiCCUXJworRb0258KpZXOwZIqKgs5MWl0mDF1legwiEA0c3GitFrQ88ftNHNJWLmSy2g0ol+/fggKCkJYWBiio6Nx/fp1AIWtL6Ojo9G4cWOEhITghx9KPx1izZo1aNy4MQIDA5GQkACr1X23RxQ2Kzpk0DaXlJW7LExISEDv3r0hk8mwfPlyJCQkIC0tDbNmzUJkZCR2796N48ePY+DAgbhy5Uqxs2avXbuGxMREnDp1Cr6+vujbty/WrFlTrAmbuzBpdJg9IYnPwje+wme5ADCs7CaCQOEpJwsWLMDGjRuhUChgs9nQuXNnfPDBB8jPz0daWhoSEhIc45csWYJhw4YV62NcFTt27MD06dNhtVrRsmVLrF+/vliTdQB45513sG3bNscu9P/+7//GK68Ufn5XrlxBQkIC8vLyYDQa8eKLL+LDDz902pEe5VqKVqtFTEyM43iuyMhIXL16FQCwadMmTJw4EQDQpk0b1K5du8TZa/PmzXj55ZdRu3ZtyGQyvPbaa0hJSXHKm3BFSosZffd9AaWl9D7N1VV1OuVkxowZjv7GO3fuxLhx4xxtaadPn46+ffs6+h6npaVh9+7dTouzUjs0li5ditjYWOTl5cFutxe5lntAQACysrKKvaakI5lLGpeUlISkpP9847viSXDloWB2BF//Cbs79oM7FYbV7ZQTLy8vx//v3bsHmUwGu93ueOzhYVkGgwEWiwV16tSpUmyPqnByLViwAJcuXcKqVatgMBiKHZ1c1gEfj44tbdzUqVMxdepUx/2yrsXtykxqLRbFzxMdhtNVx1NOli5dihUrViAnJwfJycmOA46XLFmC2NhY/POf/8Tvv/+OxMTEUvtBV0aFkmvx4sXYunUrvv/+e+j1euj1egCFnR4ezl43btyAn59fsdf6+fk5doKUNc5dKC1mDNrzOb7s8Sqsqopd705qeJ9yMnnyZEyePBmnT5/Gq6++ihdeeAHe3t5YvXo1RowYgRkzZuD27dvo1q0bIiMj0a1btwrHU5Jyb7klJSUhJSUFe/bsKTLVDho0yNH8+vjx47h58yY6dix+HN2AAQPw1Vdf4datW2CMYdWqVS55JLOzyBmDT/5tyF370M0Ke/SUE2cZOHAgwsLCSrxlZ2cXG//4F3V5Tzlp2bIl6tWrhwMHDgAonNFGjRoFAPD19UXv3r1x8OBBp72vciVXTk4Opk2bhvz8fERFRSEsLAzt2rUDACxatAg//vgjGjdujLi4OGzYsMGxp3Ds2LHYtm0bAKBhw4aYM2cOOnTogMDAQPj6+mLMmDFOeyOuxqzWYNmwWTCrNaJDcapHTznJz88HUFjif/bZZ7hy5UqlTzl5uFPh8dvjJSFQeMrJ8ePHcfHiRQAo85STCxcuOP5/5coVnDp1Cs2aNQNQ+Du5a9cuAEBBQQH27dtXaulaGeUqC+vXr1/qNlLt2rWRlpZW4nOffvppkfvjxo3DuHHjKhhi9aSymDBy+yf4LDYBFpV7JZgrn3Ly66+/IiYmBhkZGQCAWbNm4fLly1CpVFAqlVi+fDmaNm0KAFi/fj3eeOMNfPTRR7BYLOjXrx8GDhxYtQ/nEXTKCSf9k/YWSS465aT6olNOXIxFpcGa/pNEh0EEomMLOVGbTZi0cSHU5pL75RL3R8nFiV0mw10vX9irwVWKCB9UFnJiVamREhNf5eVUtgkAcZ7KNsOg5OJEYzZiyufvY8mr78Ck1lZ6OXK5HCqVynEps+pwvT538rA/l0qlqvABvZRcnNhkclwMaA6brOqVt5+fH7KysvDbb785ITJSUQ87S1YUJRcnVpUa33RzzhEoarUajRo1op7IAlBPZBekMRnw9pp3sGDM+zBpdE5ZpqtdUZaUjX5anNgUShwOi4LNjS5lTSqGfvKcWJUqpLWPFR0GEYhmLk40JgM+THoNGpNBdChEEEouTqxKFb6OegVWpUp0KEQQKgs5sSmUOBweJToMIhDNXJxoTQ+wYsEIaE0PRIdCBKHk4sSsVOPTl9+AWUmn+EsVlYWc2BVKnGraTnQYRCCauTjRGR9gbWJ/6IxUFkoVJRcnJpUaC+PnwkRXfpIsKgs5sSuUyHzOeRc7IdUPzVyc6AwFSJ3RCzpDgehQiCCUXJyYNFrMmLoKJk3lz+Ui1RuVhZzY5Qpk1XlOdBhEIJq5ONEZCrB9cicqCyWMkosTo0aHuDlbYHTSuVyk+ilXck2ePBkBAQGQyWQ4d+4cACA/P7/INb2DgoKgVCpLPBX9wIED0Ov1RcYbDO59tDiTyfBA6wFG17yQrHJtcw0cOBBvvfVWkQYLXl5ejksGA4UdUA4ePIhnnnmmxGU0a9asSFsZd6czPsCmmdEYvGg3DDoP0eEQAcqVXJ07d37imLVr15ba3U+KDFp9YWJp9aJDIYI4ZZvryJEjyMvLQ58+fUodk5mZiYiICLRp0wYrV650xmpdmowx6I0FkNEFZSTLKcmVnJyMkSNHFmsy/lBERARycnJw8uRJfPXVV1i1ahU2bdpU4tikpCTUr1/fcauubVu1JgPWzR4ALZ2JLFlVTq6CggKkpqYiPr70q8vWrFkTtWrVAlDYFWLo0KE4dOhQiWOnTp2KnJwcx62kDu3VgUHngdilh2h7S8KqnFxffvklWrRo4Wj+XJLc3FxHk+d79+5hx44dTu0964rkdhv8cq9BbreJDoUIUq7kmjhxoqMP0QsvvIBGjRo5nluzZk2JHSIf7Sq5ZcsWhIaGomXLloiMjESPHj0wevRoJ70F16QxGf+6QI1RdChEEGp+x0nssh+K3K9s8zvi2sr6/aQjNDiR26xocu0c5Dar6FCIIJRcnGgsZsxKfhcai1l0KEQQOiqeE4NWj9HztooOgwhEMxcncpsV4ReOUVkoYZRcnKitZoz9ajnUVioLpYrKQk6MGj0mvr1BdBhEIJq5OFHYrOhwaj8UVBZKFiUXJ0qrBf32p0JptYgOhQhCZSEnJo0OM6auEh0GEYhmLk6UVgt6/ridZi4Jo+TiRGGzokMGbXNJGZWFnJg0OsyekCQ6DCIQzVycKC1m9N33BZR0+JNkUXJxomB2BF//CQpmFx0KEYTKQk5Mai0Wxc8THQYRiGYuTpQWM4buTKayUMIouTiRMwaf/NuQu/a5qIQjKgs5Mas1WDZslugwiEA0c3GispgwZusyqCwm0aEQQSi5COGEykJOLCoN1vSfJDoMIhDNXJyozSZM2rgQajOVhVJFycWJXSbDXS9f2KmFkGRRWciJVaVGSkzpl/gm7o9mLk40ZiNmJidCY6Yr7kpVpTtLAkBAQACCg4Md3SJTU1NLXcb8+fMRGBiIwMBAJCYmVj1yF2eTyXExoDlsMvr+kqpKd5Z8aPPmzQgJCSnz9enp6UhJScGZM2egVCrRoUMHdOzYEb169apc1NWAVaXGN92GiA6DCFSur9XOnTujfv36lV5Jamoq4uLi4OHhAY1Gg/j4eKSkpFR6edWBxmTAnJVToaH+XJJV5Zpl+PDhCA0NxdixY3Hnzp0Sx2RlZcHf399xPyAgAFlZWVVdtUuzKZQ4HBYFm4L2GUlVlX7y6enp8PPzg8ViwT/+8Q+MGjUKO3fuLHGs7JFd0mU1VklKSkJS0n/O4K0unSUf72oCpQpp7WPFBENcQpVmLj8/PwCASqXClClTSu0W6efnh+vXrzvu37hxw/Hax7lLZ0mNyfBXfy4qC6Wq0slVUFCA/Px8x/2UlJRSu0UOGjQI69evR0FBAUwmE5KTkzFkiHtv7FuVKnwd9QqsSpXoUIggle4seevWLURFRaFFixYIDQ3FwYMH8dlnnzleExMTgxMnTgAAunbtisGDByM0NBRNmzZFz549ER0dzecduQibQonD4bTNJWXUWdJJHt/m0poe4KOPxmPatNUwavTUWdJNUWdJAcxKNT59+Q2YlWrRoRBBqGbhxK5Q4lTTdqLDIALRzMWJzvgAaxP7Q2d8IDoUIgglFycmlRoL4+fCpKKyUKqoLOTErlAi87myj7kk7o1mLk50hgKkzugFnaFAdChEEEouTkwaLWZMXQWTRis6FCIIlYWc2OUKZNV5TnQYRCBKLk50hgJsmhmNwYt2w6DzKPZHZvqjsvujspATo0aHuDlbYNToRIdCBKHk4oTJZHig9QCjqz9JFiUXJzrjA2yaGU1/RJYwSi5ODFp94faWVi86FCIIJRcnMsagNxZA5tonHRCOKLk40ZoMWDd7ALR0JrJk0a54Tgw6D8QuLfmyB0QaaObiRG63wS/3GuR2m+hQiCCUXJxoTMa/LlBDl7OWKioLOTHoPPDKh9+JDoMIRDMXJ3KbFU2unYPcZhUdChGEkosTjcWMWcnvQmMxiw6FCEJlIScGrR6j520VHQYRiGYuTuQ2K8IvHKOyUMIouThRW80Y+9VyqK1UFkoVlYWcGDV6THx7g+gwiECV7ixpNBrRr18/BAUFISwsDNHR0UWaLTzqwIED0Ov1jg6UYWFhMBjc+7Aghc2KDqf2Q0FloWSVK7kGDhyIH374oUiPLQBISEhAZmYmMjIy0KdPHyQkJJS6jGbNmiEjI8Nx0+nc+yRCpdWCfvtTobRaRIdCBKl0Z0mtVouYmBhH363IyEhcvXrV+RFWUyaN7q8L1Lj3lwgpndN2aCxduhSxsaU3e8vMzERERATatGmDlStXOmu1LktptaDnj9tp5pIwp+zQWLBgAS5duoRVq1aV+HxERARycnJQq1Yt5OTkICYmBj4+Phg8eHCxsdW1s+TjFDYrOmTsx8FWL1CPLomq8sy1ePFibN26Fbt27YJeX/JZtzVr1kStWrUAFLZcGTp0aKldKN2ls6RJo8PsCUlUFkpYlZIrKSkJKSkp2LNnD7y8vEodl5ubC7vdDgC4d+8eduzYUWoXSnehtJjRd98XUNLhT5JV6c6SOTk5mDZtGvLz8xEVFYWwsDC0a/efljljx47Ftm3bAABbtmxBaGgoWrZsicjISPTo0QOjR4/m845chILZEXz9JyiYXXQoRBDqLOkkj1/080nooqDugTpLCqC0mDF0ZzKVhRJGycWJnDH45N+G3LULA8IRHVvIiVmtwbJhs0SHQQSimYsTlcWEMVuXQWUxiQ6FCELJRQgnVBZyYlFpsKb/JNFhEIFo5uJEbTZh0saFUJupLJQqSi5O7DIZ7nr5wk4thCSLykJOrCo1UmLiRYdBBKKZixON2YiZyYnQmOmKu1JFycWJTSbHxYDmsMnoI5YqKgs5sarU+KbbENFhEIHoa5UTjcmAOSunQkP9uSSLkosTm0KJw2FRsCmoOJAq+slzYlWqkNa+9GuKEPdHMxcnGpPhr/5cVBZKFSUXJ1alCl9HvUIXp5EwKgs5sSmUOBweJToMIhDNXJxoTQ+wYsEIaE0PRIdCBKHk4sSsVOPTl9+AWakWHQoRhMpCTuwKJU41bffkgcRt0czFic74AGsT+0NnpLJQqii5ODGp1FgYPxcmFZWFUkVlISd2hRKZz4WIDoMIRDMXJzpDAVJn9ILOUCA6FCIIJRcnJo32r/5cWtGhEEEq3bYVAC5duoT27dsjKCgIbdu2xfnz50tdxvz58xEYGIjAwEAkJiZWPXIXZ5crkFXnOdjlCtGhEEGq1LZ1/PjxSEhIwM8//4y33noLY8aMKfH16enpSElJwZkzZ3D+/Hns2rUL3333XdWjd2E6QwG2T+5EZaGEVbpt6+3bt3Hy5Em8+uqrAIABAwbg2rVrJTYdT01NRVxcHDw8PKDRaBAfH4+UlJSqR+/CjBod4uZsgZH6c0lWpbe5srOzUbduXSiVhTscZTIZ/Pz8kJWVVWxsVlZWkVkvICCgxHFAYc+v+t4CZ7YAABBOSURBVPXrO27VtbMkk8nwQOsBRld/kqwq7dCQPfaLU1Y3okfHljXOXTpL6owPsGlmNP0RWcIqnVwNGjRATk4OrFYrgMKEyc7Ohp+fX7Gxfn5+RcrFGzdulDjOnRi0egxetBsGbcmtbIn7q3Ry+fr6Ijw8HJ9//jmAwu6RAQEBCAgIKDZ20KBBWL9+PQoKCmAymZCcnIwhQ9z74i0yxqA3FkBGLYQkq9JtWwFg9erVWL16NYKCgrBw4UKsWbPG8ZqYmBicOHECANC1a1cMHjwYoaGhaNq0KXr27Ino6GgOb8d1aE0GrJs9AFo6E1myqG2rk1DbVmmitq0CyO02+OVeg9xuEx0KEYSSixONyfjXBWroctZSRUfFc2LQeeCVD937KBRSNpq5OJHbrGhy7RzkNqvoUIgglFycaCxmzEp+FxqLWXQoRBAqCzkxaPUYPW+r6DCIQDRzcSK3WRF+4RiVhRJGycWJ2mrG2K+WQ22lslCqqCyspCf90dio0WPi2xueUjTEFdHMxYnCZkWHU/uhoLJQsii5OFFaLei3PxVKq0V0KEQQKgs5MWl0mDF1legwiEA0c3GitFrQ88ftNHNJGCUXJwqbFR0yaJtLyqgs5MSk0WH2hCTRYRCBaObiRGkxo+++L6Ckw58ki5KLEwWzI/j6T1Awu+hQiCBUFnJiUmuxKH6e6DCIQDRzcaK0mDF0ZzKVhRJGycWJnDH45N+G3LUvUUI4orKQE7Nag2XDZokOgwhEMxcnKosJY7Yug8piEh0KEYSSixBOqCzkxKLSYE3/SaLDIALRzMWJ2mzCpI0LoTZTWShVVZ658vPz0bVrV8f9Bw8e4OrVq7h9+zaeeeYZx+MHDhxATEwMgoKCHI8dOXIEOp179q+yy2S46+ULO7UQkqwqJ5eXlxcyMjIc9xcvXoyDBw8WSayHmjVr5rh+vLuzqtRIiYkXHQYRyOll4dq1a0tt3yolGrMRM5MToTHTFXelyqnJdeTIEeTl5aFPnz4lPp+ZmYmIiAi0adMGK1eudOaqXY5NJsfFgOawyWizVqqcurcwOTkZI0eOdLRyfVRERARycnJQq1Yt5OTkICYmBj4+Phg8eHCRcUlJSUhK+s+pGtW1batVpcY33dy7Bxkpm9O+VgsKCpCamor4+JK3M2rWrIlatWoBKGy7MnToUBw6dKjYOHdp26oxGTBn5VRoqD+XZDktub788ku0aNECwcHBJT6fm5sLu73w9It79+5hx44dCA8Pd9bqXY5NocThsCjYFPSnRKlyWnKtWbOm2I6MsWPHYtu2bQAK27qGhoaiZcuWiIyMRI8ePTB69Ghnrd7lWJUqpLWPhVWpEh0KEYQ6S1bSky4KqjEZMH/F3/GPif8PJk3xv+VRZ0n3QJ0lBbAqVfg66hWauSSMNgg4sSmUOBweJToMIhDNXJxoTQ+wYsEIaE0PRIdCBKHk4sSsVOPTl9+AWakWHQoRhMpCTuwKJU41bSc6DCIQzVyc6IwPsDaxP3RGKgulipKLE5NKjYXxc2FSUVkoVVQWcmJXKJH5XIjoMIhANHNxojMUIHVGL+gMBaJDIYJQcnFi0mgxY+oqmDRa0aEQQags5MQuVyCrznOiwyAC0czFic5QgO2TO1FZKGGUXJwYNTrEzdkCYwkH7RJpoOTihMlkeKD1AKOrP0kWJRcnOuMDbJoZTX9EljDaocGJQavH4EW7YdDqnzj28XPD6Fwv90AzFycyxqA3FkDm2ueiEo4ouTjRmgxYN3sAtHSBGsmispATg84DsUuLX92KSAfNXJzI7Tb45V6D3G4THQoRhJKLE43JiA+TXoPGRJezlioqCzkx6DzwyoffiQ6DCEQzFydymxVNrp2D3GYVHQoRhJKLE43FjFnJ70JjMYsOhQhCZSEnBq0eo+dtFR0GEcgpM1dAQACCg4MRFhaGsLAwpKamljhu/vz5CAwMRGBgIBITE52xapclt1kRfuEYlYUS5rSZa/PmzQgJKf209vT0dKSkpODMmTNQKpXo0KEDOnbsiF69ejkrBJeitpox9qvlmDZtNYzUjEGSnto2V2pqKuLi4uDh4QGNRoP4+HikpKQ8rdU/dUaNHhPf3gCj5snHFhL35LTkGj58OEJDQzF27FjcuXOn2PNZWVnw9/d33A8ICEBWVpazVu9yFDYrOpzaDwWVhZLllORKT0/H6dOncfLkSXh7e2PUqFEljpM9cm5Tac1VkpKSUL9+fcfNlTpLxi77wXF7EqXVgn77U6G0Wp5CZMQVOSW5/Pz8AAAqlQpTpkwpsWOkn58frl+/7rh/48YNx+se5S6dJU0a3V8XqKEzkaWqyslVUFCA/Px8x/2UlJQSO0YOGjQI69evR0FBAUwmE5KTkzFkiPv2DFZaLej543aauSSsysl169YtREVFoUWLFggNDcXBgwfx2WefAQBiYmJw4sQJAEDXrl0xePBghIaGomnTpujZsyeio6OrunqXpbBZ0SGDtrmkjDpLVkB5trXK69GzjelM5OqLOksKoLSY0XffF1DS4U+SRcnFiYLZEXz9JyiYXXQoRBA6dIATk1qLRfHzRIdBBKKZixOlxYyhO5OpLJQwSi5O5IzBJ/825K69v4hwRGUhJ2a1BsuGzRIdBhGIZi5OVBYTxmxdBpXFJDoUIgglFyGcUFnIiUWlwZr+k0SHQQSimYsTtdmESRsXQm2mslCqKLk4sctkuOvlCzu1EJIsKgs5sarUSImJFx0GEYhmLk40ZiNmJidCY6Yr7koVJRcnNpkcFwOawyajj1iqqCzkxKpS45tu7nsyKHky+lrlRGMyYM7KqdBQfy7JouTixKZQ4nBYFGx0zULJop88J1alCmntY0WHQQSimYsTjcnwV38uKgulipKLE6tSha+jXoFVqRIdChGEykJObAolDodHiQ6DCEQzFyda0wOsWDACWtMD0aEQQWjm4sSsVOPTl9+AWaku8XlnXqaNuCZKLk7sCiVONW0nOgwiEJWFnOiMD7A2sT90RioLparKyWU0GtGvXz8EBQUhLCwM0dHRRRouPHTgwAHo9XpH98mwsDAYDO67m9qkUmNh/FyYVCWXhcT9OaUsTEhIQO/evSGTybB8+XIkJCQgLS2t2LhmzZo5rh3v7uwKJTKfK73TJnF/VZ65tFotYmJiHL23IiMjcfXq1SoHVt3pDAVIndELOkOB6FCIIE7f5lq6dCliY0s+7CczMxMRERFo06YNVq5c6exVuxSTRvtXfy6t6FCIIE7dW7hgwQJcunQJq1atKvZcREQEcnJyUKtWLeTk5CAmJgY+Pj4YPHhwkXFJSUlISkpy3HelzpIVYZcrkFXnOdFhEIGcNnMtXrwYW7duxa5du6DXF2+yXbNmTdSqVQtAYduVoUOHltiB0l06S+oMBdg+uROVhRLmlORKSkpCSkoK9uzZAy8vrxLH5Obmwm4v7Phx79497Nixo8QOlO7CqNEhbs4WGKltq2RVOblycnIwbdo05OfnIyoqCmFhYWjXrvCPp2PHjsW2bdsAAFu2bEFoaChatmyJyMhI9OjRA6NHj67q6l0Wk8nwQOsBRld/kizqLFkBFTlkSWcowKaZ0Ri8aDcMOo8KrYc6S1Yf1FlSAINWX5hY2uLbn0QaKLk4kTEGvbEAMtcuDAhHdOAuJ1qTAetmD6hUWUg42/hK6c8NS3Xaaii5ODHoPBC7tPifGoh0UFnIidxug1/uNcjtNtGhEEEouTjRmIx/XaCGLmctVVQWcmLQeeCVD78THQYRiGYuTuQ2K5pcOwe5zSo6FCIIJRcnGosZs5LfhcZiFh0KEYTKQk4MWj1Gz9sqOgwiEM1cnMhtVoRfOEZloYRRcnGitpox9qvlUFupLJQqKgs5MWr0mPj2BtFhEIFo5uJEYbOiw6n9UFBZKFmUXJworRb0258KpdUiOhQiCJWFnJg0OsyYWvxaIkQ6JJVcj5/syPOkRKXVgm7/2o19baOrdxuhR48gf9IR448fbe7EI8yrtOyyjoLniMpCThQ2Kzpk0DaXlElq5nqaTBodZk9IevJA4rZo5uJEaTGj774voKTDnySLkosTBbMj+PpPUDC76FCIIFQWcmJSa7Eofp7oMIhANHNxorSYMXRnMpWFEkbJxYmcMfjk34acrv4kWVQWcmJWa7Bs2CzRYRCBaObiRGUxYczWZVBZTKJDIYI4JbkuXbqE9u3bIygoCG3btsX58+dLHDd//nwEBgYiMDAQiYmJzlg1IS7LKck1fvx4JCQk4Oeff8Zbb72FMWPGFBuTnp6OlJQUnDlzBufPn8euXbvw3XfuewEXi0qDNf0nwaLSiA6FCFLl5Lp9+zZOnjyJV199FQAwYMAAXLt2rVjT8dTUVMTFxcHDwwMajQbx8fFISUmp6updltpswqSNC6E2U1koVVXeoZGdnY26detCqSxclEwmg5+fH7KyshAQEOAYl5WVhS5dujjuBwQEYPPmzcWW93hnyZs3b6J+/fqlrv/+/fuVbpBXf1GlXlYu9+/fxzFPT+D/jarwa3nHVemGgm+V/nOo6vgKx1XRWCq53CfFdefOnVKfc8reQtljPahK60r06LjSxkydOhVTp04t97pdqcXQoyiuinHHuKpcFjZo0AA5OTmwWguP/maMITs7G35+fkXG+fn5FSkVb9y4UWwMIe6kysnl6+uL8PBwfP755wAKO0gGBAQUKQkBYNCgQVi/fj0KCgpgMpmQnJyMIUOGVHX1hLgsxXvvvfdeVRfy/PPPIzExER988AGOHz+OtWvXwtfXFzExMWjSpAnq1q2LgIAA3L17F+PHj8eKFSvQt29fjB8/3glvoXD9rojiqhh3i8vl27YSUl3RERqEcELJRQgnbpFccXFxqF+/PsLCwhAWFoYZM2YIi6W8h4I9bQEBAQgODnZ8RqmpTrx4TAVMnjwZAQEBkMlkOHfunONx0Z9baXFV6XNjbmDUqFFs2bJlosNgjDEWFRXF1q5dyxhj7Msvv2SRkZFiA/qLv78/O3v2rOgw2MGDB1l2dnaxeER/bqXFVZXPzS1mLldR3kPBpKxz587Fjrhxhc+tpLiqym2SKykpCS1atECfPn2QkZEhJIayDgVzBcOHD0doaCjGjh1b5mE7T5u7fm7VIrk6deoEHx+fEm/Z2dl4//33cfnyZZw5cwZjxoxB7969cf/+fSGxlvdQsKctPT0dp0+fxsmTJ+Ht7Y1Royp+zCNPbvm5OadidS1BQUHsxIkTT329t27dYjVr1mQWi4Uxxpjdbme1a9dm165de+qxlOXXX39lnp6eQmN4dFvGlT63sraxKvq5VYuZ60kePbDy6NGjyMvLQ6NGjZ56HOU9FOxpKygoQH5+vuN+SkoKwsPDBUZUlNt+bs7KeJG6d+/OQkJCWMuWLVlkZCTbt2+fsFguXrzIIiMjWePGjVmrVq3YuXPnhMXy0JUrV1hYWBgLDQ1lISEh7KWXXhI2m06YMIHVq1ePKRQKVrt2bRYYGMgYE/+5lRRXVT83OvyJEE7coiwkxBVRchHCCSUXIZxQchHCCSUXIZxQchHCCSUXIZxQchHCCSUXIZz8f17mRNkCCKXsAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 216x400 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAEWCAYAAABPON1ZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAVUklEQVR4nO3de5gldX3n8fcHRgVkWBBGDcjQiMZoWMFkUAmbxACJwCDG6JMg4moSd2ISV3BxdTTGeyIx0Rg3JkpUdBWi4SJeZr1gDKvso8ThFkDQdWG4gwhyDSLgd/+oauZM0z1zpqdPn990v1/P08/MqapT9a3L+Zxf/arOOakqJEnt2mbcBUiSNs6glqTGGdSS1DiDWpIaZ1BLUuMMaklqnEE9gyQTSSrJkv7xOUleMe66tHFJ1iU5dNx1zKckT0lyYZK7krw6yceSvHNEy/pikpeNYt6zkWRJ/zqdGHcto7Tgg7p/4d6b5O4kNyc5OcmOY6rlyUl+nOSTM4w/uT/onjQw7FVJ1ia5L8nHpkw/+WZy98Dfn454NabWvFOS9yW5pl/+9/vHu/XjNwjOqW+AW4t0/iLJrf3fu5Nkhml/LcklSW7vp/1Mkj0Gxj8myaeT/LD/OyXJTltQ3uuAc6pqaVW9fwvms4Ekb516rFbV4VX18blaxsCynpPkp1OO5c/P9XK2Vgs+qHvPq6odgV8ADgDeNKY6PgB8e7oRSf4TsM80o24A3gl8dCPz3bmqduz/3rHlZQ4nySOBfwZ+HjgM2An4JeBW4JnzVcc8WQX8JrAf8HTgSOAPZpj2O8Bzq2pnYHfg/wJ/PzD+ncAuwBPp9vnjgLduQW17AZdtwfNbccPAcbxjVT1v3AW1YrEENQBVdT3wRWBfmLa197AWxDD6553et5LuSnJBkv2mTHM0cDtdsE19/hLgfwCvmqbmM6vqLLrwmxNJXp7kyr7Wq5K8ZJaz+s/AcuAFVfWdqvppVf2gqt5RVf8rySf68Z/vW0ivA77eP/f2ftiBU2rbvT8DeszAsGf0Lc9HJNknydf6lupka3TnGdZzgy6AvtV23ZRlnZHkln47vHoj6/oy4D1VdV1/HL0HePl0E1bVzVV1w8CgB4EnDTzeGzirqu6sqjuAz9C92W22JF8Dfg342357/uw00/yX/kzntiSfS7L7wLi/SXJtkjuTnJ/kl/vhhwFvBH6nn+/F/fCHugD74+jcJH+V5Ef9Njx8YN57J/l6f5x9NckHZvn6OjDJt/ozlBuTvD/JI2aY9sgkl/fLvC7JawbGHZXk4n4+5ybZd3NrGZdFFdRJ9gSOAC4cweyfD5wGPAY4FThr8mDqT2vfDpwww3NfA3y9qv5tlsu+uj8oT07f5TCTJI8G3g8cXlVL6VrAF81yuYcCX6qqu6cbWVUvBa6hP6OpqncDv9KPnjwL+OaU59wAfBN44cDgY4DTq+p+IMC76FqqTwX2ZBat0STbAJ8HLgb2AA4Bjk/y3Bme8vP9tJMuZiPhmmR5ktuBe4HXAu8eGP0B4MgkuyTZhW5dv7i56wBQVQcD3wBe1W/P702p42C67fXbwM8AVwOfGpjk28D+rD9uT0uyXVV9Cfhz4NP9fDdoeAx4FvBdYLd+HT+SPNQldCrwr8CudPvopbNZR+AB4Lh+GQfRnb3NdDZzMvD7/bH9dOB/AyQ5APgH4BV9PR8FPtufFTZvsQT1Wf2L5ly6HffnI1jG+VU1GSbvBbYDnt2Pewfwkaq6duqT+jePPwDePItl/pCuK2cv4BeBpcApQzzvp8C+SbavqhuraranzbsCN87yuRtzKvBi6PqGgaP7YVTV96vq7Kq6r6puodvWvzqLZRwALKuqt1fVT6rqSroX8tEzTL8jcMfA4zuAHQdCaQNVdU3f9bEbXVfbFQOjLwAeSXeWdCtdi/vvZrEOw3gJ8NGquqCq7gPeAByY/uJbVX2yqm6tqgeq6j3Ao4CnbMb8r66qf6iqB4GP070ZPC7Jcrpt/OZ++54LfG4T89q9b+1O/v12X+O3q+q8vsYrgZOYeZ/fDzwtydKquq2qLuiHrwL+rp/Xg1U12ZV4wGas69gslqD+zarauar2qqo/qqp7R7CMh0K4qn4KXEd34O1P1/L86xme9z7g7f0p8Gapqruram1/AN9M13XyG9nIhamqugf4HeCVwI1J1iT5uemmzYYXdpZPM8mtdC/MuXY6XZjsTtcCL7pWI0kem+RTSa5PcifwSbow3Fx7MSUY6E71HzfD9HfT9cFP2gm4uzbxrWZVdRtdgH026y+gngZ8j+6NdSfg//Xr8TBJPjiwD9445LoN2p2uFT1Zz910+22Pfv4n9F0Fd/Tb4D+wedvzpoF5/3v/3x375d42MAwGXiMzuKF/nU7+/VNf48/1x+lN/T5/+0ZqfAFwFHBN303zrH74XsDrp+zvn6HfDq1bLEE9k3uAHQYeP34L5rXn5H/60+on0F0IfA4wQXfg3ER3GvzCJJPv9IcAf9kfhJMH/TeTHDOLGiZDY9pW3kMTVX25qn6d7kC9gq4lOd10gxd2rplmkq8Cz+27UzZV00yPp1vu7cBX6E7XjwH+cSAQ39XP4+lVtRNwLDOv78b277XAVVOCYWlVHTHDvC6ju5A4aT+Gv4C3BHgs64N+P+BDVXVPH5wfpOuSe5iqeuXAPpjNmeANdCEFPNT1tStwfd8f/Xq67bxLfwZwB+u355Z8teaNwGOSDG7/PWeaeBM+BFwKPKnf529mhn3et7yPotveX2B9N8+1wNum7O8dJt8MWrfYg/oi4Oh0F6lWAC/agnn9YpLf6ltNxwP3Ad+iO03bh64fcH+6F+UaYLIv9GfpXriT4wGeR3eBafI+0e2AbYFtk2yX9fd2PyvdPbTbJNmVru/5nI21zpM8rr+o8ui+xrvpTr1n4xN0L4Az+lbPNkl2TfLGJJPBczPd3Q2TbqHrenkiG3cq3cXKF/b/n7S0r/n2dLe8/feNzOMi4Ih0t8M9nm6/TPpX4M4kr0+yfZJtk+zb92VO538C/y3JHn1L/wTgY9NN2B8Hk/tlGV33zIV96xq6fuFX9Mvdnu60/OLp5jUHTgV+N8n+SR5F1+13XlWto9uWD9DtkyVJ3syGZw03AxN9w2OzVNXVwFrgrUkeme6i8Wzv4lhK9wZyT5KnMkP/dL89j0myU98FeRfrj+2TgD9OckA6OyZ53iYaGc1Y7EH9p3Qh+iPgbWwYCJvrs3RdCj+iu2jyW1V1f1X9e1XdNPlHFzI/7vtX6e+SGBwP8MOB7pk30V2QWk3XeryX9bcXPhH4Et0BeSld8L54E3VuQxcyNwC30fX1/dFsVrjv8zyUrlV+NnAnXQDuBpzXT/Yu4E396eZr+1PhPwP+Tz/s2dPMGrr+zCcDN1fVYIi9je42yzvo3vDO3EiJn6ALwHV0LfRPD9T+IF1w7A9cRdff/2G6U//pfIju4uMldNt6TT8MeKib6Jf7h3uwfr9cQvfG9IKBef0e3VnWdcD1dPvx5RtZj1mrqn+mO87PoGvl7sP6fvgv013E/B5d98iP2bB74rT+31sHzgA3x0uAA+m6Wt5Jt/3vm8V8TqC76+Yuum3+6Y1M+zK6i+t3Ar9PfwGzqs4D/pDuNskf0a3zsbOoZSyyiS42DSHJW+lOy7aaHS/NtySfBq6oqreMu5atzWJvUUsakb6bYZ++C+gwultYzxp3XVujrepjvBpekmnvbaa7f/ob81qMFqvH03VN7UrXzfOHVTWKzzAseHZ9SFLj7PqQpMaNpOtjt912q4mJiVHMWpIWpPPPP/+HVbVsunEjCeqJiQnWrl07illL0oKU5OqZxtn1IUmNM6glqXEGtSQ1zqCWpMYZ1JLUOINakhpnUEtS4wxqSWqcQS1JjfPb87RVm1i9Ztrh605cOc+VSKNji1qSGmdQS1LjDGpJapxBLUmNM6glqXEGtSQ1zqCWpMYZ1JLUOINakhpnUEtS4wxqSWqcQS1JjTOoJalxBrUkNc6glqTGGdSS1DiDWpIaZ1BLUuMMaklq3FBBneQ1SS5LcmmSf0yy3agLkyR1NhnUSfYAXg2sqKp9gW2Bo0ddmCSpM2zXxxJg+yRLgB2AG0ZXkiRp0CaDuqquB/4KuAa4Ebijqr4ydbokq5KsTbL2lltumftKJWmRGqbrYxfg+cDewO7Ao5McO3W6qjqpqlZU1Yply5bNfaWStEgN0/VxKHBVVd1SVfcDZwK/NNqyJEmThgnqa4BnJ9khSYBDgMtHW5YkadIwfdTnAacDFwCX9M85acR1SZJ6S4aZqKreArxlxLVIkqbhJxMlqXEGtSQ1zqCWpMYZ1JLUOINakhpnUEtS4wxqSWqcQS1JjTOoJalxBrUkNc6glqTGGdSS1DiDWpIaZ1BLUuMMaklqnEEtSY0zqCWpcUP9wos01yZWr5l2+LoTVzp/aQpb1JLUOINakhpnUEtS4wxqSWqcQS1JjTOoJalxBrUkNc6glqTGGdSS1DiDWpIaZ1BLUuMMaklqnEEtSY0zqCWpcQa1JDXOoJakxhnUktQ4g1qSGmdQS1LjhgrqJDsnOT3JFUkuT3LgqAuTJHWG/XHbvwG+VFUvSvJIYIcR1iRJGrDJoE6yE/ArwMsBquonwE9GW5YkadIwLeonArcAJyfZDzgfOK6q7hmcKMkqYBXA8uXL57pOaaswsXrNtMPXnbhynivRQjJMH/US4BeAv6+qZwD3AKunTlRVJ1XViqpasWzZsjkuU5IWr2GC+jrguqo6r398Ol1wS5LmwSaDuqpuAq5N8pR+0CHAd0ZalSTpIcPe9fFfgVP6Oz6uBH53dCVJkgYNFdRVdRGwYsS1SJKm4ScTJalxBrUkNc6glqTGGdSS1DiDWpIaZ1BLUuMMaklqnEEtSY0zqCWpcQa1JDXOoJakxhnUktQ4g1qSGmdQS1LjDGpJapxBLUmNM6glqXHD/hSXtKBNrF4z7fB1J66c50qkh7NFLUmNM6glqXEGtSQ1zqCWpMYZ1JLUOINakhpnUEtS4wxqSWqcQS1JjTOoJalxBrUkNc6glqTGGdSS1DiDWpIaZ1BLUuMMaklqnEEtSY0zqCWpcQa1JDVu6KBOsm2SC5N8YZQFSZI2tDkt6uOAy0dViCRpekMFdZInACuBD4+2HEnSVEuGnO59wOuApTNNkGQVsApg+fLlW16ZxmZi9Zpph687ceU8VyIJhmhRJzkS+EFVnb+x6arqpKpaUVUrli1bNmcFStJiN0zXx0HAUUnWAZ8CDk7yyZFWJUl6yCaDuqreUFVPqKoJ4Gjga1V17MgrkyQB3kctSc0b9mIiAFV1DnDOSCqRJE3LFrUkNc6glqTGGdSS1DiDWpIaZ1BLUuMMaklqnEEtSY0zqCWpcQa1JDXOoJakxhnUktQ4g1qSGmdQS1LjDGpJapxBLUmNM6glqXEGtSQ1brN+4UUatYnVa5qa/6jrmY2Zalp34sp5rkTzxRa1JDXOoJakxhnUktQ4g1qSGmdQS1LjDGpJapxBLUmNM6glqXEGtSQ1zqCWpMYZ1JLUOINakhpnUEtS4wxqSWqcQS1JjTOoJalxBrUkNc6glqTGGdSS1LhNBnWSPZP8S5LLk1yW5Lj5KEyS1Bnmx20fAE6oqguSLAXOT3J2VX1nxLVJkhiiRV1VN1bVBf3/7wIuB/YYdWGSpM4wLeqHJJkAngGcN824VcAqgOXLl89BaZorE6vXTDt83YkrxzKfxWic2879tvUb+mJikh2BM4Djq+rOqeOr6qSqWlFVK5YtWzaXNUrSojZUUCd5BF1In1JVZ462JEnSoGHu+gjwEeDyqnrv6EuSJA0apkV9EPBS4OAkF/V/R4y4LklSb5MXE6vqXCDzUIskaRp+MlGSGmdQS1LjDGpJapxBLUmNM6glqXEGtSQ1zqCWpMYZ1JLUOINakhpnUEtS4wxqSWqcQS1JjTOoJalxBrUkNc6glqTGGdSS1DiDWpIat8lfeFF7JlavmXb4uhNXznMlnZnq0abN5bZzPyxctqglqXEGtSQ1zqCWpMYZ1JLUOINakhpnUEtS4wxqSWqcQS1JjTOoJalxBrUkNc6glqTGGdSS1DiDWpIaZ1BLUuMMaklqnEEtSY0zqCWpcQa1JDXOoJakxg0V1EkOS/LdJN9PsnrURUmS1ttkUCfZFvgAcDjwNODFSZ426sIkSZ1hWtTPBL5fVVdW1U+ATwHPH21ZkqRJS4aYZg/g2oHH1wHPmjpRklXAqv7h3Um+u+XlzavdgB+Ou4gtkb/Y7Om3+nWeBde5t7nHy1Zma9zPe800YpigzjTD6mEDqk4CTtqMopqSZG1VrRh3HfPJdV4cXOet3zBdH9cBew48fgJww2jKkSRNNUxQfxt4cpK9kzwSOBr43GjLkiRN2mTXR1U9kORVwJeBbYGPVtVlI69s/m213TZbwHVeHFznrVyqHtbdLElqiJ9MlKTGGdSS1DiDehpJXpukkuw27lpGLclfJrkiyb8l+UySncdd0ygstq9BSLJnkn9JcnmSy5IcN+6a5kuSbZNcmOQL465lrhjUUyTZE/h14Jpx1zJPzgb2raqnA98D3jDmeubcIv0ahAeAE6rqqcCzgT9eBOs86Tjg8nEXMZcM6of7a+B1TPOhnoWoqr5SVQ/0D79Fd5/8QrPovgahqm6sqgv6/99FF1x7jLeq0UvyBGAl8OFx1zKXDOoBSY4Crq+qi8ddy5j8HvDFcRcxAtN9DcKCD61JSSaAZwDnjbeSefE+uobWT8ddyFwa5iPkC0qSrwKPn2bUnwBvBH5jfisavY2tc1V9tp/mT+hOl0+Zz9rmyVBfg7AQJdkROAM4vqruHHc9o5TkSOAHVXV+kueMu565tOiCuqoOnW54kv8I7A1cnAS6LoALkjyzqm6axxLn3EzrPCnJy4AjgUNqYd5Yvyi/BiHJI+hC+pSqOnPc9cyDg4CjkhwBbAfslOSTVXXsmOvaYn7gZQZJ1gErqmpr+wauzZLkMOC9wK9W1S3jrmcUkiyhu1B6CHA93dciHLNAP2ELQLrWxseB26rq+HHXM9/6FvVrq+rIcdcyF+yj1t8CS4Gzk1yU5IPjLmiu9RdLJ78G4XLgnxZySPcOAl4KHNzv14v6lqa2QraoJalxtqglqXEGtSQ1zqCWpMYZ1JLUOINakhpnUGtBSfJgfyvapUlOS7LDJqb/WJIX9f8/J8mC+UFULRwGtRaae6tq/6raF/gJ8MpxFyRtKYNaC9k3gCclmUhy6eTA/vvG3zq+sqTNY1BrQeo/Nn44cMm4a5G2lEGthWb7JBcBa+l+/OEjY65H2mKL7tvztODdW1X7Dw5I8gAbNkq2m9+SpC1ji1qLwc3AY5PsmuRRdF/pKm01bFFrwauq+5O8ne4XTq4CrhhzSdJm8dvzJKlxdn1IUuMMaklqnEEtSY0zqCWpcQa1JDXOoJakxhnUktS4/w9njgouICx4+gAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"l = []\n",
"\n",
"if load:\n",
" CLs_values = -1*np.array(CLs_list)\n",
"\n",
"if not os.path.exists('data/CLs/plots'):\n",
" os.mkdir('data/CLs/plots')\n",
" print(\"Directory \" , 'data/CLs/plots' , \" Created \")\n",
"\n",
"print(np.shape(CLs_values))\n",
"\n",
"figure(num=None, figsize=(2.7, 5), dpi=80, facecolor='w', edgecolor='k')\n",
"\n",
"for step in range(1,ste):\n",
" plt.clf()\n",
" plt.title('Branching ratio: {:.1E}'.format(BR_steps[step]))\n",
" plt.hist(CLs_values[0], bins = 40, range = (-5, 15), label = 'Ctt = 0', alpha = 0.8)\n",
" plt.hist(CLs_values[step], bins = 40, range = (-5, 15), label = 'Ctt = {:.2f}'.format(Ctt_steps[step]), alpha = 0.7)\n",
"# plt.hist(CLs_values[0][np.where(np.array(CLs_values[0]) > np.mean(CLs_values[0]))[0]], bins = 40, range = (-5, 15), alpha = 0.9)\n",
"# plt.hist(CLs_values[step][np.where(np.array(CLs_values[step]) < np.mean(CLs_values[0]))[0]], bins = 40, range = (-5, 15), alpha = 0.9)\n",
" plt.axvline(x=np.mean(CLs_values[0]),color='red', linewidth=1.0, linestyle = 'dotted')\n",
" plt.legend()\n",
" plt.savefig('data/CLs/plots/CLs-BR({:.1E}).png'.format(BR_steps[step]))\n",
" \n",
" l.append(len(np.where(np.array(CLs_values[step]) < np.mean(CLs_values[0]))[0])) \n",
"\n",
"figure() \n",
" \n",
"for step in range(2*ste):\n",
" if step%2 == 0:\n",
" floaty = True\n",
" else:\n",
" floaty = False\n",
" for key in pull_dic.keys():\n",
" if not os.path.exists('data/CLs/plots/{}'.format(key)):\n",
" os.mkdir('data/CLs/plots/{}'.format(key))\n",
" plt.clf()\n",
" plt.title('Pull {} - Ctt value {:.2f} - floating {}'.format(key, Ctt_steps[int(step/2)], floaty))\n",
" plt.hist(pull_dic[key][step], bins = 50, range = (-5,5))\n",
" plt.xlabel('Pull')\n",
" plt.savefig('data/CLs/plots/{}/{:.2f}Ctt{}s{}f.png'.format(key, Ctt_steps[int(step/2)], step, floaty))"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"BR: 0.0003\n",
"0.0\n",
"\n",
"BR: 0.0004\n",
"0.0\n",
"\n",
"BR: 0.0005\n",
"0.0\n",
"\n",
"BR: 0.0006\n",
"0.0\n",
"\n"
]
}
],
"source": [
"for s in range(len(l)):\n",
" print('BR: {:.4f}'.format(BR_steps[s+1]))\n",
" print(l[s]/len(np.where(np.array(CLs_values[0]) < np.mean(CLs_values[0]))[0]))\n",
" print()\n",
"# print(l[s], len(CLs_values[s]))\n",
"# print()\n",
"# print(len(CLs_values[0])/2)\n",
"# print(len(np.where(np.array(CLs_values[0]) < np.mean(CLs_values[0]))[0]))\n",
"# plt.clf()\n",
"# # plt.title('Pull {} - Ctt value {:.2f} - floating {}'.format(key, Ctt_steps[int(step/2)], floaty))\n",
"# plt.hist(CLs_values[0], bins = 150, range = (-25, 50), label = 'Ctt fixed to 0')\n",
"# plt.axvline(x=np.mean(CLs_values[0]),color='red', linewidth=1.0, linestyle = 'dotted')\n",
"# plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [],
"source": [
"# for step in range(2*ste):\n",
"# for key in pull_dic.keys():\n",
"# print(pull_dic[key][step])"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [],
"source": [
"# for param in total_f_fit.get_dependents():\n",
"# if param.floating:\n",
"# print(params[param]['value'])"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"39 s\n"
]
}
],
"source": [
"print(display_time(int(time.time()-start)))"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [],
"source": [
"# variab['mi'] =! mi"
]
}
],
"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
}
