{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numpy backup for scaling in the beginning"
   ]
  },
  {
   "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:53: 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 numpy as np\n",
    "from pdg_const import pdg\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import pickle as pkl\n",
    "import sys\n",
    "import time\n",
    "from helperfunctions import display_time, prepare_plot\n",
    "import cmath as c\n",
    "import scipy.integrate as integrate\n",
    "from scipy.optimize import fminbound\n",
    "from array import array as arr\n",
    "import collections\n",
    "from itertools import compress\n",
    "import tensorflow as tf\n",
    "import zfit\n",
    "from zfit import ztf\n",
    "from IPython.display import clear_output\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def formfactor( q2, subscript): #returns real value\n",
    "    #check if subscript is viable\n",
    "\n",
    "    if subscript != \"0\" and subscript != \"+\" and subscript != \"T\":\n",
    "        raise ValueError('Wrong subscript entered, choose either 0, + or T')\n",
    "\n",
    "    #get constants\n",
    "\n",
    "    mK = pdg['Ks_M']\n",
    "    mbstar0 = pdg[\"mbstar0\"]\n",
    "    mbstar = pdg[\"mbstar\"]\n",
    "    b0 = pdg[\"b0\"]\n",
    "    bplus = pdg[\"bplus\"]\n",
    "    bT = pdg[\"bT\"]\n",
    "\n",
    "    mmu = pdg['muon_M']\n",
    "    mb = pdg['bquark_M']\n",
    "    ms = pdg['squark_M']\n",
    "    mB = 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)*(np.sqrt(mB)-np.sqrt(mK))**2\n",
    "\n",
    "    z_oben = np.sqrt(tpos - q2) - np.sqrt(tpos - tzero)\n",
    "    z_unten = np.sqrt(tpos - q2) + np.sqrt(tpos - tzero)\n",
    "    z = z_oben/ z_unten\n",
    "\n",
    "    #calculate f0\n",
    "\n",
    "    if subscript == \"0\":\n",
    "        prefactor = 1/(1 - q2/(mbstar0**2))\n",
    "        _sum = 0\n",
    "\n",
    "        for i in range(N):\n",
    "            _sum += b0[i]*(np.power(z,i))\n",
    "\n",
    "        return 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",
    "            b = bT\n",
    "        else:\n",
    "            b = bplus\n",
    "\n",
    "        for i in range(N):\n",
    "            _sum += b[i] * (np.power(z, i) - ((-1)**(i-N)) * (i/N) * np.power(z, N))\n",
    "\n",
    "        return prefactor * _sum\n",
    "\n",
    "def resonance(q, _mass, width, phase, scale):\n",
    "\n",
    "    q2 = np.power(q, 2)\n",
    "\n",
    "    mmu = pdg['muon_M']\n",
    "\n",
    "    p = 0.5 * np.sqrt(q2 - 4*(mmu**2))\n",
    "\n",
    "    p0 =  0.5 * np.sqrt(_mass**2 - 4*mmu**2)\n",
    "\n",
    "    gamma_j = p/ q2 * _mass * width / p0\n",
    "\n",
    "    #Calculate the resonance\n",
    "\n",
    "    _top = _mass * width\n",
    "\n",
    "    _bottom = np.vectorize(complex)(_mass**2 - q2, -_mass*gamma_j)\n",
    "\n",
    "    com = _top/_bottom\n",
    "\n",
    "    #Rotate by the phase\n",
    "\n",
    "    r = scale*np.abs(com)\n",
    "\n",
    "    _phase = np.angle(com)\n",
    "\n",
    "    _phase += phase\n",
    "\n",
    "    com = r * np.exp(np.vectorize(complex)(0.0 , _phase))\n",
    "\n",
    "    return com\n",
    "\n",
    "\n",
    "def axiv_nonres(q):\n",
    "\n",
    "    GF = pdg['GF']\n",
    "    alpha_ew = pdg['alpha_ew']\n",
    "    Vtb = pdg['Vtb']\n",
    "    Vts = pdg['Vts']\n",
    "    C10eff = pdg['C10eff']\n",
    "\n",
    "    mmu = pdg['muon_M']\n",
    "    mb = pdg['bquark_M']\n",
    "    ms = pdg['squark_M']\n",
    "    mK = pdg['Ks_M']\n",
    "    mB = pdg['Bplus_M']\n",
    "\n",
    "    q2 = np.power(q, 2)\n",
    "\n",
    "    #Some helperfunctions\n",
    "\n",
    "    beta = np.sqrt(np.abs(1. - 4. * mmu**2. / q2))\n",
    "\n",
    "    kabs = np.sqrt(mB**2. +np.power(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. * (np.abs(Vtb*Vts))**2. * kabs * beta / (128. * np.pi**5.)\n",
    "\n",
    "    #left term in bracket\n",
    "\n",
    "    bracket_left = 2./3. * kabs**2. * beta**2. *np.abs(np.vectorize(complex)(C10eff, 0.0)*formfactor(q2, \"+\"))**2.\n",
    "\n",
    "    #middle term in bracket\n",
    "\n",
    "    _top = 4. * mmu**2. * (mB**2. - mK**2.) * (mB**2. - mK**2.)\n",
    "\n",
    "    _under = q2 * mB**2.\n",
    "\n",
    "    bracket_middle = _top/_under *np.power(np.abs(np.vectorize(complex)(C10eff, 0.0) * formfactor(q2, \"0\")), 2)\n",
    "\n",
    "    #Note sqrt(q2) comes from derivation as we use q2 and plot q\n",
    "\n",
    "    return prefactor1 * (bracket_left + bracket_middle) * 2 *np.sqrt(q2)\n",
    "\n",
    "def vec(q, funcs):\n",
    "    \n",
    "    q2 = np.power(q, 2)\n",
    "\n",
    "    GF = pdg['GF']\n",
    "    alpha_ew = pdg['alpha_ew']\n",
    "    Vtb = pdg['Vtb']\n",
    "    Vts = pdg['Vts']\n",
    "    C7eff = pdg['C7eff']\n",
    "\n",
    "    mmu = pdg['muon_M']\n",
    "    mb = pdg['bquark_M']\n",
    "    ms = pdg['squark_M']\n",
    "    mK = pdg['Ks_M']\n",
    "    mB = pdg['Bplus_M']\n",
    "\n",
    "    #Some helperfunctions\n",
    "\n",
    "    beta = np.sqrt(np.abs(1. - 4. * mmu**2. / q2))\n",
    "\n",
    "    kabs = np.sqrt(mB**2. + np.power(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. * (np.abs(Vtb*Vts))**2 * kabs * beta / (128. * np.pi**5.)\n",
    "\n",
    "    #right term in bracket\n",
    "\n",
    "    prefactor2 = kabs**2 * (1. - 1./3. * beta**2)\n",
    "\n",
    "    abs_bracket = np.abs(c9eff(q, funcs) * formfactor(q2, \"+\") + np.vectorize(complex)(2.0 * C7eff * (mb + ms)/(mB + mK), 0.0) * formfactor(q2, \"T\"))**2\n",
    "\n",
    "    bracket_right = prefactor2 * abs_bracket\n",
    "\n",
    "    #Note sqrt(q2) comes from derivation as we use q2 and plot q\n",
    "\n",
    "    return prefactor1 * bracket_right * 2 * np.sqrt(q2)\n",
    "\n",
    "def c9eff(q, funcs):\n",
    "\n",
    "    C9eff_nr = np.vectorize(complex)(pdg['C9eff'], 0.0)\n",
    "\n",
    "    c9 = C9eff_nr\n",
    "\n",
    "    c9 = c9 + funcs\n",
    "\n",
    "    return c9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "jpsi_mass, jpsi_width, jpsi_phase, jpsi_scale = pdg[\"jpsi\"]\n",
    "psi2s_mass, psi2s_width, psi2s_phase, psi2s_scale = pdg[\"psi2s\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def t_f(x):\n",
    "        \n",
    "    def jpsi_res(q):\n",
    "        return resonance(q, _mass = jpsi_mass, scale = jpsi_scale, phase = jpsi_phase, width = jpsi_width)\n",
    "\n",
    "    def psi2s_res(q):\n",
    "        return resonance(q, _mass = psi2s_mass, scale = psi2s_scale, phase = psi2s_phase, width = psi2s_width)\n",
    "\n",
    "    funcs = jpsi_res(x) + psi2s_res(x)\n",
    "\n",
    "    vec_f = vec(x, funcs)\n",
    "\n",
    "    axiv_nr = axiv_nonres(x)\n",
    "\n",
    "    tot = vec_f + axiv_nr\n",
    "    \n",
    "    return tot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Area: 7.126305680132978e-06, rel. err.: 0.1894614654746863%\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\ipykernel_launcher.py:4: IntegrationWarning: The integral is probably divergent, or slowly convergent.\n",
      "  after removing the cwd from sys.path.\n"
     ]
    }
   ],
   "source": [
    "x_min = 2*pdg['muon_M']\n",
    "x_max = (pdg[\"Bplus_M\"]-pdg[\"Ks_M\"]-0.1)\n",
    "\n",
    "result, err = integrate.quad(lambda x: t_f(x), 3600, 3800, limit = 100000000)\n",
    "print(\"Area: {0}, rel. err.: {1}%\".format(result, err/result))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tot:  Area: 0.0013270139056058355, rel. err.: 0.001434977783938156%\n",
      "jpsi: Area: 6.527491388105805e-06, rel. err.: 0.0012105326520570792% (before scaling)\n"
     ]
    }
   ],
   "source": [
    "print(\"tot:  Area: 0.0013270139056058355, rel. err.: 0.001434977783938156%\")\n",
    "print(\"jpsi: Area: 6.527491388105805e-06, rel. err.: 0.0012105326520570792% (before scaling)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEDCAYAAADdpATdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAFvBJREFUeJzt3X+QXeV93/H3V9qVFoT4KeESJCMR5AjVo2JHwYzF2EpMYhkIJBlaxFgTGhgz0xa3tZN2cN3ihPzRieMhradksAY8IW7CDzskURjZmImhaTKAEUFghFAsED+2IpEsQBgL/Vjtt3/cs+Jqtcte2Gd179nzfs3s7DnnPvfscx/pfu5zn/Pc50ZmIklqlhndroAk6dgz/CWpgQx/SWogw1+SGsjwl6QGMvwlqYG6Gv4R8fWI2BkRTxc635cjYnNEbImIr0ZElDivJE033e75/xGwusSJIuKjwEpgOfBB4OeAj5c4tyRNN10N/8z8G+DV9mMR8dMR8Z2IeDwi/m9ELO30dMAAMAuYDfQD/1S0wpI0TXS75z+WdcBnM/Nngd8C/rCTO2Xmw8CDwCvVz/2ZuWXKailJNdbX7Qq0i4gTgI8C32wbrp9d3fZrwE1j3O3/ZeYnI+Ic4FxgQXX8gYj4WPXuQpLUpqfCn9Y7kdcz87zRN2TmvcC973DfXwUeycw3ASLi28AFgOEvSaP01LBPZr4BbI+IfwkQLf+iw7u/BHw8Ivoiop/WxV6HfSRpDBOG/0TTMauA/mpEbIuIpyLiw53+8Yi4E3gY+JmIGIyIa4FPA9dGxJPAZuDyDk/3LeA54AfAk8CTmflXndZFkpokJlrSOSI+BrwJ/HFmfnCM2y8GPgtcDHwE+J+Z+ZEpqKskqZAJe/5jTccc5XJaLwyZmY8AJ0fEGaUqKEkqr8QF3zOBl9v2B6tjr4wuGBHXAdcBzJkz52eXLu10Cr+kXvf63gO8/Npb9M+YwdIz5na7OtPW448//qPMnD/Z85QI/7GWUBhzLCkz19Gax8+KFSty48aNBf68pF5w798P8vl7nuSMkwZ4+Auf6HZ1pq2IeLHEeUrM9hkEFrbtLwB2FDivpBrym2HroUT4rwd+vZr1cwGwJzOPGvKRJPWOCYd9qumYq4B5ETEIfInWujlk5q3ABlozfbYBe4HfmKrKSupd9vjrZcLwz8yrJrg9gX9XrEaS9B4dPHiQwcFB9u3b1+2qTNrAwAALFiygv79/Ss7fa8s7SNJ7Njg4yNy5c1m0aBF1/jqPzGT37t0MDg6yePHiKfkbPbW8gyRNxr59+zjttNNqHfwAEcFpp502pe9gDH9JRfTKkH/dg3/EVD8Ow1+SGsjwl6Rj6KGHHuLSSy8FYP/+/Vx00UWcd9553H333ce0Hl7wlVRU9swAUO974oknOHjwIJs2bTrmf9uev6QiJlohuAleeOEFli5dytVXX83y5cu54oor2Lt3L9/5zndYunQpF154Iffe2/pOqp07d7J27Vo2bdrEeeedx3PPPXdM62rPX1JRMeZyX8fe7/zVZp7Z8UbRcy77qRP50i//83css3XrVm6//XZWrlzJNddcw80338zXvvY1vve973HOOedw5ZVXAnD66adz22238ZWvfIX77ruvaD07Yc9fUlFNH/ZZuHAhK1euBGDt2rVs3LiRxYsXs2TJEiKCtWvXdrmGLfb8JU1LE/XQp8roKZp79uzpyemn9vwlFdHs/v7bXnrpJR5++GEA7rzzTi666CK2b99+eEz/zjvv7Gb1DjP8Jamgc889lzvuuIPly5fz6quv8rnPfY5169ZxySWXcOGFF3LWWWd1u4qAwz6SVNSMGTO49dZbjzi2evVqnn322aPKrlq1ilWrVh2jmh3Jnr8kNZDhL6kMB/1ZtGgRTz/9dLer0RHDX9K0Ml0+bDbVj8PwlzRtDAwMsHv37tq/AIys5z8wMDBlf8MLvpKK6mbuLliwgMHBQXbt2tW9ShQy8k1eU8Xwl1REL3yyt7+/f8q++Wq6cdhHUlE9+GFWjcHwl1RUzYfbG8Pwl6QGMvwlqYEMf0lFONxTL4a/JDWQ4S9JDWT4SyrK0Z96MPwlFWHo14vhL6koP+NVD4a/pKJ8B1APhr8kNZDhL6kI5/nXi+EvSQ1k+EtSAxn+ktRAHYV/RKyOiK0RsS0ibhjj9vdHxIMR8UREPBURF5evqqRe1gtf5qLOTRj+ETETuAX4FLAMuCoilo0q9l+BezLzQ8Aa4A9LV1SSVE4nPf/zgW2Z+XxmHgDuAi4fVSaBE6vtk4Ad5aooSSqtk/A/E3i5bX+wOtbut4G1ETEIbAA+O9aJIuK6iNgYERunwxcsSzqaUz7roZPwH+vT2qP/ea8C/igzFwAXA9+IiKPOnZnrMnNFZq6YP3/+u6+tpJ5l6NdLJ+E/CCxs21/A0cM61wL3AGTmw8AAMK9EBSXVi1/gXg+dhP9jwJKIWBwRs2hd0F0/qsxLwCcAIuJcWuHvuI7UQL4DqIcJwz8zh4DrgfuBLbRm9WyOiJsi4rKq2G8Cn4mIJ4E7gX+d6X8BSepVfZ0UyswNtC7kth+7sW37GWBl2apJqhN7e/XiJ3wlqYEMf0lqIMNfUmEOANWB4S9JDWT4SyrDCX61YvhLUgMZ/pLUQIa/JDWQ4S+pCEf868Xwl6QGMvwlqYEMf0lFOeOzHgx/SUUY+vVi+Esqyi9zqQfDX1JRvgOoB8NfkhrI8JdUhF/eVy+GvyQ1kOEvSQ1k+EtSAxn+kopwxL9eDH9JaiDDX5IayPCXVJTDP/Vg+Esqwmn+9WL4SyrKpX3qwfCXVJRvAOrB8JekBjL8JRVhj79eDH9JaiDDX5IayPCXpAYy/CUV5br+9WD4SyrC0K+XjsI/IlZHxNaI2BYRN4xT5l9FxDMRsTki/rRsNSXVRfgN7rXQN1GBiJgJ3AL8IjAIPBYR6zPzmbYyS4AvACsz87WIOH2qKiypt/kOoB466fmfD2zLzOcz8wBwF3D5qDKfAW7JzNcAMnNn2WpKkkrqJPzPBF5u2x+sjrX7APCBiPi7iHgkIlaPdaKIuC4iNkbExl27dr23GkuSJq2T8B9rAG/0+7o+YAmwCrgKuC0iTj7qTpnrMnNFZq6YP3/+u62rpB42MtrjmH89dBL+g8DCtv0FwI4xyvxlZh7MzO3AVlovBpKkHtRJ+D8GLImIxRExC1gDrB9V5i+AnweIiHm0hoGeL1lRSb0tXd2nViYM/8wcAq4H7ge2APdk5uaIuCkiLquK3Q/sjohngAeB/5SZu6eq0pKkyZlwqidAZm4ANow6dmPbdgKfr34kNZAzPOvFT/hKKmIk+73cWw+GvyQ1kOEvqQiHferF8JdUlNP868Hwl1SEUz3rxfCXVITDPvVi+EtSAxn+ktRAhr+kIt5ex98rvnVg+EtSAxn+korwgm+9GP6SinKefz0Y/pKKsONfL4a/pCIc9qkXw19SUY761IPhL6kIl3eoF8NfkhrI8JdUhGP+9WL4SyrC7K8Xw1+SGsjwl1SG4z61YvhLKiJH/VZvM/wlqYEMf0lFOOpTL4a/pKJ8EagHw19SEX7Ct14Mf0lF2OOvF8NfUmG+CtSB4S+pCCO/Xgx/SWogw19SESNj/o7914PhL6kIZ/vUi+EvSQ1k+EsqI4/4pR5n+EsqwtCvl47CPyJWR8TWiNgWETe8Q7krIiIjYkW5KkqSSpsw/CNiJnAL8ClgGXBVRCwbo9xc4N8Dj5aupKTel9U0n3S6Ty100vM/H9iWmc9n5gHgLuDyMcr9LvBlYF/B+kmSpkAn4X8m8HLb/mB17LCI+BCwMDPve6cTRcR1EbExIjbu2rXrXVdWUu+yw18vnYR/jHHs8D9zRMwA/gD4zYlOlJnrMnNFZq6YP39+57WU1PP8Jq966ST8B4GFbfsLgB1t+3OBDwIPRcQLwAXAei/6SlLv6iT8HwOWRMTiiJgFrAHWj9yYmXsyc15mLsrMRcAjwGWZuXFKaiypJznsUy8Thn9mDgHXA/cDW4B7MnNzRNwUEZdNdQUl1YsvAvXQ10mhzNwAbBh17MZxyq6afLUk1Y1r+9SLn/CVVIQ9/nox/CUVNeyrQC0Y/pLKMvtrwfCXVMTh5R26XA91xvCXVJTDPvVg+Esq4vAnfM3+WjD8JRXllM96MPwlFTHS4x82+2vB8JdURPo9jrVi+Esqygu+9WD4Syoi7fjXiuEvqSi/xrEeDH9JRYxEvhd868Hwl1REe4ff3n/vM/wlFWf29z7DX1IhOcaWepXhL6mI4eG2bbv+Pc/wl1RE+7IOZn/vM/wlFdE+y8f1fXqf4S+piPahHnv+vc/wl1TGEVM9u1cNdcbwl1REe8/fC769z/CXVMSRY/7qdYa/pCLaA99P+PY+w19SEUcO+3SxIuqI4S+piCN6+4Z/zzP8JRXRnv1e8O19hr+kIo6Y59/Feqgzhr+kIoZd0rlWDH9JRRw57NO9eqgzhr+kIvKIYR/Tv9cZ/pKKOHKef9eqoQ4Z/pKKcGG3ejH8JRUx7FTPWuko/CNidURsjYhtEXHDGLd/PiKeiYinIuKvI+Ks8lWV1MvSqZ61MmH4R8RM4BbgU8Ay4KqIWDaq2BPAisxcDnwL+HLpikrqbelUz1rppOd/PrAtM5/PzAPAXcDl7QUy88HM3FvtPgIsKFtNSb3OMf966ST8zwRebtsfrI6N51rg22PdEBHXRcTGiNi4a9euzmspqecd2fPvXj3UmU7CP8Y4NuY/bUSsBVYAvz/W7Zm5LjNXZOaK+fPnd15LST3PL3Opl74OygwCC9v2FwA7RheKiIuALwIfz8z9ZaonqS5c1LNeOun5PwYsiYjFETELWAOsby8QER8CvgZclpk7y1dTUq9r/1SvF3x734Thn5lDwPXA/cAW4J7M3BwRN0XEZVWx3wdOAL4ZEZsiYv04p5M0TQ27tk+tdDLsQ2ZuADaMOnZj2/ZFheslqWaOHOc3/Xudn/CVVISretaL4S+piHSef60Y/pKKOORUz1ox/CUVMXTo7cA/5LhPzzP8JRVxaDiZ3deKlCHDv+cZ/pKKOJTJQP/M1vbwcJdro4kY/pKKaO/5Hzxkz7/XGf6Sihg6lMzub0WKY/69z/CXVESr598a9nHMv/cZ/pKKGGq/4HvIMf9eZ/hLKmK47YKvPf/eZ/hLKmLo0HBbz9/w73WGv6Qijpzn77BPrzP8JRUx1HbB19k+vc/wl1TEcL491dNhn95n+EsqYmg4GXCqZ20Y/pImbXg4yYRZjvnXhuEvadL2D7XC/vjZVc/fYZ+eZ/hLmrT9Q4cAOGFW65thveDb+wx/SZM20vOfM7sV/gf8hG/PM/wlTdr+g62wnzvQV+0f6mZ11AHDX9KkjQz7HDdrJsf1z2TvAcO/1xn+kiZtZNhndt9Mjp81k732/Hue4S9p0vZVYT/QP4PjZs3kLXv+Pc/wlzRp+w4e2fM3/Huf4S9p0t7cfxBoXfA9blafwz41YPhLmrQ33hoC4MTj+jm+fyZvHRjqco00EcNf0qTteavV8z9xoI85s/v48T7Dv9cZ/pIm7Y19B5kRMGdWH/NOmMXunxzodpU0AcNf0qTt/skBTjqunxkzgvlzZ7P7zf0u8dDjDH9Jk/bK62/xUycfB8D8ubMZTnjV3n9PM/wlTdqO1/cdDv/3nThQHXurm1XSBAx/SZNyYGiY7T/6CYvnzQHgZ943F4Bn//GNblZLEzD8JU3K5h17OHBomOULTgLg/acezwmz+9j08p4u10zvxPCXNCn3PfUKfTOClT89D4AZM4KPf2A+3938j4eXfVDv6Sj8I2J1RGyNiG0RccMYt8+OiLur2x+NiEWlKyqp93x/+6v870de5NLlZ3DKnFmHj3/6I+9n908O8N83bHHWT4/qm6hARMwEbgF+ERgEHouI9Zn5TFuxa4HXMvOciFgD/B5w5VRUWM2TOXZ4jHOY8aJm3PO8q3O/u7qMp9T5381jHb/seDe0lmp+62DrZ9/BYfYeGGLXj/fz4u69PPbCq/ztth9x1qnH818uOfeIu370nHlcs3IxX/+77fyff9jFJ859H0tOP4F/dtIAJx7Xz4kDfczum0nfzKBvxgz6ZwZ9M2fQNyOYEQFABAQQEUR13ojWviZvwvAHzge2ZebzABFxF3A50B7+lwO/XW1/C/hfERE53rNtEm7/2+3c/N2tpU97TL2bJ2GpMBjvhp4KG9XK2fPncP3Pn8NnPnY2Jw70H3X7f7v0XFYsOoU/efRFvvHIixwYKvvtXiMvDq3tqF4o4PBLRdvtnZxrwjIdnq2Tc33pl5dx5c+9v6PzTZVOwv9M4OW2/UHgI+OVycyhiNgDnAb8qL1QRFwHXFftvhkRu0eXabB52BZgO7Tr6bZ4EXgQ+K1j8+d6ui3erTW/C2ve213nAWeVqEMn4T/W69jovlsnZcjMdcC6w3eK2JiZKzqow7RnW7TYDm+zLd5mW7RU7bCoxLk6ueA7CCxs218A7BivTET0AScBr5aooCSpvE7C/zFgSUQsjohZtN6trB9VZj1wdbV9BfC9qRjvlySVMeGwTzWGfz1wPzAT+Hpmbo6Im4CNmbkeuB34RkRso9Xj73Q4a93ERRrDtmixHd5mW7zNtmgp1g5hB12SmsdP+EpSAxn+ktRAXQv/iZaMqLuI+HpE7IyIp9uOnRoRD0TED6vfp1THIyK+WrXFUxHx4bb7XF2V/2FEXD3W3+p1EbEwIh6MiC0RsTki/kN1vFHtEREDEfH9iHiyaoffqY4vrpZF+WG1TMqs6vi4y6ZExBeq41sj4pPdeUSTFxEzI+KJiLiv2m9cW0TECxHxg4jYFBEbq2NT/9zIzGP+Q+vC8XPA2cAs4ElgWTfqMoWP8WPAh4Gn2459Gbih2r4B+L1q+2Lg27Q+L3EB8Gh1/FTg+er3KdX2Kd1+bO+hLc4APlxtzwX+AVjWtPaoHs8J1XY/8Gj1+O4B1lTHbwX+TbX9b4Fbq+01wN3V9rLqOTMbWFw9l2Z2+/G9xzb5PPCnwH3VfuPaAngBmDfq2JQ/N7rV8z+8ZERmHgBGloyYNjLzbzj6sw6XA3dU23cAv9J2/I+z5RHg5Ig4A/gk8EBmvpqZrwEPAKunvvZlZeYrmfn31faPgS20PhXeqPaoHs+b1W5/9ZPAL9BaFgWOboeR9vkW8ImIiOr4XZm5PzO3A9toPadqJSIWAJcAt1X7QUPbYgxT/tzoVviPtWTEmV2qy7H0vsx8BVqBCJxeHR+vPaZdO1Vv1z9Eq9fbuPaohjk2ATtpPUGfA17PzKGqSPtjOmLZFGBk2ZTat0PlfwD/GRhZ9Oc0mtkWCXw3Ih6P1hI4cAyeG50s7zAVOloOokHGa49p1U4RcQLwZ8B/zMw3YvwVsKZte2TmIeC8iDgZ+HPg3LGKVb+nbTtExKXAzsx8PCJWjRweo+i0bwtgZWbuiIjTgQci4tl3KFusHbrV8+9kyYjp6J+qt2hUv3dWx8drj2nTThHRTyv4/yQz760ON7Y9MvN14CFa47YnR2tZFDjyMY23bMp0aIeVwGUR8QKtYd9foPVOoHFtkZk7qt87aXUIzucYPDe6Ff6dLBkxHbUvg3E18Jdtx3+9upJ/AbCneqt3P/BLEXFKdbX/l6pjtVKNzd4ObMnMm9tualR7RMT8qsdPRBwHXETr+seDtJZFgaPbYaxlU9YDa6oZMIuBJcD3j82jKCMzv5CZC7K1SNkaWo/t0zSsLSJiTkTMHdmm9X/6aY7Fc6OLV7gvpjXr4zngi92qxxQ+vjuBV4CDtF6Vr6U1RvnXwA+r36dWZYPWF+Y8B/wAWNF2nmtoXcTaBvxGtx/Xe2yLC2m9BX0K2FT9XNy09gCWA09U7fA0cGN1/GxagbUN+CYwuzo+UO1vq24/u+1cX6zaZyvwqW4/tkm2yyrenu3TqLaoHu+T1c/mkSw8Fs8Nl3eQpAbyE76S1ECGvyQ1kOEvSQ1k+EtSAxn+ktRAhr8kNZDhL0kN9P8Bmz5lnu9Pn6cAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "test_q = np.linspace(x_min, x_max, 2000000)\n",
    "calcs_test = t_f(test_q)\n",
    "plt.clf()\n",
    "# plt.plot(x_part, calcs, '.')\n",
    "plt.plot(test_q, calcs_test, label = 'pdf')\n",
    "# plt.plot(test_q, res_y, label = 'res')\n",
    "plt.legend()\n",
    "plt.ylim(0.0, 1e-8)\n",
    "# plt.yscale('log')\n",
    "# plt.xlim(3080, 3110)\n",
    "plt.savefig('np-test.png')\n",
    "# print(jpsi_width)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}