diff --git a/.ipynb_checkpoints/raremodel-nb-checkpoint.ipynb b/.ipynb_checkpoints/raremodel-nb-checkpoint.ipynb index b3e49f8..4688758 100644 --- a/.ipynb_checkpoints/raremodel-nb-checkpoint.ipynb +++ b/.ipynb_checkpoints/raremodel-nb-checkpoint.ipynb @@ -781,200 +781,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Tensorflow scaling" - ] - }, - { - "cell_type": "code", - "execution_count": 16, - "metadata": {}, - "outputs": [], - "source": [ - "# def scaling_func(x):\n", - "\n", - "# funcs = resonance(x, _mass = ztf.constant(jpsi_mass), scale = ztf.constant(jpsi_scale), phase = ztf.constant(jpsi_phase), width = ztf.constant(jpsi_width)) + resonance(x, _mass = ztf.constant(psi2s_mass), scale = ztf.constant(psi2s_scale), phase = ztf.constant(psi2s_phase), width = ztf.constant(psi2s_width))\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\n", - "\n", - "\n", - "# def s_func(x):\n", - " \n", - "# q = ztf.constant(x)\n", - " \n", - "# return zfit.run(scaling_func(q))\n", - " \n", - "\n", - "# print(integrate.quad(s_func, x_min, x_max, limit = 50))" - ] - }, - { - "cell_type": "code", - "execution_count": 17, - "metadata": {}, - "outputs": [], - "source": [ - "# factor_jpsi = pdg[\"NR_auc\"]*pdg[\"jpsi_BR\"]/(pdg[\"NR_BR\"]*pdg[\"jpsi_auc\"])\n", - "# factor_jpsi = pdg[\"NR_auc\"]*pdg[\"jpsi_BR\"]/(pdg[\"NR_BR\"]*inte_fl)\n", - "# print(np.sqrt(factor_jpsi)*jpsi_scale)\n", - "# print(np.sqrt(factor_jpsi))\n", - "# # print(psi2s_scale)\n", - "# factor_psi2s = pdg[\"NR_auc\"]*pdg[\"psi2s_BR\"]/(pdg[\"NR_BR\"]*pdg[\"psi2s_auc\"])\n", - "# factor_psi2s = pdg[\"NR_auc\"]*pdg[\"psi2s_BR\"]/(pdg[\"NR_BR\"]*inte_fl)\n", - "# print(np.sqrt(factor_psi2s)*psi2s_scale)\n", - "# print(np.sqrt(factor_psi2s))" - ] - }, - { - "cell_type": "code", - "execution_count": 18, - "metadata": {}, - "outputs": [], - "source": [ - "# def _t_f(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", - "# funcs = psi2s_res(xq) + jpsi_res(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 t_f(x):\n", - "# _ = np.array(x)\n", - "# probs = zfit.run(_t_f(_))\n", - "# return probs" - ] - }, - { - "cell_type": "code", - "execution_count": 19, - "metadata": {}, - "outputs": [], - "source": [ - "# print(36000*(1+ pdg[\"jpsi_BR\"]/pdg[\"NR_BR\"] + pdg[\"psi2s_BR\"]/pdg[\"NR_BR\"]))" - ] - }, - { - "cell_type": "code", - "execution_count": 20, - "metadata": {}, - "outputs": [], - "source": [ - "# start = time.time()\n", - "\n", - "# result, err = integrate.quad(lambda x: t_f(x), x_min, x_max, limit = 5)\n", - "# print(result, \"{0:.2f} %\".format(err/result))\n", - "# print(\"Time:\", time.time()-start)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ "# Sampling\n", - "## One sample\n", - "! total_f.sample() always returns the same set !" - ] - }, - { - "cell_type": "code", - "execution_count": 21, - "metadata": {}, - "outputs": [], - "source": [ - "# nevents = int(pdg[\"number_of_decays\"])\n", - "# event_stack = 5000\n", - "\n", - "# calls = int(nevents/event_stack + 1)\n", - "\n", - "# total_samp = []\n", - "\n", - "# start = time.time()\n", - "\n", - "# samp = total_f.sample(n=event_stack)\n", - "# s = samp.unstack_x()\n", - "\n", - "# for call in range(calls):\n", - "\n", - "# sam = zfit.run(s)\n", - "# clear_output(wait=True)\n", - " \n", - "# # if call != 0:\n", - "# # print(np.sum(_last_sam-sam))\n", - " \n", - "# # _last_sam = sam\n", - " \n", - "# c = call + 1 \n", - "# print(\"{0}/{1}\".format(c, calls))\n", - "# print(\"Time taken: {}\".format(display_time(int(time.time() - start))))\n", - "# print(\"Projected time left: {}\".format(display_time(int((time.time() - start)/c*(calls-c)))))\n", - " \n", - "# with open(\"data/zfit_toys/toy_1/{}.pkl\".format(call), \"wb\") as f:\n", - "# pkl.dump(sam, f, pkl.HIGHEST_PROTOCOL)" - ] - }, - { - "cell_type": "code", - "execution_count": 22, - "metadata": {}, - "outputs": [], - "source": [ - "# print(\"Time to generate full toy: {} s\".format(int(time.time()-start)))\n", - "\n", - "# total_samp = []\n", - "\n", - "# for call in range(calls):\n", - "# with open(r\"data/zfit_toys/toy_1/{}.pkl\".format(call), \"rb\") as input_file:\n", - "# sam = pkl.load(input_file)\n", - "# total_samp = np.append(total_samp, sam)\n", - "\n", - "# total_samp = total_samp.astype('float64')\n", - "\n", - "# data2 = zfit.data.Data.from_numpy(array=total_samp[:int(nevents)], obs=obs)\n", - "\n", - "# print(total_samp[:nevents].shape)" - ] - }, - { - "cell_type": "code", - "execution_count": 23, - "metadata": {}, - "outputs": [], - "source": [ - "# bins = int((x_max-x_min)/7)\n", - "\n", - "# # calcs = zfit.run(total_test_tf(samp))\n", - "\n", - "# plt.hist(total_samp[:event_stack], bins = bins, range = (x_min,x_max))\n", - "\n", - "# # plt.plot(sam, calcs, '.')\n", - "# # plt.plot(test_q, calcs_test)\n", - "# plt.ylim(0, 20)\n", - "# # plt.xlim(3000, 3750)\n", - "\n", - "# plt.savefig('test2.png')\n", - "# 1-(0.21+0.62)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ "## Toys" ] }, @@ -1091,7 +898,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 28, "metadata": { "scrolled": false }, @@ -1100,9 +907,9 @@ "name": "stdout", "output_type": "stream", "text": [ - "3/6 of Toy 1/1\n", - "Time taken: 3 min, 27 s\n", - "Projected time left: 3 min, 27 s\n" + "6/6 of Toy 1/1\n", + "Time taken: 1 min, 22 s\n", + "Projected time left: \n" ] } ], @@ -1149,7 +956,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 29, "metadata": {}, "outputs": [], "source": [ @@ -1166,9 +973,18 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 30, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Time to generate full toy: 82 s\n", + "(5404696,)\n" + ] + } + ], "source": [ "print(\"Time to generate full toy: {} s\".format(int(time.time()-start)))\n", "\n", @@ -1190,9 +1006,29 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 31, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(5404696,)\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAD8CAYAAACVZ8iyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAXBklEQVR4nO3df7BcZZ3n8fd3IZoVMkOIFzZFGG6gApIZJLBXpGShVFB+WUF0Zgw1pdGhyMpg1TiO7kaxZtmyqGJ01V1qdrHCQgFVCDj8iHFgdkyxOjilIDcSQjBGEiYzXkkld8KAoRCXH9/9o8+FJvS9t2/36e7T3e9XVVef8/TTp5/z3NvPp8+PPh2ZiSRpuP2bXjdAktR7hoEkyTCQJBkGkiQMA0kShoEkiSbCICKOjojvRcS2iHg8Iv60KD88IjZGxBPF/cKiPCLi2ojYERFbIuLUTq+EJKk9zWwZvAT8eWaeCJwOXBERy4G1wP2ZuQy4v5gHOB9YVtzWANeV3mpJUqlmDYPM3J2ZPymm9wPbgKOAi4Cbi2o3Ax8spi8CbsmaB4HDImJx6S2XJJXm4LlUjohR4BTgIeDIzNwNtcCIiCOKakcBv6h72kRRtvuAZa2htuXAIYcc8u/f9ra3tdB8abg89stnZ61z0lG/3YWWTK9RG3vZpsd++WzP+6RTNm3a9C+ZOVLGspoOg4g4FLgL+HRm/ioipq3aoOwN17zIzHXAOoCxsbEcHx9vtinS0Bpde++sdcavubALLZleozb2sk2ja+/teZ90SkT8U1nLaupsooiYRy0Ibs3Mu4viPVO7f4r7vUX5BHB03dOXAE+V01xJUic0czZRADcA2zLza3UPbQBWF9OrgW/XlX+sOKvodODZqd1JkqRqamY30RnAR4HHImJzUfYF4BrgWxFxKfDPwB8Uj90HXADsAJ4HPlFqiyVJpZs1DDLzH2h8HADg7Ab1E7iizXbx4osvMjExwQsvvNDuoipt/vz5LFmyhHnz5vW6KZKG2JzOJuqmiYkJFixYwOjoKDMcrO5rmcm+ffuYmJhg6dKlvW6OpCFW2ctRvPDCCyxatGhggwAgIli0aNHAb/1Iqr7KhgEw0EEwZRjWUVL1VToMJEndUdljBgdq5ss2c7GrhS+hXHXVVRx66KF89rOfbfj4+vXrOf7441m+fHm7zZOkrnLLoETr16/npz/9aa+bIUlzZhjM4uqrr+aEE07gnHPOYfv27QBcf/31vOMd7+Dkk0/mwx/+MM8//zw//OEP2bBhA5/73OdYsWIFO3fubFhPkqrIMJjBpk2buP3223nkkUe4++67efjhhwH40Ic+xMMPP8yjjz7KiSeeyA033MC73vUuVq5cyVe+8hU2b97Mcccd17CeJFVR3xwz6IUf/OAHXHzxxbzlLW8BYOXKlQBs3bqVL37xizzzzDM899xznHvuuQ2f32w9Seo1w2AWjU79/PjHP8769es5+eSTuemmm/j+97/f8LnN1pOkXnM30QzOOuss7rnnHn7961+zf/9+vvOd7wCwf/9+Fi9ezIsvvsitt976av0FCxawf//+V+enqydJVdM3WwatnArarlNPPZWPfOQjrFixgmOOOYYzzzwTgC996Uu8853v5JhjjuGkk056NQBWrVrFZZddxrXXXsudd945bT1Jqpq+CYNeufLKK7nyyivfUH755Ze/oeyMM8543amll19+ecN6klQ17iaSJBkGkqSKh0HtpxEG2zCso6Tqq2wYzJ8/n3379g30YDn1ewbz58/vdVMkDbnKHkBesmQJExMTTE5O9ropHTX1S2eS1EuzhkFE3Ah8ANibmb9XlN0BnFBUOQx4JjNXRMQosA3YXjz2YGZ+spWGzZs3z1//kgbE6Np7e3J6uJrXzJbBTcBfAbdMFWTmR6amI+KrwLN19Xdm5oqyGihJ6rxZwyAzHyg+8b9B1K7V8IfAe8ttliSpm9o9gHwmsCczn6grWxoRj0TE30fEmW0uX5LUBe2GwSXAbXXzu4HfycxTgM8A34yI32r0xIhYExHjETE+6AeJJfVW2b+UOIhaDoOIOBj4EHDHVFlm/iYz9xXTm4CdwPGNnp+Z6zJzLDPHRkZGWm2GJKkE7WwZnAP8LDMnpgoiYiQiDiqmjwWWAU+210RJUqfNGgYRcRvwI+CEiJiIiEuLh1bx+l1EAGcBWyLiUeBO4JOZ+XSZDZYkla+Zs4kumab84w3K7gLuar9ZkqRuquzlKCRJ3WMYSJIMA6kfjK69t+nTIz2NUq0wDCRJhoEkyTCQJGEYSJIwDCRJGAaSJAwDSRKGgSQJw0CShGEgqQf8lnT1GAaSJMNAkmQYSJIwDCRJGAaSJAwDSRJNhEFE3BgReyNia13ZVRHxy4jYXNwuqHvs8xGxIyK2R8S5nWq4JKk8zWwZ3ASc16D865m5orjdBxARy4FVwO8Wz/lfEXFQWY2VJHXGrGGQmQ8ATze5vIuA2zPzN5n5j8AO4LQ22idJ6oJ2jhl8KiK2FLuRFhZlRwG/qKszUZS9QUSsiYjxiBifnJxsoxmSpHa1GgbXAccBK4DdwFeL8mhQNxstIDPXZeZYZo6NjIy02AxJUhlaCoPM3JOZL2fmK8D1vLYraAI4uq7qEuCp9pooSeq0lsIgIhbXzV4MTJ1ptAFYFRFvjoilwDLgx+01UZLUaQfPViEibgPeDbw1IiaA/wK8OyJWUNsFtAv4jwCZ+XhEfAv4KfAScEVmvtyZpkuqGq9G2r9mDYPMvKRB8Q0z1L8auLqdRklSGQyn5vkNZEl9bboB3yCYG8NAkmQYSJIMA0kShoEkCcNAkoRhIGlIeHbRzAwDSV3loFxNhoEkyTCQNDzcKpmeYSBJMgwkSYaBJAnDQJKEYSBJwjCQNGQ8o6gxw0CSZBhIVecnWXXDrGEQETdGxN6I2FpX9pWI+FlEbImIeyLisKJ8NCJ+HRGbi9s3Otl4SVI5mtkyuAk474CyjcDvZebbgZ8Dn697bGdmrihunyynmZKkTpo1DDLzAeDpA8q+m5kvFbMPAks60DZJUpeUcczgj4G/rZtfGhGPRMTfR8SZ0z0pItZExHhEjE9OTpbQDGnweLxgZvZPedoKg4i4EngJuLUo2g38TmaeAnwG+GZE/Faj52bmuswcy8yxkZGRdpohSWpTy2EQEauBDwB/lJkJkJm/ycx9xfQmYCdwfBkNlSR1TkthEBHnAf8ZWJmZz9eVj0TEQcX0scAy4MkyGipJ6pxmTi29DfgRcEJETETEpcBfAQuAjQecQnoWsCUiHgXuBD6ZmU83XLAk4X7/qjh4tgqZeUmD4humqXsXcFe7jZI0+EbX3suuay7sdTNU8BvIkoaOWyNvZBhIkgwDSZJhIFVWO7syhm03SP36Dtu6l8UwkCQZBpJ6z0/zvWcYSCqFA3p/MwwkSYaBJMkwkFQh7mrqHcNAkmQYSJIMA2lg9dsul35r76AxDKQB5gDbvmHpQ8NAkmQYSKq2uXwyH5ZP8Z1gGEgaOIbC3BkGFeU/swZNp/+nvcpre5oKg4i4MSL2RsTWurLDI2JjRDxR3C8syiMiro2IHRGxJSJO7VTjJbWuCgNgFdqgmma3DG4CzjugbC1wf2YuA+4v5gHOB5YVtzXAde03U5LKZRC9XlNhkJkPAE8fUHwRcHMxfTPwwbryW7LmQeCwiFhcRmMlVU/Zg+rU8hysu6udYwZHZuZugOL+iKL8KOAXdfUmirLXiYg1ETEeEeOTk5NtNEPSsDAgOqcTB5CjQVm+oSBzXWaOZebYyMhIB5oh9Yd+HOD6sc2NjK69t+l1GZR1nk47YbBnavdPcb+3KJ8Ajq6rtwR4qo3XkYbKXAaoVpffK2W+9qAPzt3WThhsAFYX06uBb9eVf6w4q+h04Nmp3UmSesdPwJpJs6eW3gb8CDghIiYi4lLgGuB9EfEE8L5iHuA+4ElgB3A98Celt1pS6VoJgV4Fh4FVvoObqZSZl0zz0NkN6iZwRTuNkiR1l99AlnrET7fTs2+6zzCQ1BcMiM4yDKSKqOJgZ5uGh2Eg9UC3B7R2X2/YB+BhWH/DQOqCqp7bX/VBrurtGySGgaSmOTgPbh8YBlIFDOoA02+avUjeIP69DANJorkBfhBDYIphIGlG0w2AgzwwDiPDQFLlGTydZxhImpNhGZiHZT2nGAaS1KJBCgzDQJJkGEh6zSB90m3VsPaBYSCpbcM6gA4Sw0CSg7kMA2kYOfjrQIaBJKm5n71sJCJOAO6oKzoW+AvgMOAyYLIo/0Jm3tdyCyV1VLNbCW5NDLaWwyAztwMrACLiIOCXwD3AJ4CvZ+Z/K6WFkqSOazkMDnA2sDMz/ykiSlqkNHyq8um7Ku1Q95R1zGAVcFvd/KciYktE3BgRCxs9ISLWRMR4RIxPTk42qiL1naoPolVvn3qn7TCIiDcBK4G/LoquA46jtgtpN/DVRs/LzHWZOZaZYyMjI+02Q6qcAwdeB+LB0uxvH/SLMrYMzgd+kpl7ADJzT2a+nJmvANcDp5XwGpKkDiojDC6hbhdRRCyue+xiYGsJryFJ6qC2wiAi3gK8D7i7rvjLEfFYRGwB3gP8WTuvIZWhm5vyc9k9NCi7GNT/2gqDzHw+Mxdl5rN1ZR/NzJMy8+2ZuTIzd7ffTKm6HNB1oH78n/AbyFKX9eNAoZkNwt/UMJAkGQaSJMNA6qhB2H2g4WAYaCDNdRAeXXtvS8/RYBumv7FhIDVpmAYGDR/DQJojQ0GDyDCQJBkGklSGft9iNAykktQPBv0+MGj4GAaSJMNAg6edT+V+otewMgykBmb74RJ/uEaDpqzfQJZ6ruwBuZnl+UU1DQq3DKQmOIhr0BkGGnoO9JJhoCFmCKgT+vX/yjDQQGj1Ddivb1xVUz//P7V9ADkidgH7gZeBlzJzLCIOB+4ARoFdwB9m5r+2+1oaXqNr72XXNReWtixJr1fWlsF7MnNFZo4V82uB+zNzGXB/MS+9joOyVB2d2k10EXBzMX0z8MEOvY4kqQRlhEEC342ITRGxpig7MjN3AxT3Rxz4pIhYExHjETE+OTlZQjM06NrdknBLRJpeGV86OyMzn4qII4CNEfGzZp6UmeuAdQBjY2NZQjskSS1qe8sgM58q7vcC9wCnAXsiYjFAcb+33deRWuUWhTS7tsIgIg6JiAVT08D7ga3ABmB1UW018O12Xkf9r5XfGJ5pWfX3ktrX7m6iI4F7ImJqWd/MzP8TEQ8D34qIS4F/Bv6gzdeRKsEA0qBqa8sgM5/MzJOL2+9m5tVF+b7MPDszlxX3T5fTXKl5DtzqtX76H/QbyOor/fTmkvqJYaCWdWJgbnQ8YLbXmelxf3dAvdJv/2uGgbpqLgO3pO4xDNRRDvBSfzAMVIpBGvQHaV2kZhkGFeRgJKnbDANVggEo9ZZhoJaU/W3iuT5ueEjlMgwkSYaByjOXawb5yV6qFsNATSnzQnONlt2L50rd0C//o4aBGpppX73f6pUGj2EwZDo9cJexi8hLVGvQ9MP/smEwhOZy3R9Jw8Ew6BAHWftA6ieGgWY9l99BXRp8hoHewMFfGj6GgSSp9d9AjoijgVuAfwe8AqzLzP8REVcBlwGTRdUvZOZ97Ta0n/TDJ+tutrEf+kMadu1sGbwE/HlmngicDlwREcuLx76emSuK28AFQT8MbvXfB+iH9krqrZa3DDJzN7C7mN4fEduAo8pqmCSpe0o5ZhARo8ApwENF0aciYktE3BgRC8t4jWHRyqf4mbYAerVV4NaI1F/aDoOIOBS4C/h0Zv4KuA44DlhBbcvhq9M8b01EjEfE+OTkZKMqqoB2foxeUv9oKwwiYh61ILg1M+8GyMw9mflyZr4CXA+c1ui5mbkuM8cyc2xkZKSdZlRKK4Nn2QNqq1sXnbwYnaRqazkMIiKAG4Btmfm1uvLFddUuBra23rz+5sAqqV+0s2VwBvBR4L0Rsbm4XQB8OSIei4gtwHuAPyujoVVRlQG+Ku2QNBjaOZvoH4Bo8NBAnkq665oLm6rXyrIlqdeG9hvIVRiEywiPZr9LUIX1lVRdQxsG3eSBWUlVZxj0wGxXCZ2tXivLlqSZGAYzmGlgbXbQ7cTg7IAvqWyGQZ1Wd+d0MhjKvL6QISL1VpXfg4ZBRbhLSFIvGQaSpOEIg3Y+QfvpW9IwGIowmE67A71BIWlQDHUYzKTTF5QzSKThVNXvHQ1VGDR7fn8rdTr5x/UXy6T+V/X378CGwUwDf/1jjQbaufzRqv4HlqRmDGwYgAO1JDVroMOgnsEgSdMbqDCo0m8AS9JMqjY2DVQYSJJa0/KP21RV1dJWkvqBWwaSJMNAktTBMIiI8yJie0TsiIi1nXodcNeQJLWrI2EQEQcB/xM4H1gOXBIRyzvxWpKk9nVqy+A0YEdmPpmZ/w+4HbioQ68FuHUgSe3o1NlERwG/qJufAN5ZXyEi1gBritnnImIf8C8dak+/eSv2xRT7osZ+eM3A9EX8ZVtPfytwTDkt6VwYRIOyfN1M5jpg3atPiBjPzLEOtaev2BevsS9q7IfX2Bc1RT+MlrW8Tu0mmgCOrptfAjzVodeSJLWpU2HwMLAsIpZGxJuAVcCGDr2WJKlNHdlNlJkvRcSngL8DDgJuzMzHZ3naulkeHyb2xWvsixr74TX2RU2p/RCZOXstSdJA8xvIkiTDQJJUkTDo5qUreiUiboyIvRGxta7s8IjYGBFPFPcLi/KIiGuL/tgSEafWPWd1Uf+JiFjdi3VpR0QcHRHfi4htEfF4RPxpUT5UfRER8yPixxHxaNEP/7UoXxoRDxXrdEdxAgYR8eZifkfx+Gjdsj5flG+PiHN7s0bti4iDIuKRiPibYn4o+yIidkXEYxGxOSLGi7LOvz8ys6c3ageYdwLHAm8CHgWW97pdHVjPs4BTga11ZV8G1hbTa4G/LKYvAP6W2vc1TgceKsoPB54s7hcW0wt7vW5z7IfFwKnF9ALg59QuWTJUfVGsz6HF9DzgoWL9vgWsKsq/AVxeTP8J8I1iehVwRzG9vHjPvBlYWryXDur1+rXYJ58Bvgn8TTE/lH0B7ALeekBZx98fVdgy6PqlK3ohMx8Anj6g+CLg5mL6ZuCDdeW3ZM2DwGERsRg4F9iYmU9n5r8CG4HzOt/68mTm7sz8STG9H9hG7RvrQ9UXxfo8V8zOK24JvBe4syg/sB+m+udO4OyIiKL89sz8TWb+I7CD2nuqr0TEEuBC4H8X88GQ9sU0Ov7+qEIYNLp0xVE9aku3HZmZu6E2SAJHFOXT9clA9VWxeX8KtU/FQ9cXxW6RzcBeam/WncAzmflSUaV+nV5d3+LxZ4FFDEA/FP478J+AV4r5RQxvXyTw3YjYFLXL9kAX3h9V+KWzWS9dMYSm65OB6auIOBS4C/h0Zv6q9sGucdUGZQPRF5n5MrAiIg4D7gFObFStuB/YfoiIDwB7M3NTRLx7qrhB1YHvi8IZmflURBwBbIyIn81Qt7S+qMKWwTBfumJPsUlHcb+3KJ+uTwairyJiHrUguDUz7y6Kh7IvADLzGeD71Pb5HhYRUx/S6tfp1fUtHv9tarsdB6EfzgBWRsQuaruJ30ttS2EY+4LMfKq430vtQ8JpdOH9UYUwGOZLV2wApo7yrwa+XVf+seJMgdOBZ4tNw78D3h8RC4uzCd5flPWNYt/uDcC2zPxa3UND1RcRMVJsERAR/xY4h9rxk+8Bv19UO7Afpvrn94H/m7UjhRuAVcUZNkuBZcCPu7MW5cjMz2fmkqxddG0VtXX7I4awLyLikIhYMDVN7f96K914f/T6yHndEfGfU9tnemWv29OhdbwN2A28SC21L6W2n/N+4Ini/vCiblD7caCdwGPAWN1y/pjagbEdwCd6vV4t9MN/oLa5ugXYXNwuGLa+AN4OPFL0w1bgL4ryY6kNYDuAvwbeXJTPL+Z3FI8fW7esK4v+2Q6c3+t1a7Nf3s1rZxMNXV8U6/xocXt8ajzsxvvDy1FIkiqxm0iS1GOGgSTJMJAkGQaSJAwDSRKGgSQJw0CSBPx/Vnu3Jq1jkDcAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], "source": [ "plt.clf()\n", "\n", @@ -1217,7 +1053,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 32, "metadata": {}, "outputs": [], "source": [ @@ -1240,7 +1076,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 33, "metadata": {}, "outputs": [], "source": [ @@ -1249,7 +1085,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 34, "metadata": {}, "outputs": [], "source": [ @@ -1265,17 +1101,73 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 35, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.6628998341322463\n", + "6.258685216704922\n", + "2.249780125960397\n", + "-4.790281653148225\n", + "-0.546187997552245\n", + "-4.053408742543857\n", + "3.3921906479017068\n", + "1.380997849250071\n", + "0.5628459148842957\n", + "------------------------------------------------------------------\n", + "| FCN = -7.175E+05 | Ncalls=253 (253 total) |\n", + "| EDM = 2.21E-05 (Goal: 5E-06) | up = 0.5 |\n", + "------------------------------------------------------------------\n", + "| Valid Min. | Valid Param. | Above EDM | Reached call limit |\n", + "------------------------------------------------------------------\n", + "| True | True | False | False |\n", + "------------------------------------------------------------------\n", + "| Hesse failed | Has cov. | Accurate | Pos. def. | Forced |\n", + "------------------------------------------------------------------\n", + "| False | True | True | True | False |\n", + "------------------------------------------------------------------\n", + "Function minimum: -717521.3025826928\n", + "---------------------------------------------------------------------------------------------\n", + "| | Name | Value | Hesse Err | Minos Err- | Minos Err+ | Limit- | Limit+ | Fixed |\n", + "---------------------------------------------------------------------------------------------\n", + "| 0 | jpsi_p | 1.486 | 0.016 | | |-6.28319 | 6.28319 | |\n", + "| 1 | p3770_p | 2.55 | 0.09 | | |-6.28319 | 6.28319 | |\n", + "| 2 | psi2s_p | 1.344 | 0.027 | | |-6.28319 | 6.28319 | |\n", + "| 3 | omega_p | -5.80 | 0.22 | | |-6.28319 | 6.28319 | |\n", + "| 4 | p4040_p | -3.37 | 0.16 | | |-6.28319 | 6.28319 | |\n", + "| 5 | rho_p | -0.19 | 0.30 | | |-6.28319 | 6.28319 | |\n", + "| 6 | phi_p | -5.59 | 0.15 | | |-6.28319 | 6.28319 | |\n", + "| 7 | p4160_p | 4.04 | 0.08 | | |-6.28319 | 6.28319 | |\n", + "| 8 | p4415_p | -3.13 | 0.13 | | |-6.28319 | 6.28319 | |\n", + "---------------------------------------------------------------------------------------------\n", + "-------------------------------------------------------------------------------------\n", + "| | jpsi_p p3770_p psi2s_p omega_p p4040_p rho_p phi_p p4160_p p4415_p |\n", + "-------------------------------------------------------------------------------------\n", + "| jpsi_p | 1.000 -0.133 0.219 0.015 -0.176 -0.110 0.023 -0.065 -0.143 |\n", + "| p3770_p | -0.133 1.000 -0.527 -0.007 0.014 0.005 -0.011 0.001 0.014 |\n", + "| psi2s_p | 0.219 -0.527 1.000 0.016 -0.312 -0.060 0.025 -0.006 -0.140 |\n", + "| omega_p | 0.015 -0.007 0.016 1.000 -0.007 -0.255 0.034 -0.011 -0.012 |\n", + "| p4040_p | -0.176 0.014 -0.312 -0.007 1.000 0.001 -0.011 -0.381 -0.104 |\n", + "| rho_p | -0.110 0.005 -0.060 -0.255 0.001 1.000 -0.246 0.026 0.017 |\n", + "| phi_p | 0.023 -0.011 0.025 0.034 -0.011 -0.246 1.000 -0.018 -0.018 |\n", + "| p4160_p | -0.065 0.001 -0.006 -0.011 -0.381 0.026 -0.018 1.000 0.164 |\n", + "| p4415_p | -0.143 0.014 -0.140 -0.012 -0.104 0.017 -0.018 0.164 1.000 |\n", + "-------------------------------------------------------------------------------------\n", + "Hesse errors: OrderedDict([(, {'error': 0.015523017204917444}), (, {'error': 0.0949040578205631}), (, {'error': 0.027368429790692872}), (, {'error': 0.21984151021717713}), (, {'error': 0.1583791551448439}), (, {'error': 0.2962129131334841}), (, {'error': 0.15237884170408078}), (, {'error': 0.08086187469444539}), (, {'error': 0.12776706213127875})])\n" + ] + } + ], "source": [ "start = time.time()\n", "\n", "for param in total_f.get_dependents():\n", " param.randomize()\n", " \n", - "for param in total_f.get_dependents():\n", - " print(zfit.run(param))\n", + "# for param in total_f.get_dependents():\n", + "# print(zfit.run(param))\n", " \n", "nll = zfit.loss.UnbinnedNLL(model=total_f, data=data2, fit_range = (x_min, x_max))\n", "\n", @@ -1295,9 +1187,17 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 36, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Time taken for fitting: 2 min, 6 s\n" + ] + } + ], "source": [ "print(\"Time taken for fitting: {}\".format(display_time(int(time.time()-start))))\n", "\n", @@ -1309,16 +1209,29 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 60, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ8AAAD8CAYAAACo9anUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO29eXgc1ZX3/znd2mxJlmTZsmXLtrzINgYbYwtjtsQBAgYTTBJITMKESUh4kwm/rPML8GSbZMLMMG8Sz5tJMoEJeYdsGEIWHGLWACEQbLwvMl7kBVveZO2b1VJ33/ePqm61Wl3dLVlSV0vn8zx+VH3r1r23ylJ9+5x77rlijEFRFEVRhhNPqgegKIqijD5UfBRFUZRhR8VHURRFGXZUfBRFUZRhR8VHURRFGXZUfBRFUZRhJynxEZGVIrJfRKpF5P4Y57NF5An7/CYRKY8494Bdvl9EbkjUpojMtNs4aLeZZZe/S0S2iYhfRG6L6v8uu/5BEbmr/49BURRFGU4Sio+IeIEfATcCC4A7RGRBVLW7gUZjzBxgLfCQfe0CYA1wIbAS+LGIeBO0+RCw1hhTATTabQMcA/4e+HXU+MYD3wQuA5YB3xSRomQfgKIoijL8JGP5LAOqjTGHjTFdwDpgdVSd1cBj9vFTwLUiInb5OmOMzxhzBKi224vZpn3NNXYb2G3eCmCMOWqM2QUEo/q+AXjRGNNgjGkEXsQSOkVRFMWlZCRRZypwPOJzDZaVEbOOMcYvIs1AsV2+MeraqfZxrDaLgSZjjD9G/f6Mr881InIPcA9Abm7u0vnz5ydoVlGUoDFUnWxh8rgcJuZnh8vPtvo43dLJhVMK8EgKBwi0dnZztL6DORPz6AoEOdbQQUVJPjmZOqU92GzdurXOGDNxMNpKRnxi/WpF5+RxquNUHuu3Il79eCR1jTHmEeARgMrKSrNly5YEzSqK0tLZzaJ/eoGvrbqAT149K1z+yGuH+JcN+/jrt24gNzuZ18jQ8eLeM3zq51tYd+9VHG/s4B9+tY2nvnA18yePS+m4RiIi8s5gtZXMV4MaYFrE5zLgpFMdEckACoCGONc6ldcBhXYbTn0NZHyKogwAYzu5PdL7O17oc9AFuSFD+SlFCFthLhiWkoBkxGczUGFHoWVhBRCsj6qzHghFmd0GvGys34j1wBo7Gm4mUAG85dSmfc0rdhvYbT6dYHzPA9eLSJEdaHC9XaYoynkSsN/i0a41CYvPcI+oL72H4B5RVOKTUHzs+Zd7sV7obwNPGmOqROTbInKLXe1RoFhEqoEvAffb11YBTwJ7geeAzxpjAk5t2m3dB3zJbqvYbhsRuVREaoDbgYdFpMruowH4ZyxB2wx82y5TFOU8Cdjq4vVEWz7WTzdkxQ8NQS2f9CIpZ60xZgOwIarsGxHHnViiEOvaB4EHk2nTLj+MFQ0XXb4Zy6UWq4+fAT+LexOKovSbkLh4+oiP9TngBtOHkHUmYYssFeLT3d1NTU0NnZ2dw9/5IJOTk0NZWRmZmZlD1kdqZwoVRXE1PW63KPHxuMftFoywfEKjNAnjlAafmpoa8vPzKS8vD4tgOmKMob6+npqaGmbOnDlk/WgsoqIojoRe7N4+AQfWT1e53RA8nt5lw0lnZyfFxcVpLTxgzecVFxcPuQWn4qMoiiPBYE8kWSQeVwUc9IxRUhxwkO7CE2I47kPFR1EURxIFHLghqqzH8ukRydSPSkmEio+iKI6ExCVafMRN63zsn9Ir4CD143Irr776KjfffDMAPp+P6667jsWLF/PEE08M6zg04EBRFEeC4QWcDotMozMtpoDIRabhgAPVnqTYvn073d3d7NixY9j7VstHURRHnAIOvJ7Q+dS/5SPdbiFRTP2oUsPRo0eZP38+d911F4sWLeK2226jo6OD5557jvnz53PVVVfxu9/9DoDa2lruvPNOduzYweLFizl06NCwjlUtH0VRHAnN+URnOHBVeh16rLOQRgZTHAnxrT9Wsfdky6C2uWDKOL75vgsT1tu/fz+PPvooV155JZ/4xCf4/ve/z8MPP8zLL7/MnDlz+PCHPwxASUkJP/3pT/nud7/LM888M6hjTQa1fBRFcSQsPo5zPsM+pD5owEFvpk2bxpVXXgnAnXfeyZYtW5g5cyYVFRWICHfeeWeKR2ihlo+iKI6YdFrn44JQ6xDJWChDRfT8XHNzsytDwNXyURTFkXCGg6g3RTi9jgvEJxiRhSH8jk39sFLGsWPHePPNNwF4/PHHue666zhy5Eh4Tufxxx9P5fDCqPgoiuJI0Cm9jpui3eyfVmLR0R1wAHDBBRfw2GOPsWjRIhoaGvjiF7/II488wqpVq7jqqquYMWNGqocIqNtNUZQ4BINO4mOfd4HlY2JYPm4YV6rweDz85Cc/6VW2cuVK9u3b16fuihUrWLFixTCNrDdq+SiK4ohzhoPUZY+OJhT04BHRLRXSCBUfRVEciXyxR+Jx0TqfYMQi09G+mVx5eTl79uxJ9TCSQsVHURRHeuZ8epe7Kr1OZLRbikOt3RD9NxgMx32o+CiK4ohTbjdXLTKNmPPxpFB9cnJyqK+vT3sBCu3nk5OTM6T9aMCBoiiOBMJbKkSl13HRItNI16CEy4Z/YGVlZdTU1HD27Nlh73uwCe1kOpSo+CiK4oiz5WOfd4H6hOd8SG0gRGZm5pDu/DnSULeboiiOhNbxRGc4cGN6HQ21Ti9UfBRFcSTQK5KsBzel1wlbPhFvs9SPSkmEio+iKI4YJ7ebxz3pdSItHzetP1Lio+KjKIojAdvt5phexwUv+cg5H3GRRabER8VHURRHAmHLp3e5q9Lr2D97WT6pG46SJCo+iqI4ErmGJpIe91bqX/ORGQ404CB9UPFRFMWRgGNiURdltY6xzke1x/2o+CiK4khoTic64MBNFkYwGGn5qNstXVDxURTFkcgXeySuSq9j/4xc5+MGd6ASHxUfRVEcCTiEWoc+uynazSOpzXCg9A8VH0VRHHHez8f66QbLJySAkuLcbkr/UPFRFMWRkPhkeHq/KtyVXseE3W09brfUjUdJDhUfRVEc8SewfNwwt2JMj7tN1/mkDyo+iqI4ErBjqTMc9vMJuMD0CRrTZ7M7dbu5HxUfRVEccbZ83ON2C5oeN6AnbJKlcEBKUqj4KIriSCAQmvOJTixq/XSDhWGMCQcaaMBB+pCU+IjIShHZLyLVInJ/jPPZIvKEfX6TiJRHnHvALt8vIjckalNEZtptHLTbzIrXh4hkishjIrJbRN4WkQcG+jAURelNIsvHDXM+ltvNGk8Kd9FW+klC8RERL/Aj4EZgAXCHiCyIqnY30GiMmQOsBR6yr10ArAEuBFYCPxYRb4I2HwLWGmMqgEa7bcc+gNuBbGPMQmAp8L8ixU9RlIETCFrzKdHbaLvJ7WYFHFjHus4nfUjG8lkGVBtjDhtjuoB1wOqoOquBx+zjp4BrxfptXQ2sM8b4jDFHgGq7vZht2tdcY7eB3eatCfowQK6IZABjgC6gJeknoCiKIwFj+oRZg/vW+YQtn3BZ6selxCcZ8ZkKHI/4XGOXxaxjjPEDzUBxnGudyouBJruN6L6c+ngKaAdOAceA7xpjGqJvQkTuEZEtIrLl7NmzSdy2oiiBoOnjcoOIdT4uMH2Cvdb5aKh1upCM+PT9zev7f+tUZ7DK4/WxDAgAU4CZwJdFZFafisY8YoypNMZUTpw4MUZTiqJE4w+YPsEG0DMH5IZQa2uRadScj1o+ricZ8akBpkV8LgNOOtWx3V8FQEOca53K64BCu43ovpz6+AjwnDGm2xhTC7wBVCZxX4qiJCAQDOL1xhCf0DofF7zjDT1uQN1SIX1IRnw2AxV2FFoWVgDB+qg664G77OPbgJeN9dVjPbDGjlSbCVQAbzm1aV/zit0GdptPJ+jjGHCNWOQCy4F9yT8CRVGc8AcdLB+vu9xufTIcqPq4noxEFYwxfhG5F3ge8AI/M8ZUici3gS3GmPXAo8AvRKQayxpZY19bJSJPAnsBP/BZY0wAIFabdpf3AetE5DvAdrttnPrAipr7v8AerC8+/9cYs2vAT0RRlDBOcz4hQfK7Qnzo43ZzwbCUBCQUHwBjzAZgQ1TZNyKOO7FCnmNd+yDwYDJt2uWHseZxostj9mGMaXPqW1GU88OyfGJFu4XmfFK/lanRgIO0RDMcKIriSCLLJ5B67em1zkcDDtIHFR9FURxxmvPxeNxj+fTKcGCXqfa4HxUfRVEcCQSDMS0fsKwft8z59N1SIfXjUuKj4qMoiiP+QGy3G1hrfQIuMDEisxlowEH6oOKjKIojTnM+YIuPSxb6hGIiNLdb+qDioyiKI1ZuN/dbPp6oxKea2839qPgoiuJIPMsnwyOuSK/TK7Fo7KEqLkTFR1EUR6zcbrFfE17XBBz0bCbncVHCUyU+Kj6KojiSaM7HDS95Q4/F46Z9hpT4qPgoiuKIPxgkI0ZiUYAMj8cVlo/pldvNKnPDXJQSHxUfRVEciWf5eDzu2FIhGIyc8xE8ohkO0gEVH0VRHHHKcACW5eMK8YnI7QaWELlhXEp8VHwURXEk4TofF7zkrTmfnjF6XBICrsRHxUdRFEecslqDtaGc3wW53UxEtBtY43JDIIQSHxUfRVEcCQRNOIloNJblM8wDikG0deb1iEa7pQEqPoqiOOIPBp3nfLziiqzWAUMvgRRxRyCEEh8VH0VRHAnESSzqEZcsMg0aIqPBLcsn9eNS4qPioyiKI91BQ6Y39msiwyUv+T5uN3HHuJT4qPgoiuJIdyBIlsMiU49H8Lsgq3UgKrGoiDvmopT4qPgoiuJItz8Y1/Jxw9xKsE/AgeZ2SwdUfBRFcaQ7YMhwEB+3bKkQMH3dbm4YlxIfFR9FUWJijKErjtvNLYtMg8HebjePS+ailPio+CiKEpOQsMRzu7llzifS8vHoItO0QMVHUZSYdNvCkpnh7HZzg4URiEgsCiF3YAoHpCSFio+iKDHpskPGnCwf12wmFzREDtEjGnCQDqj4KIoSk+6w+DjN+Xhc8ZLvE3DgEotMiY+Kj6IoMelOYPlkuMjy6RVwoFsqpAUqPoqixCQUTOAkPm55yccMOFDLx/Wo+CiKEpOuBG43tywy9QcMXtGs1umGio+iKDEJud2ynAIOvC5xu/WxfDSrdTqg4qMoSky6/dYL3DHDgbhkS4Wo9Dq6yDQ9UPFRFCUmidxurgm1Nr03vPO6ZC5KiY+Kj6IoMUnkdsv0uuMlHwj2nvNRyyc9UPFRFCUm/gQZDjK9nrBApZI+bjcBF3gDlQQkJT4islJE9otItYjcH+N8tog8YZ/fJCLlEecesMv3i8gNidoUkZl2GwftNrOS6GORiLwpIlUisltEcgbyMBRF6SHROh9LfAwmxVZG0MRKr6OWj9tJKD4i4gV+BNwILADuEJEFUdXuBhqNMXOAtcBD9rULgDXAhcBK4Mci4k3Q5kPAWmNMBdBotx2vjwzgl8CnjTEXAiuA7n4+B0VRokg055NlW0TdKU6kFuiTXscd7kAlPslYPsuAamPMYWNMF7AOWB1VZzXwmH38FHCtiIhdvs4Y4zPGHAGq7fZitmlfc43dBnabtybo43pglzFmJ4Axpt4YE0j+ESiKEotkMhxE1ksVgeiAA4+k3BpTEpOM+EwFjkd8rrHLYtYxxviBZqA4zrVO5cVAk91GdF9OfcwFjIg8LyLbROQrsW5CRO4RkS0isuXs2bNJ3LaijG6ScbtF1ksVweiAA91MLi1IRnxi2dzR/7NOdQarPF4fGcBVwEftn+8XkWv7VDTmEWNMpTGmcuLEiTGaUhQlkvCWCg5ut0y3uN1ipNdxQRyEkoBkxKcGmBbxuQw46VTHnoMpABriXOtUXgcU2m1E9xWvj78YY+qMMR3ABmBJEvelKEocElk+oR1OU2n5GGMwfQIOdEuFdCAZ8dkMVNhRaFlYAQTro+qsB+6yj28DXjaW03U9sMaOVJsJVABvObVpX/OK3QZ2m08n6ON5YJGIjLVF6d3A3uQfgeLzB6g62ZzqYSguw9dtiUp2nFBrSK34hAILNLFo+pFQfOz5lXuxXvJvA08aY6pE5Nsicotd7VGgWESqgS8B99vXVgFPYonBc8BnjTEBpzbttu4DvmS3VWy3Ha+PRuD7WIK2A9hmjPnTQB/IaORrv9/Dqh+8zpmWzlQPRXERnX4rbicn0xvzfIYbxMfEEB8NtU4LMhJXAWPMBix3VmTZNyKOO4HbHa59EHgwmTbt8sNY0XDR5fH6+CVWuLUyALYdawSgtbObSeN0iZRiEbJ8nDIchNxuXf7UvehDi0l7ud1E1O2WBmiGAyWMfllUIvH5g2RleHqFMUfiCrdb2PLpKdMtFdIDFR9FUWLS2R1wnO+BHvHxpzCXTSAQEp+ecYpuqZAWqPgoYST2F1xllOLzB8nOiD3fAz3ik0q3W3ewbxYGrwYcpAUqPoqixMTXHSAn0/kVkZWR+lDrUPLTjAjLx+uSHVaV+Kj4KIoSE8vycX5FhF74qRSfUN8Z3qjEoio+rkfFR1GUmPj8Accwa3BHwEFoM7tIt5tbtnpQ4qPioyhKTDq741s+IbdbVwrT6/hDlk+E2y3DJTusKvFR8VEUJSY+fyCpgAN/St1ufS0fr1fCc0GKe1HxURQlJp3dwbgBB+5wu/W1fDI9npSGfyvJoeKjhNHoVCWSZC2fVLrdQpZPZMBBhtdaZKpZDtyNio+iKDHx+RNZPnaotT+VodZ9M2+HLTK1flyNio8SRheZKpFYGQ7SI9otIyIFUOhY533cjYqPoigx8fmDZLt8zqcrZPlEROVlhAMhVHzcjIqPoigx6exOtM4nlNU69RkOMiMDDkLuQHW7uRoVH0VR+hAMGivaLc46HxEhO8NDpwvmfHoFHHjU8kkHVHwURenDuW5rI7nc7Phbfo3J8tJp100F3TEyHGS4YHtvJTEqPoqi9KG9yw/A2ATik5ORWvGJleEgJESa5cDdqPgoYXSdjxLiXJdt+WQ5z/kA5GR66OxO/ZxPbLebWj5uRsVHUZQ+tPss8RmblcDyyUy12y3WOp+Q202/TbkZFR9FUfrQYbvdcrPjWz7Zmd4UBxzEWueT+h1WlcSo+ChhdJGpEqK9K0nLJ8OTWssnHO0Wuc5HLZ90QMVHUZQ+dPiSs3xyMr34Uio+lsBkxUivo3M+7kbFR0G/HyrRtIcDDhLN+aQ24MDnt8YZue9QOL2ORru5GhUfRVH6EJrzGZsw2s0bXhOUCnz+IFleD57IOZ8Upv3pDgT50MNv8p9/PjjsfacbKj4KOtWjRJN0tFuK1/n4Yuy2Gl7nk4I5nwNnWnnrSAPfe/HAsPedbqj4KGF0nY8SoqPLjwhxt1SAkNstlZZPoE/y01RGu9W2+MLHje1dw95/OqHio+icj9KHjq4AuVkZSIIQyJwUh1r7/ME+2z6kcp1PbWtn+PhIffuw959OqPgoYVSElBCtnd0JI93AWufT5Q+mbNdQS3yiLB9v6iyfMxGWz5nmzjg1FRUfJTzno243JUTLOT8FYzIT1gu55Xwpsn46uwNkZzpYPn7rF/q5Pae599fbMMPwC17b2hleL3dSxScuKj5KGKO2j2LTfK47KfEZY7/4UzXvE8vyCe1B1GmHYX/197t5ZtcpGoZhDuZMi4+5JfnkZHo43XxuyPtLZ1R8FJUcpQ8tnd2My0ksPqFQ7I5UiU93oI/4hD777PVH9bbonG4ZekukttVHybhsSgvGqOWTABUfJYy63ZQQyVo+edlWnbZO/1APKSbWVt+93W4hyye0ADVEXdvQWz61LZ2U5OcwaVy2zvkkQMVH0TkfpQ8t57oZl4z45FjrgNp83UM9pJjEDDjwCB4hnHkhZJ3Vtfr6XD+YBIOGs60+Jo3LZvK4nGGxtNKZ+CvIlFGFzvkoYL1EW31+xuUkfj3k2RFxbb5Uzfn0dbtZ23t7w5ZPfk4GHV0B6tuHVnwaOrrwBw0l+dkEjKG2xYcxJmG4+mhFxUdRyVF60erzYwzJWT6pdrt1913nA1YUXigCL7TotH6I3W6hBaaTxuVggK5AkIb2Lorzsoe033QlKbebiKwUkf0iUi0i98c4ny0iT9jnN4lIecS5B+zy/SJyQ6I2RWSm3cZBu82sRH3Y56eLSJuI/GN/H4JioW43BSyXGyQpPil2u3V0+WPmn8uOSPsTsoCGes4ntMC0xHa7wfAEOaQrCcVHRLzAj4AbgQXAHSKyIKra3UCjMWYOsBZ4yL52AbAGuBBYCfxYRLwJ2nwIWGuMqQAa7bYd+4hgLfBssjeuKEpsmm3xSSrgICskPqlxu7V3BcjN7uvAyY6wfDrDUW9D63YLWT4l+TlMKrDE54yKjyPJWD7LgGpjzGFjTBewDlgdVWc18Jh9/BRwrViOztXAOmOMzxhzBKi224vZpn3NNXYb2G3emqAPRORW4DBQlfytKyE04ECJpKnDEp/CJMQnlAUhFW637kCQLn8wPO8USU5My2eIxce2fCbmR1g+zUPbZzqTjPhMBY5HfK6xy2LWMcb4gWagOM61TuXFQJPdRnRfMfsQkVzgPuBb8W5CRO4RkS0isuXs2bMJbnl0oZqjRBJ6SU/ITzxXkeH1MCbTS3vX8ItPuy+07YOz5RMImnCOt6Ge8znT4qNgTCY5mV4m5mcjom63eCQjPrFCNaLfV051Bqs8Xh/fwnLTtcU431PRmEeMMZXGmMqJEyfGqzpq0Wg3BSLEJze5ifLc7AxaU2D5hDa8y4vhdsvJ8OLrDoatHhFLfIYyxU5tayeTxlnPLNPrYUKervWJRzLRbjXAtIjPZcBJhzo1IpIBFAANCa6NVV4HFIpIhm3dRNZ36uMy4DYR+XegEAiKSKcx5odJ3JsSgbrdFLAm5jO9wrgxyQXD5udk0OZLoeUTw+2WnemhzecPZzmYUjCGE03n7BDyxO7EgXCmxUdJfk74s671iU8yls9moMKOQsvCCiBYH1VnPXCXfXwb8LKxvmKsB9bYkWozgQrgLac27WtesdvAbvPpeH0YY642xpQbY8qB/wD+RYVHUQZOfZuP4tzspNen5GVn0NY5/NFuIfGJGXCQ4cHXHQznd5taNAYYWtfbmZZOJo3rEZ9J43I04CAOCcXHtkDuBZ4H3gaeNMZUici3ReQWu9qjWPMv1cCXgPvta6uAJ4G9wHPAZ40xAac27bbuA75kt1Vst+3YhzJ4qOGjgJULbUJ+VtL1C8dm0nQuFeJjCUtuzDkfL53+QNjymVpoic9QBR0EgobaVh+TC3pclZMLstXyiUNSdrUxZgOwIarsGxHHncDtDtc+CDyYTJt2+WGsaLjocsc+Iur8U7zzSnyGI+W84n7qbMsnWYrGZnG8oWMIRxSbtrDlE2udj8ee8+ktPvVDJD71bT4CQROOcgPL7dbU0U1ndyCcb07pQXO7KWFUehSwXFPFeclbPuNzs4Zlu4JoQuITK+AgLzuD9i5/OOBgStjyGZpxnrIDCyYXjAmXhVxw6nqLjYqPoihhLPdR77mLRBSNzaKl0093YHg3lGvqsISkcGxfoczNzqDd5w9bPqWF1v0M1ZxPyL3Wy/IpCK31UfGJhYqPEka9bkpdm4/ugAlbCskwPteKHmvsGF7rp7GjC69HYiZAzcvOoDtgaLUDIcZmeikcmzlkWQ5ONVkbx4UEB6DUtoJONOmmcrFQ8VEiUPUZ7YRelFML+2H55FqWR2P78AYdNHZ0UzgmM2ZUXsgVF7J0sjO9FOdmDVnAwZG6dvKyM5gQ4a6cPn4sXo9w+Gz7kPSZ7qj4KGHU8lFO2uLTL8vHdnsN97xPU0cXhWNjr9kJhV+HxpTl9VCclz1kcz6H69qZOSG3lxBmZXiYPn4sh+virn8ftaj4KIoSZiDiE7Z8htvt1t5NUYz5HuixfELik53pYUJe1pBFux2pa2fWxNw+5bMm5Krl44CKjxL2tqnho5xs6iQ/O6NfWQBCkXFnh3in0GgaO7piBhtAhNstwvIpyc/hdHPnoC8paO7opqbxHHMn5fc5N7skj8N17QSC+tcVjYqPEkbdbsqxho5wNoBkmZCbTaZXONk8vBPrdW1dveZYIgntMxRp+cwpyaO9KzDoAQA7apoAuGRaYZ9zsyfm0uUPUtM4/Oug3I6KjxJGF5kqh8+2MXtiXr+u8XiEyQU5nGoavpDiLn+QujZfOKIsmtA2C2Hx8XqZP9myTPafbh3UsWw6XI/XIywsK+hzbkGpVbb7RPOg9jkSUPFR1N2mANa+N8caOmLOXSSitGAMp4bR8gkt3CwtiB2V1yfgIMPDXFt89p2H+FTXtvHwXw7xTr01jxMMGjbsPsXls4rJj+GqnDc5nyyvh101Kj7RJJe2VhnRhCweFaHRzbH6DoKGfls+AFMKcth8tHEIRhWb8KJOB/EJzVmFNnjLyvAwxuNlRvFYdhxvGlCfbT4/ax7ZSF2bj++9cIBPr5jNuJwMjtZ38MX3zo15TVaGhwtK89k5wD5HMio+Shj1uo1uDtlRWQOyfArHcKblFMGgweNJLhv2+RBKZ+Nk+YzN8pKT6aGzO0iGR/DaY1o+s5hn95wiEDThsmR5oeo0dW0+/vOOS3hx7xl+8OeDVpuzxnPzoimO1y0qK+R322rwB4JkeNXZFEKfhBJGN5Mb3Rw404oIzBqA5VNWNAZ/0HBqmPKYHbPdXk7BESISTo46JiKp5+Wzi2np9LP3ZEu/+9xxvIncLC83LSzlB3dcwlOfvpwf3HEJ//PxZXGFbNnM8bR3Bdil8z69UPFRVHIUAHbVNDNzQm7MRJ2JCLnqDtUOz4LK6to2phaOibmFdohQJFxeRPqdqyom4PUIz1Wd6nefO483sbCsICw0leXjueXiKQkzVl85ZwIi8PrBun73OZJR8VF63G2qQqOaPSeaWTi1b8RWMswpscSnepjE59DZ2Is6IynOsyyfyM3mJuRlc+WcCTy94yTBfqy96ewOsPdUC4unFfV7rONzs7hoSgF/PXi239eOZFR8lDCqPaOX2tZOTrd0Dlh8inOzKBybSfXZoRefQNBwKImQ8PF25oVoS+6DS6ZS03iOP++rTbrPt0+10B0wLJ42sOezYsJa3sgAABu+SURBVN5Etr7TOGS55dIRFR9F53oUdhyzorEWlfVdKJkMIsKciXkcPDO4a2hicbC2lY6uAItirKuJpCS/75wPwKqFpZQVjeGHLx9M2voJRasNxPIBuHnRFIIGNuzuv7tvpKLio4TdbhrtNnp583A92RkeLh7gN3uAi6YWsOdEC/4h3tdn2zuWECyZHl8IKiZZllFoQ7kQGV4Pn7+2gp01zTyx5XhSfe443sSkcdmOod2JmDc5n7mT8vjjzpMJ6xpjRsWCbxUfJYxaQKOXNw/Vc2n5eLIzBr7d85IZRZzrDpzXIs5k2HSkngl5WcwoHhu33iXTivAI3HJx3zDo25aWsXzWeB7809tU1yYe786aZi4eoFUYYvXiqWw+2sihOK7J7kCQm//zdW754Rt9RHOkoeKjhBkFX7aUGNS1+dh3upXLZxefVztLplsv563vDN1i0+5AkFf21bJiXknMfXwiKZ+Qy9avvZe7rijvc05E+N6HFpOT6eGTj20JL0aNRVNHF0fq2rk4Ru62/nB7ZRmZXuEXb77jWOe1A2epOtnC7hPN/HHnyHbRqfgoKjqjnJf2ngGsSfHzYWrhGKYWjuH16qELKX69uo6WTj/vXTApqfpFuVmOIjW1cAwP/91Salt9rHlko+N219uOWWK6dMbA5ntClOTncNPCUn67tYY2nz9mnZ3HmxCxtuP+3baa8+rP7aj4KGFUg0YnG/acZkbxWBaUjjuvdkSE6y4o4a8Hz3Kua2hcRj//21Em5GXznnklg9Le0hnjeewTyzjT3MmtP3qDXTV90+BsOdpIhkfO2+0GcNcV5bT6/DyxOfZc077TrcyckMuHKsvYeLh+2LepGE5UfJQwo2GSU+lNY3sXf6uuY+VFkxO6sZLh+gsn09kd5OV+hDEny/Zjjbyy/yx3Lp9OVsbgvbouLR/Pbz59BV6PcPtP3uS3W2t6/S28vK+WS6YXMiZr4PNhIZZML+KymeN55LVDdHb3Fej9Z1q5YPI4VtnRcc9VnT7vPt2Kio+iiUVHMb/Zehx/0PD+S6YOSnvLZxUztXAMv9h4dFDaC9HZHeCrv99DSX42n7x61qC2DbBgyjjW33sli6cV8uXf7OSzv95GXZuPLUcb2He6lVULSwetr89dW8GZFh+/iYq0a/f5eae+IxwZN6ckjz/tShwdl66o+CgqOqOUYNDwq03HuLS8iPmTz8/lFsLrET52+Qw2Hm5g6zsNg9KmPxDkK0/tYu+pFv7l/QsHlP4nGYrzsvn1p5Zz38r5vLj3DFf868vc+egmSgty+NCl0watnytmF7N0RhH/9eohuvw9YekH7DVS8ybnIyKsWljKpiMNcYMh0hkVH6UHVaFRxQt7z/BOfQd3Lp8xqO3euXwGk8Zl842nq3q9XAdCXZuPT/18C+t3nuS+lfO5LslAg4Hi9QifWTGb57/wLj52+Qzef0kZv/7U8rg55PqLiPC5ays42dzJE5uPhctDm9yFNr1btagUY+C5PSPT9abio/QsMlX1GTUEg4a1Lx5g1sTcQXUpgZVL7Vu3XEjVyRa+9ofd/cqhFuJcV4D/fu0w1699jTcO1fOdWy/iMytmD+o44zFrYh5fu3kB//qBhcyc0P8tJhLxrooJXDZzPP/nzwfDkW/7z7QyJtPLtCJr/dLcSflUlOTxp10jM+RaxUcJo/EGo4enttWw/0wrn7+2Ykj2mFl5USmfu2YOT26p4ZM/30JNY0fCa9p9fl7dX8tXntrJsgdf4sENb7OgdBx/vPeqQbfOUo2I8MBNF1DX1sUjrx0GrPxxcyfl9doPadWiUt462kDtMG1VMZzoZnKKWjyjjLo2H/+y4W0qZxTxvjiboJ0vX3zvXMbnZvFvz+3j3f/7Va6umEDljCImF4xhTKaX9i4/je1dHDrbxv4zbVSdaMYfNORmeVl5USl3LJtGZfn4IRtfqlk8rZBVi0r579cO88ElU9l+rImPXDa9V51VC0v5j5cO8uye0zEXy6YzKj5KGLV8Rj7BoOH+3+6m3efnXz+wcEh3HRUR/v7KmVy3YBK/3HiMF6pO8+r+vtsKTMjLZk5JLp961yyumF1M5YzxgxLWnA7cd8N8XtlXy3u//xpdgSArotYvVUzKZ96kfP6065SKjzLy6JnzUUY6D792mJfePsM3bl5AxaT8YemzrGgs9984n/tvnE+bz099m49z3QFyszIoHJtJfk7msIzDjUwvHsv3P7SYLz+5g+WzJ3L1nAl96qxaVMralw5wurlzwIlN3YiKjxJGF5mObJ7ecYJ/f34fqxaV8vEry1MyhrzsjCELlU5XVl40mfcuuAGPEHOh700LS/n+iwd4ds8pPn7lzBSMcGjQgANFLZ5RwLO7T/GlJ3dy2czxfO/2iwclm4EyeHg94vh/Mqckj/mT81mfxHYM6YSKj6JutxHO/7xxhH/49TYWlRXw07suJSdzdMynjCQ+uKSM7cea2He6JdVDGTSSEh8RWSki+0WkWkTuj3E+W0SesM9vEpHyiHMP2OX7ReSGRG2KyEy7jYN2m1nx+hCR94rIVhHZbf+8ZqAPY7SjXreRxbmuAA/8bhf/9Me9XHfBJH79yeXq8kpTPri0jKwMD7/edCxx5TQhofiIiBf4EXAjsAC4Q0QWRFW7G2g0xswB1gIP2dcuANYAFwIrgR+LiDdBmw8Ba40xFUCj3bZjH0Ad8D5jzELgLuAX/XsEito8I4+qk83c8sPXWbf5OJ9ZMZuf3Ll01ESQjUTG52axamEpv9t2gnaH7RjSjWQsn2VAtTHmsDGmC1gHrI6qsxp4zD5+CrhWLAfmamCdMcZnjDkCVNvtxWzTvuYauw3sNm+N14cxZrsxJuQMrQJyRCQ72QegRFo8KkLpTpvPzz8/s5dbfvgGjR3d/PwTy7hv5Xy8QxhSrQwPH71sOm0+f1JbcacDydjgU4HI9Ks1wGVOdYwxfhFpBort8o1R14bS58ZqsxhoMsb4Y9R36iNy56oPAtuNMX02wRCRe4B7AKZPnx59elQTkhx1u6Uv3YEgv91aw9qXDlDb6uOOZdP5yg3zKBybleqhKYPE0hlFzJuUz682HWPNsvR/hyUjPrG+MkW/ppzqOJXHsrji1U84DhG5EMsVd32MehhjHgEeAaisrNTXbASBoG6pkK4Eg4Y/7jrJ2hcPcLS+g4unFfJfdy5lyfTz23VTcR8iwkeXT+cbT1exq6aJRYOwuV0qSUZ8aoDIfOJlQLTdF6pTIyIZQAHQkODaWOV1QKGIZNjWT2R9pz4QkTLg98DHjDGHkrgnJYKBJH5UUsu5rgBPbT3Oz944ypG6duZNyue/P1bJdReUaBj1CObWS6byrxv28auNx1h0W3qLTzJzPpuBCjsKLQsrgGB9VJ31WJP9ALcBLxtrxeJ6YI0dqTYTqADecmrTvuYVuw3sNp+O14eIFAJ/Ah4wxrzRn5tXLAKhzeRUg1zP6eZOvvv8fi7/tz/z9aerGJeTwQ8/cgnPfv5q3rtgkgrPCGdcTiarF09h/c6TtHR2p3o450VCy8eeX7kXeB7wAj8zxlSJyLeBLcaY9cCjwC9EpBrLGlljX1slIk8CewE/8FljTAAgVpt2l/cB60TkO8B2u22c+gDuBeYAXxeRr9tl1xtjBn8f3xFKMLyTqaqPG+kOBHllXy1PbD7OK/trMcB1F0ziU1fP4tLyIhWcUcYdy6azbvNxnt5xkr9L42zfMhpTqlRWVpotW7akehiuYe5Xn6UrEOQ/77iE9108dFmOlf5x+GwbT22t4Tdbazjb6mNifja3Ly3jw5dOY0bx4O8xo6QHxhhu+sHrCPCnz101rF8+RGSrMaZyMNrSFWdK2O2mpJ5Tzed4Zucp1u88ye4TzXgErplfwocvnc575k0ckr13lPRCRPjIsml8/ekqdp9oTtvAAxUfRaPdUkxDexcbdluCs/loA8bAorICvrbqAt538RQmjRs5mYyVwWH1JVN5cMPbPP7WcRUfJf0ZjS7YVHG6uZMX957m+aozbDxcjz9omD0xly9cO5dbFk8Zkq2blZHDuJxMbl40hfU7TvC1VReQm4Zpk9JvxIqSplTXtvF81Wle2HuGncebAJg1IZdPXj2L911cyoLScRo8oCTNHcum8dTWGv6482RaLjpV8VGUISIQNOw43sRLb5/h+arTHD7bDsDFZQX8/zfM44YLJzF7Yp4KjjIglkwvYu6kPB7ffFzFR0lv1Ot2/tS3+Xjt4Fle3X+Wvxw4S1NHNxkeYfmsYv7+inLeu2ASpQVjUj1MZQQgIqy5dDrffmYve0+2sGDKuFQPqV+o+ChhdJ1P/wkGDbtPNPPq/rO8sr+WnTVNGAPFuVlcM7+EFfNKeHfFRArGjt6topWh4wNLpvJvz+7jN1uP880pF6Z6OP1CxUcJo5ZPctS3+XjjUD2v7q/lL/vPUt/ehQhcXFbIF66dy4p5E1k4tQCPZpJWhpjCsdaXnD/uPMlXb7ogrULxVXwUJQHnugJsPtrAG9V1/PVgHXtPWbtJFo7N5N1zJ/KeeSVcXTGB4jzdyUMZft6/ZCrPVZ3mr9V1vGdeSaqHkzQqPkoYtXwsAkHDnhPNvF5dx+sH69j6TiNdgSCZXmHpjCL+8fq5XDlnAovKCnWfHCXlvGdeCYVjM/n9thMqPsrwUtvSyUPP7ec7t17U790quwPB8PFo1R5jDO/Ud4TF5m+H6mjptLaUuqB0HHddMYMr50xg2czxjM3SPxnFXWRleLh5USlPba2hzedPm63S02OUSlzWvnSQ326rYfms8dxeOS3xBRH4/MHElUYYxhiO1new8XA9mw7Xs/FwA6dbOgGYUpDDyosmc1XFRK6YXcwEdaUpacD7L5nKLzce49ndp/r9DkgVKj4jgDZ7T/esjP5PNvq6A+HjkZrhwBjD4bp2W2wa2Hi4ntpWa7PbCXnZLJ81nstmFXPVnAmUF4/VdTdK2rFkehEzisfy++0nVHyU4eNclyUgOZn9c7kBdPpHntvNGMOhs228ebghbNnUtVliU5KfzfJZxVw2azzLZxUza0Kuio2S9ogIty6eyg9ePsjp5k4mF7g/H6CKzwigM8J66S8dttUEpK36+ANB9p1uZfPRBrYcbWTTkXrq2roAmDwuh6vmFHPZrGKWzypWy0YZsbzv4in8nz8f5Nk9p/j4lTNTPZyEqPiMANq7LAHpGsD8TfO59NsNsc3nZ/uxRrYcbWTrO41sP9ZIu239TSnI4eqKiSy3LZvp41VslNHBnJI85k/OZ8NuFR9lmOjwWS/eyMi1ZGnq6BEft2Y4ONV8ji1HG9lytIEt7zTy9qkWggZEYP7kcXxwaRlLZxRRWT6eqYWaukYZvdy0sJTvv3ggLVxvKj4jgJDlMyDxibB83BBvEAgaDpxpDQvNlqONnGg6B8CYTC+XTC/k3vfMobJ8PJdMLyQ/R9PWKEqIkPikg+tNxWcE0G7P23QF+q8eZ+wQY0jNlE9Hl58dx5vYerSRze80sv2dRlrt+ynJz+bS8vHcfdVMKsuLuKB0HJlplD5EUYabkOvtT7tUfJQhxhhDo+066x7AnM+x+o7BHlJcTjadY+s71lzNtmONVJ1sIRA0iMDcknzet3gKl5YXUTljPGVFY3S+RlH6yaqFpXwvDVxvKj5pTiiqC6BrAG63g7WtTBs/huMN5wbd7dYdCPL2qZaw2Gx9p5FTzZalNSbTy+JphXzm3bNZOqOIJdOLNPOzogwCNy2yxMftrjcVnzRnh70jJvTf8jnXFWDPiRZuWjjZEp/zdLw1dXSx7ZglMluONrKzponObmtMUwpyqCwfz9LphSydMZ75pfnqQlOUIWD2xPRwvan4pDkbdp8iN8tLe1eg3wEHz+45RVcgyHvml/CHHSf7bfmcaelko72I860j9Ryyd+r0eoQLp4zjjmXTw1bNFI1CU5RhIx1cbyo+acyLe8/whx0nuPvKmfx84zv9Cjho6ujiey8cYO6kPK6cMwGwIs3i0dDexV8PnuXNQ/VsOtLAkTpLbPKzM6gsL+IDS6yQ50VlBZqAU1FSSDq43vQNkYa0+fz88OVqHnntEBdNKeBL189l3ebjSVs+B8+08v89vp2zbT5+cMdyxthpeXz+vpkSjta188yuk7y8r5btx61dOsflZLBs5ng+smw6y2cVs2DKON1aQFFcRMj15uYFpyo+aURLZzePbzrGT18/wtlWH7ctLeNbt1zI2KwMMr2SMMPBwTOtPPzaYf6w/QTjxmTy049VsnRGEX5btHz2/EwwaHhh7xkeff0wm482ArCorIDPXVPBe+aXsHBqgYqNoricmxaWsvalA5xp6WTSOPe53lR8XE4waNh6rJHfbz/B+h0nafP5uXxWMQ//3VKWTC8K18v0emJaPrUtnazfeZI/7DjBnhMt5GR6uHP5DD77njlMzLe2C8jwevB6BJ8/yPGGDr7wxA62vtNIefFYvrJyHu+/ZCqlBTpnoyjpRHjB6e5T/L0LrR8VHxdyrivAxsP1/OXAWV7ce4YTTefIyfSw8sLJ3H3VLBaWFfS5JivDEw61Pt7QwYbdp3i+6nTYVbZwagFfv3kBqxdPiblHTXaGh4O1rXzwv/7Gue4A//7BRXxgydS02hNeUZQe5pTkMW9SPht2n1bxUWLT1NHFzppmdh5vYvPRBjYdaaDLHyQn08Pls4r5xxvmcv2CyeTG2aEwy+uhsb2Lz6/bzvqdVuTaRVPH8cXr5nLTwlLmlOTFHUNWhofnq84wLieD337mCuZOyh/s21QUZZi5aWEp//HnA9S2dFLiMtebis8w4/MH2HuyhZ3Hm9hxvImdNc3hqLHQKv+PLZ/Bu+dN5NLy8Unv0ZPp9fDK/rNkeIRPv3s2H1k2nWnjxyY9rhY7x9v9N16gwqMoI4SbFk5m7UsHeHbPae66ojzVw+mFis8QEgwajtS39wjN8Sb2nmqh2w6JLsnPZvG0Qm6vLGNxWSELywoGnCiz6ZyV6eAzK2bz5evn9fv6L1w3l6qTzay5ND12QVQUJTEVk/KpKMljw+5TKj4jmeZz3ew43sQ2O2/ZzuNNtHRaSTLHZnlZVFbAJ66aySXTCrl4WuGgTuLfcvEU/vx2Lf+wYs6Arv/ctRWDNhZFUdzDTQtL+cHLB6lt7aQk3z2uNzFuyKM/zFRWVpotW7YMSltVJ5t5fs9pXtl/lt0nmgHwCMydlM8l0wtZPK2QxdOKmFOSN+ThycGgwaMh0IqiRHDgTCvXr32Nf159IX93efl5tSUiW40xlYMxLrV8Bsgr+2tZ++IBdtU04xFYMr2IL143l8ryIi6eVkhenOCAoUKFR1GUaOZOymdOSR5/2n3qvMVnMFHx6SeBoOGb6/fwy43HKC8ey7duuZBbLp5CUW5WqoemKIoSk5sWlvJDl7neklrEISIrRWS/iFSLyP0xzmeLyBP2+U0iUh5x7gG7fL+I3JCoTRGZabdx0G4za6B9DAU/+cshfrnxGPe8axYvfPHd3HVFuQqPoiiu5paLSwka+M2WmlQPJUxC8RERL/Aj4EZgAXCHiCyIqnY30GiMmQOsBR6yr10ArAEuBFYCPxYRb4I2HwLWGmMqgEa77X730d8HkQyd3QF+8pdDXL9gEg/cOJ+sDF2AqSiK+5lTks+75k7kf/52lM7uvjkcU0Eyb89lQLUx5rAxpgtYB6yOqrMaeMw+fgq4VqwtKFcD64wxPmPMEaDabi9mm/Y119htYLd56wD7GHSO1rfT2unn5oun6A6biqKkFZ9dMZuzrT6+/JudqR4KkNycz1TgeMTnGuAypzrGGL+INAPFdvnGqGun2sex2iwGmowx/hj1B9JHGBG5B7jH/tgmIvVAneNdx2H1QwO5ytVMYIDPYgSiz8JCn0MPI+pZ/Bj48UcHdOkEYMZgjSMZ8Yn1FT86PtupjlN5LIsrXv2B9NG7wJhHgEdCn0Vky2CFDKY7+ix60Gdhoc+hB30WFvZzKB+s9pJxu9UAkcvey4CTTnVEJAMoABriXOtUXgcU2m1E99XfPhRFURSXkoz4bAYq7Ci0LKzJ/fVRddYDd9nHtwEvG2v16npgjR2pNhOoAN5yatO+5hW7Dew2nx5gH4qiKIpLSeh2s+dX7gWeB7zAz4wxVSLybWCLMWY98CjwCxGpxrJG1tjXVonIk8BewA981hgTAIjVpt3lfcA6EfkOsN1um4H0kYBHElcZNeiz6EGfhYU+hx70WVgM6nMYlel1FEVRlNSiC1UURVGUYUfFR1EURRl2RqX4JEoXNBIQkZ+JSK2I7IkoGy8iL9qpi14UkSK7XETkB/bz2CUiSyKuucuuf1BE7orVl5sRkWki8oqIvC0iVSLyebt8VD0LEckRkbdEZKf9HL5ll7s6ndVQYmdb2S4iz9ifR+WzEJGjIrJbRHaIyBa7bOj/Powxo+ofVoDDIWAWkAXsBBakelxDcJ/vApYAeyLK/h243z6+H3jIPr4JeBZrzdRyYJNdPh44bP8sso+LUn1v/XwOpcAS+zgfOICV0mlUPQv7fvLs40xgk31/TwJr7PKfAJ+xj/8B+Il9vAZ4wj5eYP/NZAMz7b8lb6rvb4DP5EvAr4Fn7M+j8lkAR4EJUWVD/vcxGi2fZNIFpT3GmNewogIjiUxRFJ266OfGYiPWWqtS4AbgRWNMgzGmEXgRK39e2mCMOWWM2WYftwJvY2XAGFXPwr6fNvtjpv3P4OJ0VkOJiJQBq4Cf2p9dndorBQz538doFJ9Y6YL6pOMZoUwyxpwC66UMlNjlTs9kRD0r211yCda3/lH3LGw30w6gFuvlcIgk01kBkems0vo52PwH8BUgaH9OOrUXI+9ZGOAFEdkqVhoyGIa/j9G4n09S6XhGGeeVuigdEJE84LfAF4wxLeKcGHbEPgtjrX9bLCKFwO+BC2JVs3+O2OcgIjcDtcaYrSKyIlQco+qIfxY2VxpjTopICfCiiOyLU3fQnsVotHxGczqeM7aJjP2z1i7vbxqktEJEMrGE51fGmN/ZxaPyWQAYY5qAV7F89qMxndWVwC0ichTL7X4NliU0Gp8FxpiT9s9arC8lyxiGv4/RKD7JpAsaqUSmKIpOXfQxO5JlOdBsm9rPA9eLSJEd7XK9XZY22L75R4G3jTHfjzg1qp6FiEy0LR5EZAxwHdb816hLZ2WMecAYU2asJJlrsO7to4zCZyEiuSKSHzrG+r3ew3D8faQ60iIV/7AiNg5g+by/murxDNE9Pg6cArqxvpXcjeWn/jNw0P453q4rWJv7HQJ2A5UR7XwCayK1Gvh4qu9rAM/hKizzfxeww/5302h7FsAirHRVu+yXyzfs8llYL8xq4DdAtl2eY3+uts/Pimjrq/bz2Q/cmOp7O8/nsoKeaLdR9yzse95p/6sKvQ+H4+9D0+soiqIow85odLspiqIoKUbFR1EURRl2VHwURVGUYUfFR1EURRl2VHwURVGUYUfFR1EURRl2VHwURVGUYef/AXsV0R4QSSFAAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], "source": [ "plt.clf()\n", "# plt.plot(x_part, calcs, '.')\n", "plt.plot(test_q, calcs_test, label = 'pdf')\n", "# plt.plot(test_q, res_y, label = 'res')\n", "plt.legend()\n", - "plt.ylim(0.0, 5e-6)\n", + "plt.ylim(0.0, 10e-6)\n", "# plt.yscale('log')\n", "# plt.xlim(3080, 3110)\n", "plt.savefig('test3.png')\n", @@ -1327,7 +1240,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 38, "metadata": {}, "outputs": [], "source": [ @@ -1341,7 +1254,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 39, "metadata": {}, "outputs": [], "source": [ @@ -1362,7 +1275,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 40, "metadata": {}, "outputs": [], "source": [ @@ -1372,10 +1285,60 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 41, "metadata": {}, "outputs": [], - "source": [] + "source": [ + " def jpsi_res(q):\n", + " return resonance(q, _mass = jpsi_mass, scale = jpsi_scale,\n", + " phase = jpsi_phase, width = jpsi_width)\n", + "\n", + " def psi2s_res(q):\n", + " return resonance(q, _mass = psi2s_mass, scale = psi2s_scale,\n", + " phase = psi2s_phase, width = psi2s_width)\n", + " \n", + " def p3770_res(q):\n", + " return resonance(q, _mass = p3770_mass, scale = p3770_scale,\n", + " phase = p3770_phase, width = p3770_width)\n", + " \n", + " def p4040_res(q):\n", + " return resonance(q, _mass = p4040_mass, scale = p4040_scale,\n", + " phase = p4040_phase, width = p4040_width)\n", + " \n", + " def p4160_res(q):\n", + " return resonance(q, _mass = p4160_mass, scale = p4160_scale,\n", + " phase = p4160_phase, width = p4160_width)\n", + " \n", + " def p4415_res(q):\n", + " return resonance(q, _mass = p4415_mass, scale = p4415_scale,\n", + " phase = p4415_phase, width = p4415_width)" + ] + }, + { + "cell_type": "code", + "execution_count": 59, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "2.404870995999285" + ] + }, + "execution_count": 59, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "_0 = jpsi_scale*np.cos(jpsi_phase)*jpsi_width/jpsi_mass**3 + psi2s_scale*np.cos(psi2s_phase)*psi2s_width/psi2s_mass**3\n", + "_1 = p3770_scale*np.cos(p3770_phase)*p3770_width/p3770_mass**3 + p4040_scale*np.cos(p4040_phase)*p4040_width/p4040_mass**3\n", + "_2 = p4160_scale*np.cos(p4160_phase)*p4160_width/p4160_mass**3 + p4415_scale*np.cos(p4415_phase)*p4415_width/p4415_mass**3\n", + "\n", + "R = (0.1/(1300**2) - ((_0 + _1 + _2)))\n", + "\n", + "R*10*2010**2" + ] }, { "cell_type": "code", diff --git a/__pycache__/pdg_const.cpython-37.pyc b/__pycache__/pdg_const.cpython-37.pyc index 6b27d3c..e91f5c6 100644 --- a/__pycache__/pdg_const.cpython-37.pyc +++ b/__pycache__/pdg_const.cpython-37.pyc Binary files differ diff --git a/data/zfit_toys/toy_0/0.pkl b/data/zfit_toys/toy_0/0.pkl index 7b1e788..b6df2c8 100644 --- a/data/zfit_toys/toy_0/0.pkl +++ b/data/zfit_toys/toy_0/0.pkl Binary files differ diff --git a/data/zfit_toys/toy_0/1.pkl b/data/zfit_toys/toy_0/1.pkl index 1376396..9968cdd 100644 --- a/data/zfit_toys/toy_0/1.pkl +++ b/data/zfit_toys/toy_0/1.pkl Binary files differ diff --git a/data/zfit_toys/toy_0/2.pkl b/data/zfit_toys/toy_0/2.pkl index c787017..c020748 100644 --- a/data/zfit_toys/toy_0/2.pkl +++ b/data/zfit_toys/toy_0/2.pkl Binary files differ diff --git a/data/zfit_toys/toy_0/3.pkl b/data/zfit_toys/toy_0/3.pkl index 15f1a4d..4b4828d 100644 --- a/data/zfit_toys/toy_0/3.pkl +++ b/data/zfit_toys/toy_0/3.pkl Binary files differ diff --git a/data/zfit_toys/toy_0/4.pkl b/data/zfit_toys/toy_0/4.pkl index f5ef0a4..967cb5a 100644 --- a/data/zfit_toys/toy_0/4.pkl +++ b/data/zfit_toys/toy_0/4.pkl Binary files differ diff --git a/data/zfit_toys/toy_0/5.pkl b/data/zfit_toys/toy_0/5.pkl index f027626..8c05273 100644 --- a/data/zfit_toys/toy_0/5.pkl +++ b/data/zfit_toys/toy_0/5.pkl Binary files differ diff --git a/pdg_const.py b/pdg_const.py index c14d561..52f9764 100644 --- a/pdg_const.py +++ b/pdg_const.py @@ -142,5 +142,10 @@ "p4040_BR": 4.2e-10, "p4160_BR": 2.6e-9, "p4415_BR": 6.1e-10, + +# Estimates + "Dbar_scale": 1.25, #with phase = pi + + "DDstar_scale": 2.41, #with phase = pi } diff --git a/raremodel-nb.ipynb b/raremodel-nb.ipynb index 7cf7c7e..4688758 100644 --- a/raremodel-nb.ipynb +++ b/raremodel-nb.ipynb @@ -9,9 +9,31 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "metadata": {}, - "outputs": [], + "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", @@ -42,7 +64,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ @@ -61,7 +83,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ @@ -262,7 +284,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ @@ -317,7 +339,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 5, "metadata": {}, "outputs": [], "source": [ @@ -345,7 +367,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 6, "metadata": {}, "outputs": [], "source": [ @@ -422,7 +444,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 7, "metadata": {}, "outputs": [], "source": [ @@ -450,9 +472,19 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 8, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "WARNING:tensorflow:From C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\ops\\resource_variable_ops.py:435: colocate_with (from tensorflow.python.framework.ops) is deprecated and will be removed in a future version.\n", + "Instructions for updating:\n", + "Colocations handled automatically by placer.\n" + ] + } + ], "source": [ "#rho\n", "\n", @@ -548,7 +580,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 9, "metadata": {}, "outputs": [], "source": [ @@ -579,7 +611,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 10, "metadata": {}, "outputs": [], "source": [ @@ -622,9 +654,22 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 11, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ8AAAD8CAYAAACo9anUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO2deXyc1XX3v2dmtEuWLFnyIsmWjGW82xhhCBDCjllNAiROQkqWt7Qp9H0b2rdAszSh4W1J05C0JSEkJKEkwRAgwQGHfQkJixdsvIGxLG+yjSVr32e77x/zzGg0mk2yNJvO9/PRxzN37nPuncfS/Oace+65YoxBURRFURKJLdkTUBRFUSYfKj6KoihKwlHxURRFURKOio+iKIqScFR8FEVRlISj4qMoiqIknLjER0RWi8geEWkQkTvCvJ4jIo9ar78tIjVBr91pte8Rkcti2RSRWsvGXstmdhxjLBORN0Vkl4jsEJHcsdwMRVEUJTHEFB8RsQP3AZcDi4BPi8iikG5fAtqNMfOAe4F7rGsXAWuBxcBq4IciYo9h8x7gXmNMHdBu2Y42hgP4JfDXxpjFwPmAa5T3QVEURUkg8Xg+q4AGY0yjMcYJrAPWhPRZAzxkPX4cuEhExGpfZ4wZNMbsBxose2FtWtdcaNnAsnltjDEuBbYbY94FMMa0GmM88d8CRVEUJdE44uhTCRwOet4EnBmpjzHGLSKdQJnV/lbItZXW43A2y4AOY4w7TP9IY8wHjIg8B5TjE7vvhL4JEbkZuBmgoKDg9AULFsTx1hVl8nKiZ5BjnQMsnjUFm8iI1z/sGqCle5CllcVJmN0QO450UlGUw/Qpvmh7Y0svAHPLC5I5rYxky5YtJ4wx5eNhKx7xGflbB6E1eSL1idQezuOK1j/aGA7gXOAMoA94SUS2GGNeGtbRmAeABwDq6+vN5s2bw5hTFMXPj17dxz3Pvs8b/7Ka3Cz7iNe/9/we/vPlBjb96xVIGHFKBB6v4ZR/2sBtl8znf19UB8DaB97Ea+Cxv/pIUuaUyYjIwfGyFU/YrQmoDnpeBRyN1MdagykG2qJcG6n9BFBi2QgdK9oYrxljThhj+oANwMo43peiKFHweL0A2G3hhcVmtXuTWB7Sa9WmDJ6iTQStWZn6xCM+m4A6KwstG18CwfqQPuuBm6zH1wMvG9///npgrZWpVgvUARsj2bSuecWygWXzqRhjPAcsE5F8S5Q+BuyO/xYoihIOt6Uq9ghejcP6xPckUX384hPsedlEkiqISnzEDLtZ6yu34vuQtwM/M8bsEpG7gM3GmPXAg8DDItKAzxtZa127S0QewycGbuAWfzJAOJvWkLcD60Tk28BWyzZRxmgXke/hEzQDbDDGPHNSd0VRFNweg02GPJxQbCkgPn4HJ3hNSmRIlJTUJZ41H4wxG/CFs4LbvhH0eAC4IcK1dwN3x2PTam/Elw0X2h5tjF/iS7dWFGWccHsNDlvk4IjfI/Ik8YM+UtgtGXrocrloampiYGAg8YOPM7m5uVRVVZGVlTVhY8QlPoqiTD48Xi8Oe+REAnsKeD7eMJ6PTUjKmk9TUxNFRUXU1NQkLQFjPDDG0NraSlNTE7W1tRM2jpbXURQlLG6viZhsAEPi402i+PiFT0I8n2QI4sDAAGVlZWktPOBbPysrK5twD07FR1GUsHi8JpBUEA6/+LiTuuZjhs0FfGtRyZpSuguPn0S8DxUfRVHC4vIY7FHWfPyhrmQu7qdS2E0ZHSo+Ch6vobk7/RdJlfHF4/VG9XxSKdV6ZMKBik8kXn31Va666ioABgcHufjii1mxYgWPPvpoQuehCQcK33n2fX78x0Y2f+1iphXmJHs6Sorg9pqoCQepkGodEBnd5zMmtm7disvlYtu2bQkfWz0fhRffOw5AR58zyTNRUgm3J0bCgaSA+HiHzwUm9z6fAwcOsGDBAm666SaWLVvG9ddfT19fH88++ywLFizg3HPP5cknnwSgubmZG2+8kW3btrFixQr27duX0Lmq56MoSljcXi/Z9sjfT/1eUTL3+bgt9QkOD/rK6yRrRj6+9ftd7D7aNa42F82awj9fvThmvz179vDggw9yzjnn8MUvfpHvfe97/PjHP+bll19m3rx5fOpTnwKgoqKCn/70p3z3u9/l6aefHte5xoN6PoqihMXpNjiiiE8g4SAFPB+bbXjCQTK9sWRTXV3NOeecA8CNN97I5s2bqa2tpa6uDhHhxhtvTPIMfajno4woUa4oAC6Pl+x4Npkm0c3wBFKth9p8qdbJ/a2Ox0OZKELTpDs7O1MyBVw9H0VRwuL2esmKw/Nxe5K/ydQmqRV2SyaHDh3izTffBOCRRx7h4osvZv/+/YE1nUceeSSZ0wug4qMoSlhcbhNVfBy2VNjnE2aT6SROOABYuHAhDz30EMuWLaOtrY2vfOUrPPDAA1x55ZWce+65zJkzJ9lTBDTspihKBJweL0VZkT8iUqG2myfMsQ+TfZ+PzWbj/vvvH9a2evVq3n///RF9zz//fM4///wEzWw46vkoihKWWNluthTwfAJhN1twqrXu80kHVHwURQlLvGG3ZK75BMJuWl4HgJqaGnbu3JnsacSFio+iKGFxebxkOWInHCQ1280bbs0neZ5PpoheIt6Hio8SROqlYyrJw+nxkhXXkQqJmtFIwotPctahcnNzaW1tTXsB8p/nk5ubO6HjaMKBoihhcXuih938L6Wc55OkfT5VVVU0NTXR0tKS8LHHG/9JphOJio+iKGHxhd2ieT4+9fEk0fXxmNTZ55OVlTWhJ39mGhp2UxQlLE6PF0eU83yGCosmakYjCRQW1X0+aYeKj6IoYXF5vGRHSzjwh92Suc8nXHmdSb7PJ11Q8VEUJSy+NZ9oh8n5w27JLCw6Muym+3zSAxUfRVFG4PUa3N40TTiYxPt80gkVH0VRRuCyFlPiKSyazCMVIiUcTOYjFdIFFR9FUUbgsqoWRAu7pUJtN2/EVOtkzUiJFxUfRQ/0UUbg9sT2fFJBfPyejyMk7AYaekt1VHwURRmBMw7x8b/mSuY+nzCFRQPhQNWelCYu8RGR1SKyR0QaROSOMK/niMij1utvi0hN0Gt3Wu17ROSyWDZFpNaysdeymR1tDBGpEZF+Edlm/QyvJa4oyqiJJ+yWCoVFwx+p4PtX061Tm5jiIyJ24D7gcmAR8GkRWRTS7UtAuzFmHnAvcI917SJgLbAYWA38UETsMWzeA9xrjKkD2i3bEcew2GeMWWH9/PWo7oCiKCNwuWN7Pg6/55PEXabhst0k4Pmo+KQy8Xg+q4AGY0yjMcYJrAPWhPRZAzxkPX4cuEh8vwFrgHXGmEFjzH6gwbIX1qZ1zYWWDSyb18YYQzlZ9C4qIbjjyHbze0XuZCYcmMhhN9We1CYe8akEDgc9b7LawvYxxriBTqAsyrWR2suADstG6FiRxgCoFZGtIvKaiHw0jvekKEoUnO54wm6+jw93Uj0f378adks/4iksGu63L/R/NVKfSO3hRC9a/2hjHANmG2NaReR04HcistgY0zVsgiI3AzcDzJ49O4wpRVH8uOJKOBCrb/Kz3Wwh5XUguVl4Smzi8XyagOqg51XA0Uh9RMQBFANtUa6N1H4CKLFshI4VdgwrpNcKYIzZAuwD5oe+CWPMA8aYemNMfXl5eRxvexKhf6NKCPGIj4hgt0kgRJcMvOESDmya7ZYOxCM+m4A6KwstG18CwfqQPuuBm6zH1wMvG1+S/XpgrZWpVgvUARsj2bSuecWygWXzqWhjiEi5lcCAiMy1xmiM/xYoihKK00o4yIlSWBR8GW9J9XwilNcB3eeT6sQMuxlj3CJyK/AcYAd+ZozZJSJ3AZuNMeuBB4GHRaQBn8ez1rp2l4g8BuwG3MAtxhgPQDib1pC3A+tE5NvAVss2kcYAzgPuEhE34AH+2hjTNvZboijKoF98suxR+2XbbUnNdouWcKCeT2oT12FyxpgNwIaQtm8EPR4Abohw7d3A3fHYtNob8WXDhbaHHcMY8wTwRMw3oShK3Ay6PUAcno9dkrrPx59pF67Cga75pDZa4UBRlBEMxht2s9uSuubjCXOkgv+EVc12S21UfBRFGcGgK76wW1aS13zCFRYNVF5QzyelUfFRFGUE8YfdbMnd52NGZrsFCp4mURSV2Kj4KAG0XoTiJ/6wm+BKZlVrr0FkeMKBPeD5JE8Uldio+CiKMoIh8YkVdkuu5+PyGLJswz/G7IF9Pur5pDIqPoqijGDQ5UEkenkdgCxHctd83B4vjpA5hq75GGOSmg6uhEfFR1GUEQy6veQ4bIEK0ZFw2JK7z8ftNcPSrCEo7GaJ4nee20PdV/+Q8Lkp0VHxURRlBD7xiR5yA59nlMx9Pi6Pd0QJoNATVn/06j5gKIliovnXP7zH9174ICFjpTMqPoqijGDQ7YmZbAA+zyeZC/tujxkRdrNHSLXuHZx48fF6DT9+rZH/fGmvlveJgYqPoigjGHR5ycmKQ3zsyV3zcXm9gaMd/DgibDLtHXQz0Rzp6A88bukZnPDx0hkVH0WLWisjiD/slnzPJzQpInTNx0+fc+I9n9ZeZ+Bxc5eKTzRUfBRFGUH8Ybdk13bzBo7z9hO65uOn1znxnk+wdxUsRMpIVHwURRmBP9stFllJrmrt8kTJdgvxyPoT4Pl0DwyJz4lu9XyioeKjBND1UcXPoCv+bLdk7/MJzXZzRNhkmgjxGe75qPhEQ8VHUZQRDLo9cSYcJLfCgdsbJdstRBT7XRMvPj1B4tPR55rw8dIZFR8lCHV9FB/xh92SW9vN5fGOKK/jFyP/mk+25RklUnwKcxx0Daj4REPFR1GUEQy4PGTHEXZzJLm2m9tjhh2nAEMVrv37fPzZcIkIu/UMusmyC9MKs+nqn/gEh3RGxUcJoGs+ip8+p4f8GGf5QPJPMnVFCbv513yyHAn0fAbcFOY4mJKXpZ5PDFR8FEUZQb/LQ152fPt8nMk8z8cbLuHA99wviv4EhEQlHBTmOpiSm0VXv4pPNFR8lADq+Ch++p0e8uMQnxyHL9U6WaVk3OFSrUPWfPxLUonwfLoH3RRkOyjOy6JTxScqKj6KogzD6fbi9pq4xcdrkndkddjCoiFrPv59SAlZ8xlwU5TrYEqeg64BXfOJhoqPoijD8H9I52U7YvbNttZT/IfPJZpoqdb+I7b9HlAiPJ9ep5uCHA27xYOKjxJAEw4UgD6X7xt7fJ6Pr89gAj7Yw+ELu4XfZOqxPB53AsUnOOFg0O1lIEn3JR1Q8VEUZRj+Apzxht0geZ6PL+wWfs3HLzoBzydBqdZFuQ6m5Pq8xm4NvUVExUcJYDTlQCEo7BZHqrW/CoIzlcJuMpRwYIxJuPgUZPs8H0DTraOg4qMoyjCGPJ/Yaz6BsFsSPZ/QsFvwmk9wIsREh908XkOf0+NLtbbERzPeIqPiowTQNR8FoM86eiCefT7+0jWJOqI6lEH3yEPvhtZ8zLBjFSZ6/cV/ZENhji/VGlR8oqHioyjKMPpHs+aTlbw1H2MMzjCH3gUfox3s+Uz0YXI9A2HER4uLRiQu8RGR1SKyR0QaROSOMK/niMij1utvi0hN0Gt3Wu17ROSyWDZFpNaysdeymR1rDOv12SLSIyL/MNqboPhQz0eB0SYc+LPdEi8+/soKoQVQRQSb+MJgHk/iwm7+4xQKcx2UqOcTk5jiIyJ24D7gcmAR8GkRWRTS7UtAuzFmHnAvcI917SJgLbAYWA38UETsMWzeA9xrjKkD2i3bEccI4l7gD/G+cUVRwtPn8u/ziT/bzelJfNjNn+QQrvq2w2az1nx8fbLswsAEez7dlvgU5OiaTzzE4/msAhqMMY3GGCewDlgT0mcN8JD1+HHgIhERq32dMWbQGLMfaLDshbVpXXOhZQPL5rUxxkBErgUagV3xv3XFT7JKoyipSb/Tv88njoQDf9gtCZ6PP9SXHUZ87DbBExR2K8rNos/lmdDfdb/nU5TjIMtuoyDbrmf6RCEe8akEDgc9b7LawvYxxriBTqAsyrWR2suADstG6FhhxxCRAuB24FvR3oSI3Cwim0Vkc0tLS4y3PDnRVGsFhsJu8aRaDyUcJCHsFtXzEVweb0B8CnMceLxmQk9d9a/5FOT4RLskP1s9nyjEIz4Spi30fzBSn/FqjzbGt/CF6XrCvD7U0ZgHjDH1xpj68vLyaF0VZVLT5/SQm2UbcU5OOHKy/KnWiQ+7DQbEZ6RIOuwybM2nyNr0OZHrPt1BB8kBTNHiolGJ7Vf7vI/qoOdVwNEIfZpExAEUA20xrg3XfgIoERGH5d0E9480xpnA9SLyHaAE8IrIgDHmv+N4b0oQGn1TALoHXBTlZsXVN7DmkwTPxy944cJuWXabVSDVNy+/IAy4PIFMtPEmEHazhK44z0Fnv3NCxsoE4vF8NgF1VhZaNr4EgvUhfdYDN1mPrwdeNr7g6npgrZWpVgvUARsj2bSuecWygWXzqWhjGGM+aoypMcbUAN8H/p8Kj6KMnS6rMnM8JLO8TrSwm/+cIY93uOczkenWoWE3PVYhOjF/w4wxbhG5FXgOsAM/M8bsEpG7gM3GmPXAg8DDItKAzxtZa127S0QeA3YDbuAWY4wHIJxNa8jbgXUi8m1gq2WbSGMoijK+dA+44/Z8klnVOlrYzXfO0NAaj//9TGSJnR6nm2yHLXDEQ0leNp39HRM2XroT19cbY8wGYENI2zeCHg8AN0S49m7g7nhsWu2N+LLhQtsjjhHU55vRXlcUJTbdA65AYcxYBBIOklC92Rkl2y3LbsPlHun5TOiaz4B72H0rzs/SbLcoaIUDJYCu+Sjg93ziEx8RIcdhS5Ln4xOSsGE3h+D0jFzzmVDPxzpOwU+xHqsQFRUfRVGG0T3gGvYhGotkiU9MzydozacwAZ5Pz6A7MA4QSGzQQ+XCo+KjKMowekax5gOQm2VPyHEFoQxGSTjIDmS7WWG3nASIz4Cbopyh++YXnw4Vn7Co+CgBdJOp4vEaep2euMNu4Mvu6ktCaMlfVSEnzGbYbIfP83GHJBxMZImd7giej2a8hUfFR1GUAP504dF4PvnZdvoGE39iZ9Q1HyvVOnTNx39cxETQM+gKeFgAJfmW56NJB2FR8VECaMKB4j95czSeT362fcKPKwhHtOrb2XYbLrcJs+YzcWtTPQPDPZ+p+dkAtPUOTtiY6YyKj6IoAbotzyfeVGvwFSCdSI8iEr1ODyKQG2afT5Y/7GaJT0H2xK75GGN8CQdBnk95UQ4AJ3q0ykE4VHyUAOr4KEOezyjDbknwfPqdbvKy7NjC1KDLsguDQft8shxCXpY9ULF7vBl0e3F5zDDPJzfLTlGOg5Zu9XzCoeKjKEqAjj7ft3R/yCgefJ5P4sWn1+mJeOBddkiqtcMm5E2gSPo9xqKQFPVpRTm09Kj4hEPFRwmg5/oobb0+z2dqwWg9n8SH3fqdnohnDvmz3fziY7fZJrTWWrflMRaGhCunFWZzQj2fsKj4KIoSoH0snk9OcsJuvYPuiJ6Pb5Pp0GFyDptQMoHlbtotu6H3rbwohxPq+YRFxUcJoH6P0t7rJC/LTm4cB8n5yc9yDFtfSRT9rshhN/+RCh4r1dpuE6bmZwfEdbxp7w0v2tMKczThIAIqPoqiBGjrc1JaEL/XA1CQ4xOARIfefJ5PhLCb3VfbzV/V2uEXn96JEYI2S9RC7920whw6+11JOWwv1VHxURQlQEefa1TrPQB52X7xSewHbF+0hIOQox58nk9WIDw23gQ8n4KRYTeAVvV+RqDiowTQfAOlrdc5qvUeGNpDk2jxiRV2AwKp1Q6bjakF2fS7PBNSZbq9z0W23UZByHwqLPE53jUw7mOmOyo+iqIEaO8bvfj4PZ/eBJfY6R30kBch7BYQH0to7HYJvK+JSDpo73VSkp+FyPA9R7NK8gA40tE/7mOmOyo+ShDq+kx22nudTM0fXdjNv6s/0eLTNeBiSl7kVGsY8sYcVtgNfN7deNPaOxh2raxyqiU+7So+oaj4KIoC+Ap1dg24KS3IGdV1U6xqCF0DiROfAZcHp9sbGDuUvKzh3pjdJoH1l+bu8Q+BHe0YYGZx7oj2KblZFOU61PMJg4qPEkDXfCY3/jIw06eMTnyScXSAfyz/2KH4Q4E9fvERYaYVAvuwc+zi09nvCrtm9GHXQMB+KJUleer5hEHFR1EUAJoD4jPyG3w0kiE+XbHEx/J8ugfc2ARsNqGiKAcRODpG8Wlo7mHV3S/yiR++gdszVB17wOWhrdfJrDCeD0DV1Dz1fMKg4qMEUMdnctNsZWT5w1Px4i8pk0qeT26Q+Dhsvo+5LLuN8sIcPuwcmxD8busRBt1edh/r4pkdxwLtRy1hmVmsns9oUPFRFAWA411j83zsNqEo1xHwRhLBaMJu9qCq1zNL8jg2Rs9n+5FOFs+aQk1ZPg+9cSDQvq+lF4Da8oKw11VNzad70D1hG1zTFRUfRT0eBfDtRbHbhLJRVjgAnwgkQ3ymxAi79Qy4cQSLz5TcgKcyWpra+5hTls+NZ83hnUMd7DzSCcAHx7sBqKsoDHvdKRU+UdrX0jOmcTMVFR8lgCYcTG6auwcpL8wJez5OLCayYnQ4Yno+WUGej33o/cwpy+dwW/+wNZt4MMZwpL2f6qn53HB6NblZNh5+8yAAu492UVmSF/EMpFPKfaKk4jMcFR9FUQCf5zPaTDc/U3ITKz7tvU5EIp+4mpvt+2jrGRzu+ZxSUYjT4+XwKNdgWnoGGXR7qZqaR3F+FteuqOSpd4/Q3uvkzcZWzqwtjXht1dR8sh22QHhO8aHiowTQ83wmN8c6B0a93uMn0Z5PS88gZQU5OOzhP8LygqpyB6/5BLyQ5tF5IYfbfGJVNTUfgM99ZA4DLi+f//lG2nqdXLiwIuK1dpswd1oBDaMcM9NR8VEUBWMMTe19VJfmj+n6hItP92DUrLzgIyFyHEOP51nrMg2jDIE1tfcBvrRpgMWzirl6+SzebepkwYwiLls8I+r1p1QUatgthLjER0RWi8geEWkQkTvCvJ4jIo9ar78tIjVBr91pte8Rkcti2RSRWsvGXstmdrQxRGSViGyzft4VkY+P9WZMdtTvmby09Awy4PJSPTV8unAsygqzaet14k3QmT6xxCfLbiPLWuvxl9oBn0hOn5LD+8e6RjVekxWmqwy6P9//1Aqe+PJHePzLZwdqyUWirqKQQ219CS9BlMrEFB8RsQP3AZcDi4BPi8iikG5fAtqNMfOAe4F7rGsXAWuBxcBq4IciYo9h8x7gXmNMHdBu2Y44BrATqDfGrLDG+LGIhA8EK4oSFn9YaayeT3lRDm6voSNB3s+JHiflhdHXp/zeT45j+Mfc8qoSth3uGNV4Te39lBVkDzs/yG4TTp9TGqhtF42llcUYA7tHKXqZTDyezyqgwRjTaIxxAuuANSF91gAPWY8fBy4SX3nXNcA6Y8ygMWY/0GDZC2vTuuZCywaWzWujjWGM6TPG+L9O5KJf4MeMLvlMXvxhpZMRHxgq0TORGGNiej4wtO6THSI+K+dM5UBr36gKjDa19wVCbmNhaWUxANubOsdsI9OIR3wqgcNBz5ustrB9LCHoBMqiXBupvQzoCBKT4LEijYGInCkiu4AdwF8HXR9ARG4Wkc0isrmlpSWOt60okwd/WGmsH7B+L+REz8SLT1uvE6fHGzgrJxJ+jyQ7JCR2WnUJAO8cbI97zKb2/kCywViomJLL9Ck5gb1BkegecHHFD17nF3/eP+ax0oV4xCdc0n/od+RIfcarPeo8jDFvG2MWA2cAd4rIiJQdY8wDxph6Y0x9eXl5GFOKMnk53NbHtMLsiMdSxyKRno8/TXp2DC/NvwE1J2v4AW/Lq0vIy7Lz+t74voR6vb49PlWlY/d8wOf9bG+KHu57fe8Jdh/r4pu/353x2afxiE8TUB30vAo4GqmPtd5SDLRFuTZS+wmgJGjNJnisSGMEMMa8B/QCS+J4X0oIRiOWk5Z9LT3MnRZ+h348JFJ8DrXFFyIssc7vCfV8crPsnDNvGi++1xzXB3xLzyBOj/ekPB+ApZUlNJ7opWsg8rpYcPWFpgyvBxeP+GwC6qwstGx8CQTrQ/qsB26yHl8PvGx8/6vrgbVWplotUAdsjGTTuuYVywaWzaeijWHZcACIyBzgVOBA3HdAUSY5xhg+ON7DvOljF5/CHAe5WTZaEhB2OxwQn+ieSEnA8xn5MXfRwgqOdPTz3rHumOOFplmPlTNqp2IMbD7QFrFPcN25PR/Gnls6E1N8rPWTW4HngPeAx4wxu0TkLhG5xur2IFAmIg3AbcAd1rW7gMeA3cCzwC3GGE8km5at24HbLFtllu2IYwDnAu+KyDbgt8DfGGNOjO12THLU8ZmUtPQM0tnvilibLB5EhIqi3JM6Kyde4g0RlljHZheFyUa7bPEMsuzC41ua4hjPygQ8SfFZOXsq2XYbbzVGE5/+wFrWnuOZLT5xBXiNMRuADSFt3wh6PADcEOHau4G747FptTfiy4YLbQ87hjHmYeDhmG9CUZSwNBz3bX6sqyg6KTvVpXkctryEiaShuYfaaeErSAfjr/tWEEZ8SguyuXTxDJ7c2sQ/rj512KbUUPxhvpMNu+Vm2VleXczbja0R+xzrHKBueiFZdlugYGmmohUOlADq+ExO9lplX+pOIuwGvgQAf0hsovB6De9/2M3CmVNi9i2y6r71OUeePArw2TNn09Hn4jebD4d93c+htj5mTMmNKlDxctbcMnYe7QqcsBrKsY4BZhbnMX96IXuPZ3ZFBBUfRff3THLe/7CbKbmOmKnLsaiams+JHueE7uI/0tFPz6CbBTNii88586YBcNbc8EU/PzK3jPo5U7nvlX1hj8b2c6itL2ZmXbycNbcMj9ewcf9I78ft8dLcPcCs4lzqphexr6UHT4IqRiQDFR8lgIrQ5GR7UwfLqkrw7fEeO/4P6IkMve066qsQsGBm7BDhwplTeOfrl7BmRei2RB8iwm2XzufDrgF+/FpjRDsHTvSOefNtKKfPmUpelp1X94xM827uHsRrYEZxHvMqChl0ewPJDpmIio+iKdaTmAGXhz0fdmRsk9IAACAASURBVLOsqvikbfk/oA+2TtwH5qYDbWQ7bCyKI+wGvrWdaJx9yjSuWT6L+15poKF55BrLiZ5BmrsHWRiH2MWDL827jJffH5nm7c90m1mSG0j+yOTQm4qPEkBFaPKx+1gXbq9hWVXJSduaax0jPZFHB7y9v5XTqkvGZf3Fz9evWkRBjp1bf72VPufwkOFuy9NaNCs+sYuHCxZU0NTeP6LK9bFOX1bdzOLcQPXtvRl8DIOKj6LhtknMdqvA5vLqk/d8puRmUTU1b8KKZ3b0Odl9tIsz55aNq93yohx+sPY09hzv5vYndgyrzL1xfxt2m7Ck8uTvj5/zT/Wd/fPK+8NDb/409ZnFvlNRZxbnsjeMN5YpqPgoKj6TmE0H25lZnMuMMR4iF8rCmVNGfVxBvDy/+zheAxdHObhtrJw3v5x/vGwBv3/3KN9YvxOv12CM4YXdx1k5u4QpEY7IHguVJXmcOr2Il94/Pqz9aMcA+dn2wOms8yoyO+NNjx5QAqgITS68XsNb+1r52KnlJ51s4GfhzCm89N5x+p0e8rLHLzQGsGHHMSpL8gIVosebv/7YXDr7Xdz/2j6OtPezonoqe45382+fWDruY126eDr3vdIwrDr3sc5+ZhbnBv4v6iqKeGTjIbxeg802Pv8/qYR6PooySfmguZvWXidnnzJt3GwuryrGaxj1eTmxONzWx2sftHDtabPGTShDERFuX30q37pmMW/sa+XeFz/gI3PLuO70qnEf68plM/EaeHbXh4G2Ix39VAZtZK2bXki/y8ORjsys8aaejxLIulHHZ3Lx5j7fXpOPnDJ+ayj1NaXYBN5sbB1Xuz//8wHsInzurJpxsxkOEeGms2u49rRKDrf1sXDmFOwT4HWcOr2IU8oLeGb7UT531hwAjrT3szgoscGf8dbQ3DNuqd6phHo+iorOJOXl95uZO62AypKTq1kWTHFeFotnFfPWvsglZEbLodY+fvnWQa49rZIZxeOzNhWL4rwsllQWT4jwgE/krlw2i43722juHqDf6aG11zns/2Io4y0zkw5UfJQAmX5+iDJEZ7+LN/e1csni6eNu+6N109hyqJ3Wcahw7fUavvbUTuw24f9eduo4zC51uNoKva3fdpQDrb3A8GMiSvKzKS/KydikAxUfRRMNJiGv7mnG7TVctnjGuNu+atksPF7Dhp0fxu4cgx+9to8/ftDCP125kOnjlJGXKtRNL+L0OVP55VsHAyecLg7ZT1RXUZixe31UfJQAqkGThw07jlFRlMOKcdhcGsrCmUXMn17Io5sOnZQ3/cjGQ/z7c3u4evksbjxz9jjOMHW46ewaDrT2cceTO5iS66A25EC/uopC9h7vHrb3KFNQ8VG0ssEko7VnkJfea2bNilkTksIrInzp3Fp2HunitQ/iO6o6GK/X8N3n9nDnkzv42PxyvnvDsgnLcEs2VyyZwfKqYjxew5oVlSPWmBbPKqbX6WG/FZbLJDTbTdGw2yTjt1uP4PYabqivjt15jHz8tCr++5UGvvX73ZxZWxb3np+G5m7ueGIHmw+286n6av7l2iVkOzL3O7LDbuOhL67ijX2tXLhg5ObZZVblie1NHZxSfnJHXqQamfu/qoweFaGMxxjDo5sOs6K6hPnTx6dYZjiyHTbu+cQy9p/o5e8e3YrT7Y3a/1BrH3c8sZ3V33+dvc09fPeG5fzbdUszWnj8lORnc8XSmWHr1c0rLyQ3y8b2ps4kzGxiUc9HUc2ZRLz6QUvgw32iOXveNP756kV86/e7+cSP/sxtl8zn7FOmkZtlx+3xcqC1j7caW3lm+zHebGwl22Hjs2fO5tYL6wK7/ic7DruNJbOKVXyUzEbXfjKfH7+2j5nFuVyzfFZCxvvCObXMLM7ln9fv4ou/2IxNoCDbQa/TjX8NfU5ZPrddMp9PnVGdcRlt48GyqhJ+vfEgbo8Xhz1zPEEVH0XXfCYJWw628VZjG1+7cmFCw1mrl8zkwgXTeX1vC9ubOukacFGU46CqNJ9VNaXMKcvP2ISC8WB5dTE/+7OXD473jOvRDslGxUfBH3hTEcpcjDHc/cx7VBTl8JkkpC1nO2xctHA6Fy0c/02tmY7/rKV3mzoySnwyx4dTxoyKTubz7M4PeedQB7ddMp/8bP3OmU7UlOVTWpDNloPtyZ7KuKLiowRQEcpMugdc/MvTuzl1etGEplcrE4OIcPqcqWw+0JbsqYwrKj6KphlkOPc8+z7Hugb41+uWTlihTGViOaNmKgda+2juHkj2VMYNFR9FC4pmMH/ae4JfvnWIL55Ty8rZU5M9HWWM1NeUArDlQOaE3lR8lAAqQZnFh50D/J91W6mrKOTvL52f7OkoJ8GSWcXkOGxszqB1HxUfRUUnA3G6vfztI+/Q7/LwoxtXapJBmpPtsLGiuiSj1n3iEh8RWS0ie0SkQUTuCPN6jog8ar3+tojUBL12p9W+R0Qui2VTRGotG3stm9nRxhCRS0Rki4jssP69cKw3Y7Kj4bfMwBjD7U9sZ9OBdv7tumXMq5i4MjpK4jijppSdR7voc7qTPZVxIab4iIgduA+4HFgEfFpEFoV0+xLQboyZB9wL3GNduwhYCywGVgM/FBF7DJv3APcaY+qAdst2xDGAE8DVxpilwE3Aw6O7BYpqTmbxnef28NutR/iHS+cnrJKBMvGcXjMVj9ew7VBHsqcyLsTj+awCGowxjcYYJ7AOWBPSZw3wkPX4ceAi8W1ZXgOsM8YMGmP2Aw2WvbA2rWsutGxg2bw22hjGmK3GmKNW+y4gV0S0MNQYUA1Kf/7rpb386NV9fObM2dxywbxkT0cZR1bOnooIbMyQ0Fs84lMJHA563mS1he1jjHEDnUBZlGsjtZcBHZaN0LEijRHMdcBWY8yI83tF5GYR2Swim1taRn/GSCaj4bb0xxjD957fw3+88AGfOK2Su65ZrCVrMozivCxOnV7E5gzJeItHfML9Bod+WkXqM17tMechIovxheL+Kkw/jDEPGGPqjTH15eXl4bpMelSD0hOv11c65z9fbuCT9VX8+w3LM6oApTLEqtpS3jnUjtsT/YiKdCCe39AmIHhbdBVwNFIfEXEAxUBblGsjtZ8ASiwboWNFGgMRqQJ+C/yFMWZfHO9JCUI1J33pd3r48q+28NM/7efzZ9fwb59YphtJM5gzakrpc3rYfawr2VM5aeIRn01AnZWFlo0vgWB9SJ/1+Bb7Aa4HXja+WM56YK2VqVYL1AEbI9m0rnnFsoFl86loY4hICfAMcKcx5s+jefOKjyGPR2UonTjeNcCnHniT53cf55+vXsQ3r1k8IcdiK6nDGdZm043703/dJ6b4WOsrtwLPAe8BjxljdonIXSJyjdXtQaBMRBqA24A7rGt3AY8Bu4FngVuMMZ5INi1btwO3WbbKLNsRx7DszAO+LiLbrJ+R59EqSgbxp70nuOIHr9PQ3MNPPlfPF86pTfaUlAQwoziX6tI8NmVA0kFcO8+MMRuADSFt3wh6PADcEOHau4G747FptTfiy4YLbQ87hjHm28C3Y74JJSKacJA+eLyG/3p5Lz94aS/zygv50Y0rdR/PJOOMmlJe29OCMSatk0p0VVIJoBqU2hzp6OdzD77N91/cy8dXVPLUreeo8ExCVtWU0trrpPFEb7KnclJozQ1FV3pSHGMMv9ncxF1P78YYwz3XLeWT9dVp/a1XGTv+IqOb9rdxSnlhkmczdlR8lAAqQqlHc9cAdz65g5feb+asuaX8+/XLqS7NT/a0lCRySnkBZQXZbDzQxtpViT+VdrxQ8VE03JaCeLyGX719kH9/bg9Ot5dvXLWIz59do9lsCiJCfc3UtN9squKj4LXUR0UoNdjR1MlXf7eD7U2dnDtvGnetWczcNA6vKOPPGTWlPLfrOMe7Bpg+JTfZ0xkTKj5KQHyU5NI14OJ7z3/A/7x5gNKCHH6wdgXXLJ+lazvKCIL3+1ydpsVjVXwUvJb2GF31SQpuj5dHNh3m3hc+oL3PyefOmsPfX3oqxXlZyZ6akqIsnjWF/Gw7mw6o+ChpjMeropMsXt3TzN3PvMfe5h5W1ZbyjasWsaSyONnTUlIch93GytlT2ZTG6z4qPoqSBD443s3dz7zHax+0MKcsn/tvPJ3LFk/XEJsSN/U1U/nBS3vp7HelpZes4qME0KWfiaepvY8fvLiXJ95poiDHwdeuXMjnPjKHHIc92VNT0oxVNaUYA+8cbOeCBelXUUzFR1ESQHP3APe93MCvNx5CEG46u4a/vbCO0oLsZE9NSVNOmz0Vh03YeKBNxUdJb9TxGX86+pzc/1ojv3hjPy6P4ZP1VfzthXXMKslL9tSUNCcv286SymI2p2mRURUfRZkAugZc/OLPB/jJHxvpcbq5ZvksvnLxfGqmFSR7akoGcUbNVB564yADLg+5WekVulXxUQJodeuTp73Xyc/+vJ9fvHGA7gE3lyyazt9fOp8FM6Yke2pKBnJGTSk/eX0/25s6WVVbmuzpjAoVH0UZB1q6B/np6408/NZB+pweLl8yg1sumKdp08qEEigyeqBNxUdRJhMfdg5w/2v7eGTjIVweL1cvn8UtF8xj/nQ96kCZeEoLsplXUZiWh8up+CjKGNh/opefvN7I45ub8BrDx0+r5Mvnn6I12JSEc0ZNKU9vP4rXa9Kq8KyKj6KMgi0H23ngj/t4fvdxsmw2rq+v4ssfO0WPOVCSxmmzS3hk4yEaT/Sk1eGCKj5KAM03CI/Xa3jhveP85I+NbD7YTnFeFrecP4+/OHsOFUXpWVFYyRxOqy4BYNvhThUfRckEBlwennznCD99vZHGE71UTc3jm1cv4ob6agpy9E9HSQ1OKS+kMMfBtsPtXH96VbKnEzf6F6QE0KrWPlp7BvnV24f4nzcPcKLHydLKYv7r06dx+ZIZOOy2ZE9PUYZhswnLqorZdrgj2VMZFSo+GcJYN5np3p4hdh7p5BdvHGD9u0dxur2cf2o5N583l4/MLdOCn0pKs6K6hAf+2JhWm01VfDKADTuO8Te/eofn/u48Tp0xuphv8HEKk1GH3B4vz+06zi/e2M+mA+3kZ9v5VH01N509J63i58rkZkV1CW6vYeeRzsDen1RHxScDePG94wDsONI5avFxerwTMaWUp63XySMbD/HLtw5yrHOA6tI8vnblQm6or07L8vTK5GbFbH/SQYeKj5IeuNyTy/PZfbSLh944wO+2HWHQ7eXcedP4lzVLuGBBBfY02iOhKMFUFOVSWZKXVus+Kj4ZgMvjU40s++g/PAc9nvGeTsox4PLwh53H+NVbh9h8sJ28LDvXn17FTWfXaCUCJWNYXp1eSQcqPhmA0+0TkBzH6DOx/MKViTS29PDIxkP8ZksTHX0uaqcV8NUrFvLJ+mqK8zW0pmQWK6pL2LDjQ070DDKtMCfZ04lJXJ9WIrJaRPaISIOI3BHm9RwRedR6/W0RqQl67U6rfY+IXBbLpojUWjb2Wjazo40hImUi8oqI9IjIf4/1RqQzg27fuo3DNnrxcbqH1nwyQYZcHi8bdhzjsz99iwv/4zV+/ucDnHPKNH79v87k5b//GH953lwVHiUjWVE9FYB308T7ien5iIgduA+4BGgCNonIemPM7qBuXwLajTHzRGQtcA/wKRFZBKwFFgOzgBdFZL51TSSb9wD3GmPWicj9lu0fRRoDGAC+DiyxfiYdLitpwDuGRRtXhiQcNLX3sW7jYR7dfJiW7kEqS/L4v5edyg31VVqFQJkULK0sxm4Tth3u4KKF05M9nZjEE3ZbBTQYYxoBRGQdsAYIFp81wDetx48D/y2+jRFrgHXGmEFgv4g0WPYIZ1NE3gMuBD5j9XnIsvujSGMYY3qBP4nIvFG874zC772MRXyGeT5plnHg9nh5dU8Lv954iFf2NCPAhQsq+OyZczhvfrkmECiTirxsO3UVhWxv6kz2VOIiHvGpBA4HPW8CzozUxxjjFpFOoMxqfyvk2krrcTibZUCHMcYdpn+kMU7E8R4yGr+AjMWJGXSnX8LB/hO9PLrpME+800RL9yDlRTncesE81q6aTaUeT61MYpZUFvPK+80YY1J+Y3Q84hPuHYR+RY7UJ1J7uMWJaP3jnUdERORm4GaA2bNnx3tZWjB4Ep5P94A78DiV/Z4+p5sNOz7ksU2H2XigDbtNuODUcj5ZX80FCyrI0rI3isLSymIe39LEh10DzCxO7S9i8YhPE1Ad9LwKOBqhT5OIOIBioC3GteHaTwAlIuKwvJ/g/pHGiAtjzAPAAwD19fWp/Dk7ak4m7NY7mLqejzGGd5s6eXTTYX7/7lF6Bt3UTivgH1efynUrq5g+RddyFCUY/8m5O5o6M0J8NgF1IlILHMGXQPCZkD7rgZuAN4HrgZeNMUZE1gO/FpHv4Us4qAM24vNiRti0rnnFsrHOsvlUtDHG9rYziz6nT0CCS+XES+/gkOeTKq5PW6+T3249wmObDrPneDe5WTauWDqTT9VXs6q2NOXDCYqSLBbNnIJNYOfRLi5dPCPZ04lKTPGx1lduBZ4D7MDPjDG7ROQuYLMxZj3wIPCwlVDQhk9MsPo9hi85wQ3cYozxAISzaQ15O7BORL4NbLVsE2kMy9YBYAqQLSLXApeGZONlNL1On4CMRXx6gsUniXi8htf3tvCbzU08v/tDXB7D8uoS7v74Eq5ePospuZoerSixyMu2M6+ikJ1HUj/pIK5NpsaYDcCGkLZvBD0eAG6IcO3dwN3x2LTaGxnKiAtujzZGTdQ3kOH4vZex+IHD13wS7/ocau3j8S2HeXxLE0c7B5ian8Xnzqrhk2dUsWDGlITPR1HSnSWzinm9IfXzsLTCQQbgd3g8Y1Cflp6BcZ5NbPzJA7/ZfJi397chAh+tK+erVy7i4kUV5DjSoyS8oqQiSyqLeXLrEZq7BqhI4XVRFZ80JzjU5h5D2K25a5CiHAfdExx+M8aw+WA7v9l8mGe2H6PX6aGmLJ9/uHQ+n1hZxSxNkVaUcWFplZV0cKSTi1R8lImiuXvIc/GORXy6B6mYkkN3i3tCqlof6+znyXeO8PiWJvaf6CU/286VS2dyQ301Z9RM1eQBRRlnFs2cggjsPNKV0pUOVHzSnKMdQ+IzWs/HGMOB1l5WVJewr6V33OY04PLwwu7j/GZLE3/a24LXwJm1pdxywTwuXzKDghz9tVOUiaIgx8HcaQXsSPGkA/0USHMOnBgSDY93dCUOWnoG6ehzceqMIl7d08IYHKdh7Gjq5LHNh3lq2xG6BtxUluRx6wXzuO70KuaUFZyccUVR4mZJZTFvN8a9DTIpqPikOTuOdCLiy3Qbrefz7mHfN6Ol1sa0sSQs9Ay6eWrbER7ZeIidR7rIcdhYvWQGN5xezdmnlGHT+mqKknCWVhbz1LajgfJTqYiKTxpjjOG1D1pYVVPK2/vb8IzybJ7XPmgmP9vOKuvYXfcoisPtaOrk1xsPsX7bEXqdHhbMKOKuNYtZs6JSj6FWlCSzeJbvC+XOo51ccGpFkmcTHhWfNOb1vSfYf6KXm8+by9v720bl+fQ7PazfdpQLFlSQl+1LbY61SdXjNTy360N+8nojWw91kJtl46pls/jMmbM5rbpEkwcUJUVYXOnbI7ezScVHGWc6+px87Xc7qZqaxydWVvLV3+4YVYWD+1/bR9eAmy+cXRM4hC7SqaaDbg+PbTrMT17fz6G2PmaX5vONqxZx3elV6uUoSgoyJTeL2hRPOlDxSUPae518/heb+LBzgHV/dRY5DjsOmy1uz+f1vS3c90oD1yyfRX1NadCRDMPDbsYYfrftCP/x/Ac0tfezcnYJd16+gEsXz9CzchQlxVlSWcw7B9uTPY2IqPikGY0tPfzVw1s42NbHf3/mNFbO9h2da7dJXNluG3Yc47bHtjGvopC7P+47+NVhCUmweB1u6+P2J7bzxr5WFs+awv/7+FI+WjdNQ2uKkiYsrZzC7989Sluvk9KC7GRPZwQqPmnE7989yp1P7sBhF37xhTM4+5RpgdccNonq+bR0D/KdZ9/nN1uaWDm7hAf+op4iq1inzSbYBNxW2G3zgTb+8n824/IY7v74Ej59xmzNWlOUNGPJrKFKBx+bX57k2YxExScNONrRzzfX7+L53cdZObuE//rMyhEndtrtErbCgdPt5aE3DvCfL+1lwO3hr86by1cumU9u1vD6af6w3XvHurjpZxupmJLLzz5/BrXTdH+OoqQji60tFDtVfJTR0jXg4qev7+fB1xvxGMPtqxfwvz5aG/bUznCez6t7mrnr6d00tvRy4YIKvnblQuaWF4Ydy2EXBlwe/vaRrRTmOvj1X56Z8odRKYoSmeK8LOaU5bOjKTWTDlR8UpD2Xie/evsgP3l9P539Li5fMoN/umIh1aX5Ea/xrfn4xOdQax93Pb2LF99rpnZaAT///BlcsCB6uqXdJvx64yGcbi8//Yt6FR5FyQCWVBaz7VBHsqcRFhWfFOK9Y1089MYBfrv1CINuLxctqOArl8wPHI0bDX/Y7Hdbj3DHk9uxi3Dn5Qv4wjm1ZDtGekojrxe63V4WzCjiooWpuS9AUZTRsbSymGe2H6O918nUFEs6UPFJMi3dgzyz/Sjr3z3KO9bGzU+srOKms+eM6jA1u014c18rT7zTxKqaUr6/dsWovJf2PhcAn6yv1ow2RckQ/KWzdhzp5LwUW/dR8UkCXQMunt35Ib9/9yh/bjiB18CCGUX80xUL+GR9NSX5o/+G4rAJh9r6qCzJ4+dfOIP87NH9184tL6CxpZc1K2aNemxFUVKT4Iw3FZ9JyoDLw0vvNbP+3SO8sqcFp9vL7NJ8/ub8eVyzYhbzpxedlP2WnkEAPr2qetTCA3DPdcto6R6krDA1ixAqijJ6ivOzmF2az84UrHSg4jOBuDxe/tRwgt9vO8pzuz6k1+mhvCiHz545m2uWz2LFONZD81ccuHLZ2DyXM6ziooqiZBZLK4t5tyn1kg5UfMYZYwxbD3fwu61HeHr7Mdp6nUzJdXDVslmsWTGLM+eWTUhpmq9fuYgth9p1X46iKMNYUlnMMztSL+lAxWecaO0Z5JGNh/jNliYOtvaR47Bx8aLpXLuikvPmTyPHYY9t5CS47vQqrju9akLHUBQl/fAnHew82slH61Jn3UfF5yRp7hrg3hf38sQ7TTjdXs4+pYxbL5jH6iUzAuVrFEVRksUS63iFHUdUfDICYwy/fOsg//qH93F5vNxQX80Xz6lhXsXJJQ4oiqKMJyX52VSX5rHrSFeypzIMFZ8x4PUa/um3O1i36TAfm1/OXWsWM6dM11oURUlNllYWp9zZPrG3visj+P6LH7Bu02FuueAUfv75M1R4FEVJaZZUFnOorY+2XmeypxJAxWeUfHC8m/te3ccnVlbyD5eeqkcNKIqS8pxZ69tK8XZja5JnMoSKzyh58PX9ZNttfP3KRVqGRlGUtGBZVQn52XbeTDfxEZHVIrJHRBpE5I4wr+eIyKPW62+LSE3Qa3da7XtE5LJYNkWk1rKx17KZPdYxxhtjDC+9f5xLFk1PqXx5RVGUaGTZbdTXlPLGvjQSHxGxA/cBlwOLgE+LyKKQbl8C2o0x84B7gXusaxcBa4HFwGrghyJij2HzHuBeY0wd0G7ZHvUYo70R8XC8a5ATPU7qa6ZOhHlFUZQJ49x5ZTQ097D/RG+ypwLE5/msAhqMMY3GGCewDlgT0mcN8JD1+HHgIvHFpNYA64wxg8aY/UCDZS+sTeuaCy0bWDavHeMY405rr69+WkWR1j9TFCW9uHZFJXab8B/P70n2VID4Uq0rgcNBz5uAMyP1Mca4RaQTKLPa3wq5ttJ6HM5mGdBhjHGH6T+WMQKIyM3AzdbTHhFpBU5EfNdRuPyesVyV0kxjjPciA9F74UPvwxAZdS/uA+777JgunQbMGa95xCM+4VbVTZx9IrWH87ii9R/LGMMbjHkAeMD/XEQ2G2Pqw1w76dB7MYTeCx96H4bQe+HDug8142UvnrBbE1Ad9LwKOBqpj4g4gGKgLcq1kdpPACWWjdCxRjuGoiiKkqLEIz6bgDorCy0b3+L++pA+64GbrMfXAy8bY4zVvtbKVKsF6oCNkWxa17xi2cCy+dQYx1AURVFSlJhhN2t95VbgOcAO/MwYs0tE7gI2G2PWAw8CD4tIAz5vZK117S4ReQzYDbiBW4wxHoBwNq0hbwfWici3ga2WbcYyRgweiN1l0qD3Ygi9Fz70Pgyh98LHuN4H8TkPiqIoipI4tMKBoiiKknBUfBRFUZSEMynFJ1a5oExARH4mIs0isjOorVREXrBKF70gIlOtdhGR/7Tux3YRWRl0zU1W/70iclO4sVIZEakWkVdE5D0R2SUi/8dqn1T3QkRyRWSjiLxr3YdvWe0pW85qorGqrWwVkaet55PyXojIARHZISLbRGSz1Tbxfx/GmEn1gy/BYR8wF8gG3gUWJXteE/A+zwNWAjuD2r4D3GE9vgO4x3p8BfAHfHumzgLettpLgUbr36nW46nJfm+jvA8zgZXW4yLgA3wlnSbVvbDeT6H1OAt423p/jwFrrfb7gS9bj/8GuN96vBZ41Hq8yPqbyQFqrb8le7Lf3xjvyW3Ar4GnreeT8l4AB4BpIW0T/vcxGT2feMoFpT3GmD/iywoMJrhEUWjpov8xPt7Ct9dqJnAZ8IIxps0Y0w68gK9+XtpgjDlmjHnHetwNvIevAsakuhfW++mxnmZZP4YULmc1kYhIFXAl8FPreUqX9koCE/73MRnFJ1y5oBHleDKU6caYY+D7UAYqrPZI9ySj7pUVLjkN37f+SXcvrDDTNqAZ34fDPuIsZwUEl7NK6/tg8X3gHwGv9Tzu0l5k3r0wwPMiskV8ZcggAX8fk/EY7bjK8UwyTqp0UTogIoXAE8DfGWO6JPJZTBl7L4xv/9sKESkBfgssDNfN+jdj3mdFeAAAAa1JREFU74OIXAU0G2O2iMj5/uYwXTP+XlicY4w5KiIVwAsi8n6UvuN2Lyaj5zOZy/Ect1xkrH+brfbRlkFKK0QkC5/w/MoY86TVPCnvBYAxpgN4FV/MfjKWszoHuEZEDuALu1+IzxOajPcCY8xR699mfF9KVpGAv4/JKD7xlAvKVIJLFIWWLvoLK5PlLKDTcrWfAy4VkalWtsulVlvaYMXmHwTeM8Z8L+ilSXUvRKTc8ngQkTzgYnzrX5OunJUx5k5jTJXxFclci++9fZZJeC9EpEBEivyP8f1e7yQRfx/JzrRIxg++jI0P8MW8v5rs+UzQe3wEOAa48H0r+RK+OPVLwF7r31Krr+CrtL4P2AHUB9n5Ir6F1AbgC8l+X2O4D+fic/+3A9usnysm270AluErV7Xd+nD5htU+F98HZgPwGyDHas+1njdYr88NsvVV6/7sAS5P9ns7yftyPkPZbpPuXljv+V3rZ5f/8zARfx9aXkdRFEVJOJMx7KYoiqIkGRUfRVEUJeGo+CiKoigJR8VHURRFSTgqPoqiKErCUfFRFEVREo6Kj6IoipJw/j+gX7eOJq2lrAAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], "source": [ "plt.clf()\n", "# plt.plot(x_part, calcs, '.')\n", @@ -643,7 +688,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 12, "metadata": {}, "outputs": [], "source": [ @@ -659,7 +704,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 13, "metadata": {}, "outputs": [], "source": [ @@ -675,7 +720,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 14, "metadata": {}, "outputs": [], "source": [ @@ -688,7 +733,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 15, "metadata": {}, "outputs": [], "source": [ @@ -736,206 +781,13 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Tensorflow scaling" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# def scaling_func(x):\n", - "\n", - "# funcs = resonance(x, _mass = ztf.constant(jpsi_mass), scale = ztf.constant(jpsi_scale), phase = ztf.constant(jpsi_phase), width = ztf.constant(jpsi_width)) + resonance(x, _mass = ztf.constant(psi2s_mass), scale = ztf.constant(psi2s_scale), phase = ztf.constant(psi2s_phase), width = ztf.constant(psi2s_width))\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\n", - "\n", - "\n", - "# def s_func(x):\n", - " \n", - "# q = ztf.constant(x)\n", - " \n", - "# return zfit.run(scaling_func(q))\n", - " \n", - "\n", - "# print(integrate.quad(s_func, x_min, x_max, limit = 50))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# factor_jpsi = pdg[\"NR_auc\"]*pdg[\"jpsi_BR\"]/(pdg[\"NR_BR\"]*pdg[\"jpsi_auc\"])\n", - "# factor_jpsi = pdg[\"NR_auc\"]*pdg[\"jpsi_BR\"]/(pdg[\"NR_BR\"]*inte_fl)\n", - "# print(np.sqrt(factor_jpsi)*jpsi_scale)\n", - "# print(np.sqrt(factor_jpsi))\n", - "# # print(psi2s_scale)\n", - "# factor_psi2s = pdg[\"NR_auc\"]*pdg[\"psi2s_BR\"]/(pdg[\"NR_BR\"]*pdg[\"psi2s_auc\"])\n", - "# factor_psi2s = pdg[\"NR_auc\"]*pdg[\"psi2s_BR\"]/(pdg[\"NR_BR\"]*inte_fl)\n", - "# print(np.sqrt(factor_psi2s)*psi2s_scale)\n", - "# print(np.sqrt(factor_psi2s))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# def _t_f(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", - "# funcs = psi2s_res(xq) + jpsi_res(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 t_f(x):\n", - "# _ = np.array(x)\n", - "# probs = zfit.run(_t_f(_))\n", - "# return probs" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# print(36000*(1+ pdg[\"jpsi_BR\"]/pdg[\"NR_BR\"] + pdg[\"psi2s_BR\"]/pdg[\"NR_BR\"]))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# start = time.time()\n", - "\n", - "# result, err = integrate.quad(lambda x: t_f(x), x_min, x_max, limit = 5)\n", - "# print(result, \"{0:.2f} %\".format(err/result))\n", - "# print(\"Time:\", time.time()-start)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ "# Sampling\n", - "## One sample\n", - "! total_f.sample() always returns the same set !" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# nevents = int(pdg[\"number_of_decays\"])\n", - "# event_stack = 5000\n", - "\n", - "# calls = int(nevents/event_stack + 1)\n", - "\n", - "# total_samp = []\n", - "\n", - "# start = time.time()\n", - "\n", - "# samp = total_f.sample(n=event_stack)\n", - "# s = samp.unstack_x()\n", - "\n", - "# for call in range(calls):\n", - "\n", - "# sam = zfit.run(s)\n", - "# clear_output(wait=True)\n", - " \n", - "# # if call != 0:\n", - "# # print(np.sum(_last_sam-sam))\n", - " \n", - "# # _last_sam = sam\n", - " \n", - "# c = call + 1 \n", - "# print(\"{0}/{1}\".format(c, calls))\n", - "# print(\"Time taken: {}\".format(display_time(int(time.time() - start))))\n", - "# print(\"Projected time left: {}\".format(display_time(int((time.time() - start)/c*(calls-c)))))\n", - " \n", - "# with open(\"data/zfit_toys/toy_1/{}.pkl\".format(call), \"wb\") as f:\n", - "# pkl.dump(sam, f, pkl.HIGHEST_PROTOCOL)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# print(\"Time to generate full toy: {} s\".format(int(time.time()-start)))\n", - "\n", - "# total_samp = []\n", - "\n", - "# for call in range(calls):\n", - "# with open(r\"data/zfit_toys/toy_1/{}.pkl\".format(call), \"rb\") as input_file:\n", - "# sam = pkl.load(input_file)\n", - "# total_samp = np.append(total_samp, sam)\n", - "\n", - "# total_samp = total_samp.astype('float64')\n", - "\n", - "# data2 = zfit.data.Data.from_numpy(array=total_samp[:int(nevents)], obs=obs)\n", - "\n", - "# print(total_samp[:nevents].shape)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# bins = int((x_max-x_min)/7)\n", - "\n", - "# # calcs = zfit.run(total_test_tf(samp))\n", - "\n", - "# plt.hist(total_samp[:event_stack], bins = bins, range = (x_min,x_max))\n", - "\n", - "# # plt.plot(sam, calcs, '.')\n", - "# # plt.plot(test_q, calcs_test)\n", - "# plt.ylim(0, 20)\n", - "# # plt.xlim(3000, 3750)\n", - "\n", - "# plt.savefig('test2.png')\n", - "# 1-(0.21+0.62)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ "## Toys" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 24, "metadata": {}, "outputs": [], "source": [ @@ -1019,7 +871,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 25, "metadata": {}, "outputs": [], "source": [ @@ -1028,7 +880,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 26, "metadata": {}, "outputs": [], "source": [ @@ -1037,7 +889,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 27, "metadata": {}, "outputs": [], "source": [ @@ -1046,11 +898,21 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 28, "metadata": { "scrolled": false }, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "6/6 of Toy 1/1\n", + "Time taken: 1 min, 22 s\n", + "Projected time left: \n" + ] + } + ], "source": [ "# zfit.run.numeric_checks = False \n", "\n", @@ -1094,7 +956,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 29, "metadata": {}, "outputs": [], "source": [ @@ -1111,9 +973,18 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 30, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Time to generate full toy: 82 s\n", + "(5404696,)\n" + ] + } + ], "source": [ "print(\"Time to generate full toy: {} s\".format(int(time.time()-start)))\n", "\n", @@ -1135,9 +1006,29 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 31, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(5404696,)\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAD8CAYAAACVZ8iyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAXBklEQVR4nO3df7BcZZ3n8fd3IZoVMkOIFzZFGG6gApIZJLBXpGShVFB+WUF0Zgw1pdGhyMpg1TiO7kaxZtmyqGJ01V1qdrHCQgFVCDj8iHFgdkyxOjilIDcSQjBGEiYzXkkld8KAoRCXH9/9o8+FJvS9t2/36e7T3e9XVVef8/TTp5/z3NvPp8+PPh2ZiSRpuP2bXjdAktR7hoEkyTCQJBkGkiQMA0kShoEkiSbCICKOjojvRcS2iHg8Iv60KD88IjZGxBPF/cKiPCLi2ojYERFbIuLUTq+EJKk9zWwZvAT8eWaeCJwOXBERy4G1wP2ZuQy4v5gHOB9YVtzWANeV3mpJUqlmDYPM3J2ZPymm9wPbgKOAi4Cbi2o3Ax8spi8CbsmaB4HDImJx6S2XJJXm4LlUjohR4BTgIeDIzNwNtcCIiCOKakcBv6h72kRRtvuAZa2htuXAIYcc8u/f9ra3tdB8abg89stnZ61z0lG/3YWWTK9RG3vZpsd++WzP+6RTNm3a9C+ZOVLGspoOg4g4FLgL+HRm/ioipq3aoOwN17zIzHXAOoCxsbEcHx9vtinS0Bpde++sdcavubALLZleozb2sk2ja+/teZ90SkT8U1nLaupsooiYRy0Ibs3Mu4viPVO7f4r7vUX5BHB03dOXAE+V01xJUic0czZRADcA2zLza3UPbQBWF9OrgW/XlX+sOKvodODZqd1JkqRqamY30RnAR4HHImJzUfYF4BrgWxFxKfDPwB8Uj90HXADsAJ4HPlFqiyVJpZs1DDLzH2h8HADg7Ab1E7iizXbx4osvMjExwQsvvNDuoipt/vz5LFmyhHnz5vW6KZKG2JzOJuqmiYkJFixYwOjoKDMcrO5rmcm+ffuYmJhg6dKlvW6OpCFW2ctRvPDCCyxatGhggwAgIli0aNHAb/1Iqr7KhgEw0EEwZRjWUVL1VToMJEndUdljBgdq5ss2c7GrhS+hXHXVVRx66KF89rOfbfj4+vXrOf7441m+fHm7zZOkrnLLoETr16/npz/9aa+bIUlzZhjM4uqrr+aEE07gnHPOYfv27QBcf/31vOMd7+Dkk0/mwx/+MM8//zw//OEP2bBhA5/73OdYsWIFO3fubFhPkqrIMJjBpk2buP3223nkkUe4++67efjhhwH40Ic+xMMPP8yjjz7KiSeeyA033MC73vUuVq5cyVe+8hU2b97Mcccd17CeJFVR3xwz6IUf/OAHXHzxxbzlLW8BYOXKlQBs3bqVL37xizzzzDM899xznHvuuQ2f32w9Seo1w2AWjU79/PjHP8769es5+eSTuemmm/j+97/f8LnN1pOkXnM30QzOOuss7rnnHn7961+zf/9+vvOd7wCwf/9+Fi9ezIsvvsitt976av0FCxawf//+V+enqydJVdM3WwatnArarlNPPZWPfOQjrFixgmOOOYYzzzwTgC996Uu8853v5JhjjuGkk056NQBWrVrFZZddxrXXXsudd945bT1Jqpq+CYNeufLKK7nyyivfUH755Ze/oeyMM8543amll19+ecN6klQ17iaSJBkGkqSKh0HtpxEG2zCso6Tqq2wYzJ8/n3379g30YDn1ewbz58/vdVMkDbnKHkBesmQJExMTTE5O9ropHTX1S2eS1EuzhkFE3Ah8ANibmb9XlN0BnFBUOQx4JjNXRMQosA3YXjz2YGZ+spWGzZs3z1//kgbE6Np7e3J6uJrXzJbBTcBfAbdMFWTmR6amI+KrwLN19Xdm5oqyGihJ6rxZwyAzHyg+8b9B1K7V8IfAe8ttliSpm9o9gHwmsCczn6grWxoRj0TE30fEmW0uX5LUBe2GwSXAbXXzu4HfycxTgM8A34yI32r0xIhYExHjETE+6AeJJfVW2b+UOIhaDoOIOBj4EHDHVFlm/iYz9xXTm4CdwPGNnp+Z6zJzLDPHRkZGWm2GJKkE7WwZnAP8LDMnpgoiYiQiDiqmjwWWAU+210RJUqfNGgYRcRvwI+CEiJiIiEuLh1bx+l1EAGcBWyLiUeBO4JOZ+XSZDZYkla+Zs4kumab84w3K7gLuar9ZkqRuquzlKCRJ3WMYSJIMA6kfjK69t+nTIz2NUq0wDCRJhoEkyTCQJGEYSJIwDCRJGAaSJAwDSRKGgSQJw0CShGEgqQf8lnT1GAaSJMNAkmQYSJIwDCRJGAaSJAwDSRJNhEFE3BgReyNia13ZVRHxy4jYXNwuqHvs8xGxIyK2R8S5nWq4JKk8zWwZ3ASc16D865m5orjdBxARy4FVwO8Wz/lfEXFQWY2VJHXGrGGQmQ8ATze5vIuA2zPzN5n5j8AO4LQ22idJ6oJ2jhl8KiK2FLuRFhZlRwG/qKszUZS9QUSsiYjxiBifnJxsoxmSpHa1GgbXAccBK4DdwFeL8mhQNxstIDPXZeZYZo6NjIy02AxJUhlaCoPM3JOZL2fmK8D1vLYraAI4uq7qEuCp9pooSeq0lsIgIhbXzV4MTJ1ptAFYFRFvjoilwDLgx+01UZLUaQfPViEibgPeDbw1IiaA/wK8OyJWUNsFtAv4jwCZ+XhEfAv4KfAScEVmvtyZpkuqGq9G2r9mDYPMvKRB8Q0z1L8auLqdRklSGQyn5vkNZEl9bboB3yCYG8NAkmQYSJIMA0kShoEkCcNAkoRhIGlIeHbRzAwDSV3loFxNhoEkyTCQNDzcKpmeYSBJMgwkSYaBJAnDQJKEYSBJwjCQNGQ8o6gxw0CSZBhIVecnWXXDrGEQETdGxN6I2FpX9pWI+FlEbImIeyLisKJ8NCJ+HRGbi9s3Otl4SVI5mtkyuAk474CyjcDvZebbgZ8Dn697bGdmrihunyynmZKkTpo1DDLzAeDpA8q+m5kvFbMPAks60DZJUpeUcczgj4G/rZtfGhGPRMTfR8SZ0z0pItZExHhEjE9OTpbQDGnweLxgZvZPedoKg4i4EngJuLUo2g38TmaeAnwG+GZE/Faj52bmuswcy8yxkZGRdpohSWpTy2EQEauBDwB/lJkJkJm/ycx9xfQmYCdwfBkNlSR1TkthEBHnAf8ZWJmZz9eVj0TEQcX0scAy4MkyGipJ6pxmTi29DfgRcEJETETEpcBfAQuAjQecQnoWsCUiHgXuBD6ZmU83XLAk4X7/qjh4tgqZeUmD4humqXsXcFe7jZI0+EbX3suuay7sdTNU8BvIkoaOWyNvZBhIkgwDSZJhIFVWO7syhm03SP36Dtu6l8UwkCQZBpJ6z0/zvWcYSCqFA3p/MwwkSYaBJMkwkFQh7mrqHcNAkmQYSJIMA2lg9dsul35r76AxDKQB5gDbvmHpQ8NAkmQYSKq2uXwyH5ZP8Z1gGEgaOIbC3BkGFeU/swZNp/+nvcpre5oKg4i4MSL2RsTWurLDI2JjRDxR3C8syiMiro2IHRGxJSJO7VTjJbWuCgNgFdqgmma3DG4CzjugbC1wf2YuA+4v5gHOB5YVtzXAde03U5LKZRC9XlNhkJkPAE8fUHwRcHMxfTPwwbryW7LmQeCwiFhcRmMlVU/Zg+rU8hysu6udYwZHZuZugOL+iKL8KOAXdfUmirLXiYg1ETEeEeOTk5NtNEPSsDAgOqcTB5CjQVm+oSBzXWaOZebYyMhIB5oh9Yd+HOD6sc2NjK69t+l1GZR1nk47YbBnavdPcb+3KJ8Ajq6rtwR4qo3XkYbKXAaoVpffK2W+9qAPzt3WThhsAFYX06uBb9eVf6w4q+h04Nmp3UmSesdPwJpJs6eW3gb8CDghIiYi4lLgGuB9EfEE8L5iHuA+4ElgB3A98Celt1pS6VoJgV4Fh4FVvoObqZSZl0zz0NkN6iZwRTuNkiR1l99AlnrET7fTs2+6zzCQ1BcMiM4yDKSKqOJgZ5uGh2Eg9UC3B7R2X2/YB+BhWH/DQOqCqp7bX/VBrurtGySGgaSmOTgPbh8YBlIFDOoA02+avUjeIP69DANJorkBfhBDYIphIGlG0w2AgzwwDiPDQFLlGTydZxhImpNhGZiHZT2nGAaS1KJBCgzDQJJkGEh6zSB90m3VsPaBYSCpbcM6gA4Sw0CSg7kMA2kYOfjrQIaBJKm5n71sJCJOAO6oKzoW+AvgMOAyYLIo/0Jm3tdyCyV1VLNbCW5NDLaWwyAztwMrACLiIOCXwD3AJ4CvZ+Z/K6WFkqSOazkMDnA2sDMz/ykiSlqkNHyq8um7Ku1Q95R1zGAVcFvd/KciYktE3BgRCxs9ISLWRMR4RIxPTk42qiL1naoPolVvn3qn7TCIiDcBK4G/LoquA46jtgtpN/DVRs/LzHWZOZaZYyMjI+02Q6qcAwdeB+LB0uxvH/SLMrYMzgd+kpl7ADJzT2a+nJmvANcDp5XwGpKkDiojDC6hbhdRRCyue+xiYGsJryFJ6qC2wiAi3gK8D7i7rvjLEfFYRGwB3gP8WTuvIZWhm5vyc9k9NCi7GNT/2gqDzHw+Mxdl5rN1ZR/NzJMy8+2ZuTIzd7ffTKm6HNB1oH78n/AbyFKX9eNAoZkNwt/UMJAkGQaSJMNA6qhB2H2g4WAYaCDNdRAeXXtvS8/RYBumv7FhIDVpmAYGDR/DQJojQ0GDyDCQJBkGklSGft9iNAykktQPBv0+MGj4GAaSJMNAg6edT+V+otewMgykBmb74RJ/uEaDpqzfQJZ6ruwBuZnl+UU1DQq3DKQmOIhr0BkGGnoO9JJhoCFmCKgT+vX/yjDQQGj1Ddivb1xVUz//P7V9ADkidgH7gZeBlzJzLCIOB+4ARoFdwB9m5r+2+1oaXqNr72XXNReWtixJr1fWlsF7MnNFZo4V82uB+zNzGXB/MS+9joOyVB2d2k10EXBzMX0z8MEOvY4kqQRlhEEC342ITRGxpig7MjN3AxT3Rxz4pIhYExHjETE+OTlZQjM06NrdknBLRJpeGV86OyMzn4qII4CNEfGzZp6UmeuAdQBjY2NZQjskSS1qe8sgM58q7vcC9wCnAXsiYjFAcb+33deRWuUWhTS7tsIgIg6JiAVT08D7ga3ABmB1UW018O12Xkf9r5XfGJ5pWfX3ktrX7m6iI4F7ImJqWd/MzP8TEQ8D34qIS4F/Bv6gzdeRKsEA0qBqa8sgM5/MzJOL2+9m5tVF+b7MPDszlxX3T5fTXKl5DtzqtX76H/QbyOor/fTmkvqJYaCWdWJgbnQ8YLbXmelxf3dAvdJv/2uGgbpqLgO3pO4xDNRRDvBSfzAMVIpBGvQHaV2kZhkGFeRgJKnbDANVggEo9ZZhoJaU/W3iuT5ueEjlMgwkSYaByjOXawb5yV6qFsNATSnzQnONlt2L50rd0C//o4aBGpppX73f6pUGj2EwZDo9cJexi8hLVGvQ9MP/smEwhOZy3R9Jw8Ew6BAHWftA6ieGgWY9l99BXRp8hoHewMFfGj6GgSSp9d9AjoijgVuAfwe8AqzLzP8REVcBlwGTRdUvZOZ97Ta0n/TDJ+tutrEf+kMadu1sGbwE/HlmngicDlwREcuLx76emSuK28AFQT8MbvXfB+iH9krqrZa3DDJzN7C7mN4fEduAo8pqmCSpe0o5ZhARo8ApwENF0aciYktE3BgRC8t4jWHRyqf4mbYAerVV4NaI1F/aDoOIOBS4C/h0Zv4KuA44DlhBbcvhq9M8b01EjEfE+OTkZKMqqoB2foxeUv9oKwwiYh61ILg1M+8GyMw9mflyZr4CXA+c1ui5mbkuM8cyc2xkZKSdZlRKK4Nn2QNqq1sXnbwYnaRqazkMIiKAG4Btmfm1uvLFddUuBra23rz+5sAqqV+0s2VwBvBR4L0Rsbm4XQB8OSIei4gtwHuAPyujoVVRlQG+Ku2QNBjaOZvoH4Bo8NBAnkq665oLm6rXyrIlqdeG9hvIVRiEywiPZr9LUIX1lVRdQxsG3eSBWUlVZxj0wGxXCZ2tXivLlqSZGAYzmGlgbXbQ7cTg7IAvqWyGQZ1Wd+d0MhjKvL6QISL1VpXfg4ZBRbhLSFIvGQaSpOEIg3Y+QfvpW9IwGIowmE67A71BIWlQDHUYzKTTF5QzSKThVNXvHQ1VGDR7fn8rdTr5x/UXy6T+V/X378CGwUwDf/1jjQbaufzRqv4HlqRmDGwYgAO1JDVroMOgnsEgSdMbqDCo0m8AS9JMqjY2DVQYSJJa0/KP21RV1dJWkvqBWwaSJMNAktTBMIiI8yJie0TsiIi1nXodcNeQJLWrI2EQEQcB/xM4H1gOXBIRyzvxWpKk9nVqy+A0YEdmPpmZ/w+4HbioQ68FuHUgSe3o1NlERwG/qJufAN5ZXyEi1gBritnnImIf8C8dak+/eSv2xRT7osZ+eM3A9EX8ZVtPfytwTDkt6VwYRIOyfN1M5jpg3atPiBjPzLEOtaev2BevsS9q7IfX2Bc1RT+MlrW8Tu0mmgCOrptfAjzVodeSJLWpU2HwMLAsIpZGxJuAVcCGDr2WJKlNHdlNlJkvRcSngL8DDgJuzMzHZ3naulkeHyb2xWvsixr74TX2RU2p/RCZOXstSdJA8xvIkiTDQJJUkTDo5qUreiUiboyIvRGxta7s8IjYGBFPFPcLi/KIiGuL/tgSEafWPWd1Uf+JiFjdi3VpR0QcHRHfi4htEfF4RPxpUT5UfRER8yPixxHxaNEP/7UoXxoRDxXrdEdxAgYR8eZifkfx+Gjdsj5flG+PiHN7s0bti4iDIuKRiPibYn4o+yIidkXEYxGxOSLGi7LOvz8ys6c3ageYdwLHAm8CHgWW97pdHVjPs4BTga11ZV8G1hbTa4G/LKYvAP6W2vc1TgceKsoPB54s7hcW0wt7vW5z7IfFwKnF9ALg59QuWTJUfVGsz6HF9DzgoWL9vgWsKsq/AVxeTP8J8I1iehVwRzG9vHjPvBlYWryXDur1+rXYJ58Bvgn8TTE/lH0B7ALeekBZx98fVdgy6PqlK3ohMx8Anj6g+CLg5mL6ZuCDdeW3ZM2DwGERsRg4F9iYmU9n5r8CG4HzOt/68mTm7sz8STG9H9hG7RvrQ9UXxfo8V8zOK24JvBe4syg/sB+m+udO4OyIiKL89sz8TWb+I7CD2nuqr0TEEuBC4H8X88GQ9sU0Ov7+qEIYNLp0xVE9aku3HZmZu6E2SAJHFOXT9clA9VWxeX8KtU/FQ9cXxW6RzcBeam/WncAzmflSUaV+nV5d3+LxZ4FFDEA/FP478J+AV4r5RQxvXyTw3YjYFLXL9kAX3h9V+KWzWS9dMYSm65OB6auIOBS4C/h0Zv6q9sGucdUGZQPRF5n5MrAiIg4D7gFObFStuB/YfoiIDwB7M3NTRLx7qrhB1YHvi8IZmflURBwBbIyIn81Qt7S+qMKWwTBfumJPsUlHcb+3KJ+uTwairyJiHrUguDUz7y6Kh7IvADLzGeD71Pb5HhYRUx/S6tfp1fUtHv9tarsdB6EfzgBWRsQuaruJ30ttS2EY+4LMfKq430vtQ8JpdOH9UYUwGOZLV2wApo7yrwa+XVf+seJMgdOBZ4tNw78D3h8RC4uzCd5flPWNYt/uDcC2zPxa3UND1RcRMVJsERAR/xY4h9rxk+8Bv19UO7Afpvrn94H/m7UjhRuAVcUZNkuBZcCPu7MW5cjMz2fmkqxddG0VtXX7I4awLyLikIhYMDVN7f96K914f/T6yHndEfGfU9tnemWv29OhdbwN2A28SC21L6W2n/N+4Ini/vCiblD7caCdwGPAWN1y/pjagbEdwCd6vV4t9MN/oLa5ugXYXNwuGLa+AN4OPFL0w1bgL4ryY6kNYDuAvwbeXJTPL+Z3FI8fW7esK4v+2Q6c3+t1a7Nf3s1rZxMNXV8U6/xocXt8ajzsxvvDy1FIkiqxm0iS1GOGgSTJMJAkGQaSJAwDSRKGgSQJw0CSBPx/Vnu3Jq1jkDcAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], "source": [ "plt.clf()\n", "\n", @@ -1162,7 +1053,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 32, "metadata": {}, "outputs": [], "source": [ @@ -1185,7 +1076,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 33, "metadata": {}, "outputs": [], "source": [ @@ -1194,7 +1085,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 34, "metadata": {}, "outputs": [], "source": [ @@ -1210,17 +1101,73 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 35, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.6628998341322463\n", + "6.258685216704922\n", + "2.249780125960397\n", + "-4.790281653148225\n", + "-0.546187997552245\n", + "-4.053408742543857\n", + "3.3921906479017068\n", + "1.380997849250071\n", + "0.5628459148842957\n", + "------------------------------------------------------------------\n", + "| FCN = -7.175E+05 | Ncalls=253 (253 total) |\n", + "| EDM = 2.21E-05 (Goal: 5E-06) | up = 0.5 |\n", + "------------------------------------------------------------------\n", + "| Valid Min. | Valid Param. | Above EDM | Reached call limit |\n", + "------------------------------------------------------------------\n", + "| True | True | False | False |\n", + "------------------------------------------------------------------\n", + "| Hesse failed | Has cov. | Accurate | Pos. def. | Forced |\n", + "------------------------------------------------------------------\n", + "| False | True | True | True | False |\n", + "------------------------------------------------------------------\n", + "Function minimum: -717521.3025826928\n", + "---------------------------------------------------------------------------------------------\n", + "| | Name | Value | Hesse Err | Minos Err- | Minos Err+ | Limit- | Limit+ | Fixed |\n", + "---------------------------------------------------------------------------------------------\n", + "| 0 | jpsi_p | 1.486 | 0.016 | | |-6.28319 | 6.28319 | |\n", + "| 1 | p3770_p | 2.55 | 0.09 | | |-6.28319 | 6.28319 | |\n", + "| 2 | psi2s_p | 1.344 | 0.027 | | |-6.28319 | 6.28319 | |\n", + "| 3 | omega_p | -5.80 | 0.22 | | |-6.28319 | 6.28319 | |\n", + "| 4 | p4040_p | -3.37 | 0.16 | | |-6.28319 | 6.28319 | |\n", + "| 5 | rho_p | -0.19 | 0.30 | | |-6.28319 | 6.28319 | |\n", + "| 6 | phi_p | -5.59 | 0.15 | | |-6.28319 | 6.28319 | |\n", + "| 7 | p4160_p | 4.04 | 0.08 | | |-6.28319 | 6.28319 | |\n", + "| 8 | p4415_p | -3.13 | 0.13 | | |-6.28319 | 6.28319 | |\n", + "---------------------------------------------------------------------------------------------\n", + "-------------------------------------------------------------------------------------\n", + "| | jpsi_p p3770_p psi2s_p omega_p p4040_p rho_p phi_p p4160_p p4415_p |\n", + "-------------------------------------------------------------------------------------\n", + "| jpsi_p | 1.000 -0.133 0.219 0.015 -0.176 -0.110 0.023 -0.065 -0.143 |\n", + "| p3770_p | -0.133 1.000 -0.527 -0.007 0.014 0.005 -0.011 0.001 0.014 |\n", + "| psi2s_p | 0.219 -0.527 1.000 0.016 -0.312 -0.060 0.025 -0.006 -0.140 |\n", + "| omega_p | 0.015 -0.007 0.016 1.000 -0.007 -0.255 0.034 -0.011 -0.012 |\n", + "| p4040_p | -0.176 0.014 -0.312 -0.007 1.000 0.001 -0.011 -0.381 -0.104 |\n", + "| rho_p | -0.110 0.005 -0.060 -0.255 0.001 1.000 -0.246 0.026 0.017 |\n", + "| phi_p | 0.023 -0.011 0.025 0.034 -0.011 -0.246 1.000 -0.018 -0.018 |\n", + "| p4160_p | -0.065 0.001 -0.006 -0.011 -0.381 0.026 -0.018 1.000 0.164 |\n", + "| p4415_p | -0.143 0.014 -0.140 -0.012 -0.104 0.017 -0.018 0.164 1.000 |\n", + "-------------------------------------------------------------------------------------\n", + "Hesse errors: OrderedDict([(, {'error': 0.015523017204917444}), (, {'error': 0.0949040578205631}), (, {'error': 0.027368429790692872}), (, {'error': 0.21984151021717713}), (, {'error': 0.1583791551448439}), (, {'error': 0.2962129131334841}), (, {'error': 0.15237884170408078}), (, {'error': 0.08086187469444539}), (, {'error': 0.12776706213127875})])\n" + ] + } + ], "source": [ "start = time.time()\n", "\n", "for param in total_f.get_dependents():\n", " param.randomize()\n", " \n", - "for param in total_f.get_dependents():\n", - " print(zfit.run(param))\n", + "# for param in total_f.get_dependents():\n", + "# print(zfit.run(param))\n", " \n", "nll = zfit.loss.UnbinnedNLL(model=total_f, data=data2, fit_range = (x_min, x_max))\n", "\n", @@ -1240,9 +1187,17 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 36, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Time taken for fitting: 2 min, 6 s\n" + ] + } + ], "source": [ "print(\"Time taken for fitting: {}\".format(display_time(int(time.time()-start))))\n", "\n", @@ -1254,16 +1209,29 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 60, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ8AAAD8CAYAAACo9anUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO29eXgc1ZX3/znd2mxJlmTZsmXLtrzINgYbYwtjtsQBAgYTTBJITMKESUh4kwm/rPML8GSbZMLMMG8Sz5tJMoEJeYdsGEIWHGLWACEQbLwvMl7kBVveZO2b1VJ33/ePqm61Wl3dLVlSV0vn8zx+VH3r1r23ylJ9+5x77rlijEFRFEVRhhNPqgegKIqijD5UfBRFUZRhR8VHURRFGXZUfBRFUZRhR8VHURRFGXZUfBRFUZRhJynxEZGVIrJfRKpF5P4Y57NF5An7/CYRKY8494Bdvl9EbkjUpojMtNs4aLeZZZe/S0S2iYhfRG6L6v8uu/5BEbmr/49BURRFGU4Sio+IeIEfATcCC4A7RGRBVLW7gUZjzBxgLfCQfe0CYA1wIbAS+LGIeBO0+RCw1hhTATTabQMcA/4e+HXU+MYD3wQuA5YB3xSRomQfgKIoijL8JGP5LAOqjTGHjTFdwDpgdVSd1cBj9vFTwLUiInb5OmOMzxhzBKi224vZpn3NNXYb2G3eCmCMOWqM2QUEo/q+AXjRGNNgjGkEXsQSOkVRFMWlZCRRZypwPOJzDZaVEbOOMcYvIs1AsV2+MeraqfZxrDaLgSZjjD9G/f6Mr881InIPcA9Abm7u0vnz5ydoVlGUoDFUnWxh8rgcJuZnh8vPtvo43dLJhVMK8EgKBwi0dnZztL6DORPz6AoEOdbQQUVJPjmZOqU92GzdurXOGDNxMNpKRnxi/WpF5+RxquNUHuu3Il79eCR1jTHmEeARgMrKSrNly5YEzSqK0tLZzaJ/eoGvrbqAT149K1z+yGuH+JcN+/jrt24gNzuZ18jQ8eLeM3zq51tYd+9VHG/s4B9+tY2nvnA18yePS+m4RiIi8s5gtZXMV4MaYFrE5zLgpFMdEckACoCGONc6ldcBhXYbTn0NZHyKogwAYzu5PdL7O17oc9AFuSFD+SlFCFthLhiWkoBkxGczUGFHoWVhBRCsj6qzHghFmd0GvGys34j1wBo7Gm4mUAG85dSmfc0rdhvYbT6dYHzPA9eLSJEdaHC9XaYoynkSsN/i0a41CYvPcI+oL72H4B5RVOKTUHzs+Zd7sV7obwNPGmOqROTbInKLXe1RoFhEqoEvAffb11YBTwJ7geeAzxpjAk5t2m3dB3zJbqvYbhsRuVREaoDbgYdFpMruowH4ZyxB2wx82y5TFOU8Cdjq4vVEWz7WTzdkxQ8NQS2f9CIpZ60xZgOwIarsGxHHnViiEOvaB4EHk2nTLj+MFQ0XXb4Zy6UWq4+fAT+LexOKovSbkLh4+oiP9TngBtOHkHUmYYssFeLT3d1NTU0NnZ2dw9/5IJOTk0NZWRmZmZlD1kdqZwoVRXE1PW63KPHxuMftFoywfEKjNAnjlAafmpoa8vPzKS8vD4tgOmKMob6+npqaGmbOnDlk/WgsoqIojoRe7N4+AQfWT1e53RA8nt5lw0lnZyfFxcVpLTxgzecVFxcPuQWn4qMoiiPBYE8kWSQeVwUc9IxRUhxwkO7CE2I47kPFR1EURxIFHLghqqzH8ukRydSPSkmEio+iKI6ExCVafMRN63zsn9Ir4CD143Irr776KjfffDMAPp+P6667jsWLF/PEE08M6zg04EBRFEeC4QWcDotMozMtpoDIRabhgAPVnqTYvn073d3d7NixY9j7VstHURRHnAIOvJ7Q+dS/5SPdbiFRTP2oUsPRo0eZP38+d911F4sWLeK2226jo6OD5557jvnz53PVVVfxu9/9DoDa2lruvPNOduzYweLFizl06NCwjlUtH0VRHAnN+URnOHBVeh16rLOQRgZTHAnxrT9Wsfdky6C2uWDKOL75vgsT1tu/fz+PPvooV155JZ/4xCf4/ve/z8MPP8zLL7/MnDlz+PCHPwxASUkJP/3pT/nud7/LM888M6hjTQa1fBRFcSQsPo5zPsM+pD5owEFvpk2bxpVXXgnAnXfeyZYtW5g5cyYVFRWICHfeeWeKR2ihlo+iKI6YdFrn44JQ6xDJWChDRfT8XHNzsytDwNXyURTFkXCGg6g3RTi9jgvEJxiRhSH8jk39sFLGsWPHePPNNwF4/PHHue666zhy5Eh4Tufxxx9P5fDCqPgoiuJI0Cm9jpui3eyfVmLR0R1wAHDBBRfw2GOPsWjRIhoaGvjiF7/II488wqpVq7jqqquYMWNGqocIqNtNUZQ4BINO4mOfd4HlY2JYPm4YV6rweDz85Cc/6VW2cuVK9u3b16fuihUrWLFixTCNrDdq+SiK4ohzhoPUZY+OJhT04BHRLRXSCBUfRVEciXyxR+Jx0TqfYMQi09G+mVx5eTl79uxJ9TCSQsVHURRHeuZ8epe7Kr1OZLRbikOt3RD9NxgMx32o+CiK4ohTbjdXLTKNmPPxpFB9cnJyqK+vT3sBCu3nk5OTM6T9aMCBoiiOBMJbKkSl13HRItNI16CEy4Z/YGVlZdTU1HD27Nlh73uwCe1kOpSo+CiK4oiz5WOfd4H6hOd8SG0gRGZm5pDu/DnSULeboiiOhNbxRGc4cGN6HQ21Ti9UfBRFcSTQK5KsBzel1wlbPhFvs9SPSkmEio+iKI4YJ7ebxz3pdSItHzetP1Lio+KjKIojAdvt5phexwUv+cg5H3GRRabER8VHURRHAmHLp3e5q9Lr2D97WT6pG46SJCo+iqI4ErmGJpIe91bqX/ORGQ404CB9UPFRFMWRgGNiURdltY6xzke1x/2o+CiK4khoTic64MBNFkYwGGn5qNstXVDxURTFkcgXeySuSq9j/4xc5+MGd6ASHxUfRVEcCTiEWoc+uynazSOpzXCg9A8VH0VRHHHez8f66QbLJySAkuLcbkr/UPFRFMWRkPhkeHq/KtyVXseE3W09brfUjUdJDhUfRVEc8SewfNwwt2JMj7tN1/mkDyo+iqI4ErBjqTMc9vMJuMD0CRrTZ7M7dbu5HxUfRVEccbZ83ON2C5oeN6AnbJKlcEBKUqj4KIriSCAQmvOJTixq/XSDhWGMCQcaaMBB+pCU+IjIShHZLyLVInJ/jPPZIvKEfX6TiJRHnHvALt8vIjckalNEZtptHLTbzIrXh4hkishjIrJbRN4WkQcG+jAURelNIsvHDXM+ltvNGk8Kd9FW+klC8RERL/Aj4EZgAXCHiCyIqnY30GiMmQOsBR6yr10ArAEuBFYCPxYRb4I2HwLWGmMqgEa7bcc+gNuBbGPMQmAp8L8ixU9RlIETCFrzKdHbaLvJ7WYFHFjHus4nfUjG8lkGVBtjDhtjuoB1wOqoOquBx+zjp4BrxfptXQ2sM8b4jDFHgGq7vZht2tdcY7eB3eatCfowQK6IZABjgC6gJeknoCiKIwFj+oRZg/vW+YQtn3BZ6selxCcZ8ZkKHI/4XGOXxaxjjPEDzUBxnGudyouBJruN6L6c+ngKaAdOAceA7xpjGqJvQkTuEZEtIrLl7NmzSdy2oiiBoOnjcoOIdT4uMH2Cvdb5aKh1upCM+PT9zev7f+tUZ7DK4/WxDAgAU4CZwJdFZFafisY8YoypNMZUTpw4MUZTiqJE4w+YPsEG0DMH5IZQa2uRadScj1o+ricZ8akBpkV8LgNOOtWx3V8FQEOca53K64BCu43ovpz6+AjwnDGm2xhTC7wBVCZxX4qiJCAQDOL1xhCf0DofF7zjDT1uQN1SIX1IRnw2AxV2FFoWVgDB+qg664G77OPbgJeN9dVjPbDGjlSbCVQAbzm1aV/zit0GdptPJ+jjGHCNWOQCy4F9yT8CRVGc8AcdLB+vu9xufTIcqPq4noxEFYwxfhG5F3ge8AI/M8ZUici3gS3GmPXAo8AvRKQayxpZY19bJSJPAnsBP/BZY0wAIFabdpf3AetE5DvAdrttnPrAipr7v8AerC8+/9cYs2vAT0RRlDBOcz4hQfK7Qnzo43ZzwbCUBCQUHwBjzAZgQ1TZNyKOO7FCnmNd+yDwYDJt2uWHseZxostj9mGMaXPqW1GU88OyfGJFu4XmfFK/lanRgIO0RDMcKIriSCLLJ5B67em1zkcDDtIHFR9FURxxmvPxeNxj+fTKcGCXqfa4HxUfRVEcCQSDMS0fsKwft8z59N1SIfXjUuKj4qMoiiP+QGy3G1hrfQIuMDEisxlowEH6oOKjKIojTnM+YIuPSxb6hGIiNLdb+qDioyiKI1ZuN/dbPp6oxKea2839qPgoiuJIPMsnwyOuSK/TK7Fo7KEqLkTFR1EUR6zcbrFfE17XBBz0bCbncVHCUyU+Kj6KojiSaM7HDS95Q4/F46Z9hpT4qPgoiuKIPxgkI0ZiUYAMj8cVlo/pldvNKnPDXJQSHxUfRVEciWf5eDzu2FIhGIyc8xE8ohkO0gEVH0VRHHHKcACW5eMK8YnI7QaWELlhXEp8VHwURXEk4TofF7zkrTmfnjF6XBICrsRHxUdRFEecslqDtaGc3wW53UxEtBtY43JDIIQSHxUfRVEcCQRNOIloNJblM8wDikG0deb1iEa7pQEqPoqiOOIPBp3nfLziiqzWAUMvgRRxRyCEEh8VH0VRHAnESSzqEZcsMg0aIqPBLcsn9eNS4qPioyiKI91BQ6Y39msiwyUv+T5uN3HHuJT4qPgoiuJIdyBIlsMiU49H8Lsgq3UgKrGoiDvmopT4qPgoiuJItz8Y1/Jxw9xKsE/AgeZ2SwdUfBRFcaQ7YMhwEB+3bKkQMH3dbm4YlxIfFR9FUWJijKErjtvNLYtMg8HebjePS+ailPio+CiKEpOQsMRzu7llzifS8vHoItO0QMVHUZSYdNvCkpnh7HZzg4URiEgsCiF3YAoHpCSFio+iKDHpskPGnCwf12wmFzREDtEjGnCQDqj4KIoSk+6w+DjN+Xhc8ZLvE3DgEotMiY+Kj6IoMelOYPlkuMjy6RVwoFsqpAUqPoqixCQUTOAkPm55yccMOFDLx/Wo+CiKEpOuBG43tywy9QcMXtGs1umGio+iKDEJud2ynAIOvC5xu/WxfDSrdTqg4qMoSky6/dYL3DHDgbhkS4Wo9Dq6yDQ9UPFRFCUmidxurgm1Nr03vPO6ZC5KiY+Kj6IoMUnkdsv0uuMlHwj2nvNRyyc9UPFRFCUm/gQZDjK9nrBApZI+bjcBF3gDlQQkJT4islJE9otItYjcH+N8tog8YZ/fJCLlEecesMv3i8gNidoUkZl2GwftNrOS6GORiLwpIlUisltEcgbyMBRF6SHROh9LfAwmxVZG0MRKr6OWj9tJKD4i4gV+BNwILADuEJEFUdXuBhqNMXOAtcBD9rULgDXAhcBK4Mci4k3Q5kPAWmNMBdBotx2vjwzgl8CnjTEXAiuA7n4+B0VRokg055NlW0TdKU6kFuiTXscd7kAlPslYPsuAamPMYWNMF7AOWB1VZzXwmH38FHCtiIhdvs4Y4zPGHAGq7fZitmlfc43dBnabtybo43pglzFmJ4Axpt4YE0j+ESiKEotkMhxE1ksVgeiAA4+k3BpTEpOM+EwFjkd8rrHLYtYxxviBZqA4zrVO5cVAk91GdF9OfcwFjIg8LyLbROQrsW5CRO4RkS0isuXs2bNJ3LaijG6ScbtF1ksVweiAA91MLi1IRnxi2dzR/7NOdQarPF4fGcBVwEftn+8XkWv7VDTmEWNMpTGmcuLEiTGaUhQlkvCWCg5ut0y3uN1ipNdxQRyEkoBkxKcGmBbxuQw46VTHnoMpABriXOtUXgcU2m1E9xWvj78YY+qMMR3ABmBJEvelKEocElk+oR1OU2n5GGMwfQIOdEuFdCAZ8dkMVNhRaFlYAQTro+qsB+6yj28DXjaW03U9sMaOVJsJVABvObVpX/OK3QZ2m08n6ON5YJGIjLVF6d3A3uQfgeLzB6g62ZzqYSguw9dtiUp2nFBrSK34hAILNLFo+pFQfOz5lXuxXvJvA08aY6pE5Nsicotd7VGgWESqgS8B99vXVgFPYonBc8BnjTEBpzbttu4DvmS3VWy3Ha+PRuD7WIK2A9hmjPnTQB/IaORrv9/Dqh+8zpmWzlQPRXERnX4rbicn0xvzfIYbxMfEEB8NtU4LMhJXAWPMBix3VmTZNyKOO4HbHa59EHgwmTbt8sNY0XDR5fH6+CVWuLUyALYdawSgtbObSeN0iZRiEbJ8nDIchNxuXf7UvehDi0l7ud1E1O2WBmiGAyWMfllUIvH5g2RleHqFMUfiCrdb2PLpKdMtFdIDFR9FUWLS2R1wnO+BHvHxpzCXTSAQEp+ecYpuqZAWqPgoYST2F1xllOLzB8nOiD3fAz3ik0q3W3ewbxYGrwYcpAUqPoqixMTXHSAn0/kVkZWR+lDrUPLTjAjLx+uSHVaV+Kj4KIoSE8vycX5FhF74qRSfUN8Z3qjEoio+rkfFR1GUmPj8Accwa3BHwEFoM7tIt5tbtnpQ4qPioyhKTDq741s+IbdbVwrT6/hDlk+E2y3DJTusKvFR8VEUJSY+fyCpgAN/St1ufS0fr1fCc0GKe1HxURQlJp3dwbgBB+5wu/W1fDI9npSGfyvJoeKjhNHoVCWSZC2fVLrdQpZPZMBBhtdaZKpZDtyNio+iKDHx+RNZPnaotT+VodZ9M2+HLTK1flyNio8SRheZKpFYGQ7SI9otIyIFUOhY533cjYqPoigx8fmDZLt8zqcrZPlEROVlhAMhVHzcjIqPoigx6exOtM4nlNU69RkOMiMDDkLuQHW7uRoVH0VR+hAMGivaLc46HxEhO8NDpwvmfHoFHHjU8kkHVHwURenDuW5rI7nc7Phbfo3J8tJp100F3TEyHGS4YHtvJTEqPoqi9KG9yw/A2ATik5ORWvGJleEgJESa5cDdqPgoYXSdjxLiXJdt+WQ5z/kA5GR66OxO/ZxPbLebWj5uRsVHUZQ+tPss8RmblcDyyUy12y3WOp+Q202/TbkZFR9FUfrQYbvdcrPjWz7Zmd4UBxzEWueT+h1WlcSo+ChhdJGpEqK9K0nLJ8OTWssnHO0Wuc5HLZ90QMVHUZQ+dPiSs3xyMr34Uio+lsBkxUivo3M+7kbFR0G/HyrRtIcDDhLN+aQ24MDnt8YZue9QOL2ORru5GhUfRVH6EJrzGZsw2s0bXhOUCnz+IFleD57IOZ8Upv3pDgT50MNv8p9/PjjsfacbKj4KOtWjRJN0tFuK1/n4Yuy2Gl7nk4I5nwNnWnnrSAPfe/HAsPedbqj4KGF0nY8SoqPLjwhxt1SAkNstlZZPoE/y01RGu9W2+MLHje1dw95/OqHio+icj9KHjq4AuVkZSIIQyJwUh1r7/ME+2z6kcp1PbWtn+PhIffuw959OqPgoYVSElBCtnd0JI93AWufT5Q+mbNdQS3yiLB9v6iyfMxGWz5nmzjg1FRUfJTzno243JUTLOT8FYzIT1gu55Xwpsn46uwNkZzpYPn7rF/q5Pae599fbMMPwC17b2hleL3dSxScuKj5KGKO2j2LTfK47KfEZY7/4UzXvE8vyCe1B1GmHYX/197t5ZtcpGoZhDuZMi4+5JfnkZHo43XxuyPtLZ1R8FJUcpQ8tnd2My0ksPqFQ7I5UiU93oI/4hD777PVH9bbonG4ZekukttVHybhsSgvGqOWTABUfJYy63ZQQyVo+edlWnbZO/1APKSbWVt+93W4hyye0ADVEXdvQWz61LZ2U5OcwaVy2zvkkQMVH0TkfpQ8t57oZl4z45FjrgNp83UM9pJjEDDjwCB4hnHkhZJ3Vtfr6XD+YBIOGs60+Jo3LZvK4nGGxtNKZ+CvIlFGFzvkoYL1EW31+xuUkfj3k2RFxbb5Uzfn0dbtZ23t7w5ZPfk4GHV0B6tuHVnwaOrrwBw0l+dkEjKG2xYcxJmG4+mhFxUdRyVF60erzYwzJWT6pdrt1913nA1YUXigCL7TotH6I3W6hBaaTxuVggK5AkIb2Lorzsoe033QlKbebiKwUkf0iUi0i98c4ny0iT9jnN4lIecS5B+zy/SJyQ6I2RWSm3cZBu82sRH3Y56eLSJuI/GN/H4JioW43BSyXGyQpPil2u3V0+WPmn8uOSPsTsoCGes4ntMC0xHa7wfAEOaQrCcVHRLzAj4AbgQXAHSKyIKra3UCjMWYOsBZ4yL52AbAGuBBYCfxYRLwJ2nwIWGuMqQAa7bYd+4hgLfBssjeuKEpsmm3xSSrgICskPqlxu7V3BcjN7uvAyY6wfDrDUW9D63YLWT4l+TlMKrDE54yKjyPJWD7LgGpjzGFjTBewDlgdVWc18Jh9/BRwrViOztXAOmOMzxhzBKi224vZpn3NNXYb2G3emqAPRORW4DBQlfytKyE04ECJpKnDEp/CJMQnlAUhFW637kCQLn8wPO8USU5My2eIxce2fCbmR1g+zUPbZzqTjPhMBY5HfK6xy2LWMcb4gWagOM61TuXFQJPdRnRfMfsQkVzgPuBb8W5CRO4RkS0isuXs2bMJbnl0oZqjRBJ6SU/ITzxXkeH1MCbTS3vX8ItPuy+07YOz5RMImnCOt6Ge8znT4qNgTCY5mV4m5mcjom63eCQjPrFCNaLfV051Bqs8Xh/fwnLTtcU431PRmEeMMZXGmMqJEyfGqzpq0Wg3BSLEJze5ifLc7AxaU2D5hDa8y4vhdsvJ8OLrDoatHhFLfIYyxU5tayeTxlnPLNPrYUKervWJRzLRbjXAtIjPZcBJhzo1IpIBFAANCa6NVV4HFIpIhm3dRNZ36uMy4DYR+XegEAiKSKcx5odJ3JsSgbrdFLAm5jO9wrgxyQXD5udk0OZLoeUTw+2WnemhzecPZzmYUjCGE03n7BDyxO7EgXCmxUdJfk74s671iU8yls9moMKOQsvCCiBYH1VnPXCXfXwb8LKxvmKsB9bYkWozgQrgLac27WtesdvAbvPpeH0YY642xpQbY8qB/wD+RYVHUQZOfZuP4tzspNen5GVn0NY5/NFuIfGJGXCQ4cHXHQznd5taNAYYWtfbmZZOJo3rEZ9J43I04CAOCcXHtkDuBZ4H3gaeNMZUici3ReQWu9qjWPMv1cCXgPvta6uAJ4G9wHPAZ40xAac27bbuA75kt1Vst+3YhzJ4qOGjgJULbUJ+VtL1C8dm0nQuFeJjCUtuzDkfL53+QNjymVpoic9QBR0EgobaVh+TC3pclZMLstXyiUNSdrUxZgOwIarsGxHHncDtDtc+CDyYTJt2+WGsaLjocsc+Iur8U7zzSnyGI+W84n7qbMsnWYrGZnG8oWMIRxSbtrDlE2udj8ee8+ktPvVDJD71bT4CQROOcgPL7dbU0U1ndyCcb07pQXO7KWFUehSwXFPFeclbPuNzs4Zlu4JoQuITK+AgLzuD9i5/OOBgStjyGZpxnrIDCyYXjAmXhVxw6nqLjYqPoihhLPdR77mLRBSNzaKl0093YHg3lGvqsISkcGxfoczNzqDd5w9bPqWF1v0M1ZxPyL3Wy/IpCK31UfGJhYqPEka9bkpdm4/ugAlbCskwPteKHmvsGF7rp7GjC69HYiZAzcvOoDtgaLUDIcZmeikcmzlkWQ5ONVkbx4UEB6DUtoJONOmmcrFQ8VEiUPUZ7YRelFML+2H55FqWR2P78AYdNHZ0UzgmM2ZUXsgVF7J0sjO9FOdmDVnAwZG6dvKyM5gQ4a6cPn4sXo9w+Gz7kPSZ7qj4KGHU8lFO2uLTL8vHdnsN97xPU0cXhWNjr9kJhV+HxpTl9VCclz1kcz6H69qZOSG3lxBmZXiYPn4sh+virn8ftaj4KIoSZiDiE7Z8htvt1t5NUYz5HuixfELik53pYUJe1pBFux2pa2fWxNw+5bMm5Krl44CKjxL2tqnho5xs6iQ/O6NfWQBCkXFnh3in0GgaO7piBhtAhNstwvIpyc/hdHPnoC8paO7opqbxHHMn5fc5N7skj8N17QSC+tcVjYqPEkbdbsqxho5wNoBkmZCbTaZXONk8vBPrdW1dveZYIgntMxRp+cwpyaO9KzDoAQA7apoAuGRaYZ9zsyfm0uUPUtM4/Oug3I6KjxJGF5kqh8+2MXtiXr+u8XiEyQU5nGoavpDiLn+QujZfOKIsmtA2C2Hx8XqZP9myTPafbh3UsWw6XI/XIywsK+hzbkGpVbb7RPOg9jkSUPFR1N2mANa+N8caOmLOXSSitGAMp4bR8gkt3CwtiB2V1yfgIMPDXFt89p2H+FTXtvHwXw7xTr01jxMMGjbsPsXls4rJj+GqnDc5nyyvh101Kj7RJJe2VhnRhCweFaHRzbH6DoKGfls+AFMKcth8tHEIRhWb8KJOB/EJzVmFNnjLyvAwxuNlRvFYdhxvGlCfbT4/ax7ZSF2bj++9cIBPr5jNuJwMjtZ38MX3zo15TVaGhwtK89k5wD5HMio+Shj1uo1uDtlRWQOyfArHcKblFMGgweNJLhv2+RBKZ+Nk+YzN8pKT6aGzO0iGR/DaY1o+s5hn95wiEDThsmR5oeo0dW0+/vOOS3hx7xl+8OeDVpuzxnPzoimO1y0qK+R322rwB4JkeNXZFEKfhBJGN5Mb3Rw404oIzBqA5VNWNAZ/0HBqmPKYHbPdXk7BESISTo46JiKp5+Wzi2np9LP3ZEu/+9xxvIncLC83LSzlB3dcwlOfvpwf3HEJ//PxZXGFbNnM8bR3Bdil8z69UPFRVHIUAHbVNDNzQm7MRJ2JCLnqDtUOz4LK6to2phaOibmFdohQJFxeRPqdqyom4PUIz1Wd6nefO483sbCsICw0leXjueXiKQkzVl85ZwIi8PrBun73OZJR8VF63G2qQqOaPSeaWTi1b8RWMswpscSnepjE59DZ2Is6IynOsyyfyM3mJuRlc+WcCTy94yTBfqy96ewOsPdUC4unFfV7rONzs7hoSgF/PXi239eOZFR8lDCqPaOX2tZOTrd0Dlh8inOzKBybSfXZoRefQNBwKImQ8PF25oVoS+6DS6ZS03iOP++rTbrPt0+10B0wLJ42sOezYsJa3sgAABu+SURBVN5Etr7TOGS55dIRFR9F53oUdhyzorEWlfVdKJkMIsKciXkcPDO4a2hicbC2lY6uAItirKuJpCS/75wPwKqFpZQVjeGHLx9M2voJRasNxPIBuHnRFIIGNuzuv7tvpKLio4TdbhrtNnp583A92RkeLh7gN3uAi6YWsOdEC/4h3tdn2zuWECyZHl8IKiZZllFoQ7kQGV4Pn7+2gp01zTyx5XhSfe443sSkcdmOod2JmDc5n7mT8vjjzpMJ6xpjRsWCbxUfJYxaQKOXNw/Vc2n5eLIzBr7d85IZRZzrDpzXIs5k2HSkngl5WcwoHhu33iXTivAI3HJx3zDo25aWsXzWeB7809tU1yYe786aZi4eoFUYYvXiqWw+2sihOK7J7kCQm//zdW754Rt9RHOkoeKjhBkFX7aUGNS1+dh3upXLZxefVztLplsv563vDN1i0+5AkFf21bJiXknMfXwiKZ+Qy9avvZe7rijvc05E+N6HFpOT6eGTj20JL0aNRVNHF0fq2rk4Ru62/nB7ZRmZXuEXb77jWOe1A2epOtnC7hPN/HHnyHbRqfgoKjqjnJf2ngGsSfHzYWrhGKYWjuH16qELKX69uo6WTj/vXTApqfpFuVmOIjW1cAwP/91Salt9rHlko+N219uOWWK6dMbA5ntClOTncNPCUn67tYY2nz9mnZ3HmxCxtuP+3baa8+rP7aj4KGFUg0YnG/acZkbxWBaUjjuvdkSE6y4o4a8Hz3Kua2hcRj//21Em5GXznnklg9Le0hnjeewTyzjT3MmtP3qDXTV90+BsOdpIhkfO2+0GcNcV5bT6/DyxOfZc077TrcyckMuHKsvYeLh+2LepGE5UfJQwo2GSU+lNY3sXf6uuY+VFkxO6sZLh+gsn09kd5OV+hDEny/Zjjbyy/yx3Lp9OVsbgvbouLR/Pbz59BV6PcPtP3uS3W2t6/S28vK+WS6YXMiZr4PNhIZZML+KymeN55LVDdHb3Fej9Z1q5YPI4VtnRcc9VnT7vPt2Kio+iiUVHMb/Zehx/0PD+S6YOSnvLZxUztXAMv9h4dFDaC9HZHeCrv99DSX42n7x61qC2DbBgyjjW33sli6cV8uXf7OSzv95GXZuPLUcb2He6lVULSwetr89dW8GZFh+/iYq0a/f5eae+IxwZN6ckjz/tShwdl66o+CgqOqOUYNDwq03HuLS8iPmTz8/lFsLrET52+Qw2Hm5g6zsNg9KmPxDkK0/tYu+pFv7l/QsHlP4nGYrzsvn1p5Zz38r5vLj3DFf868vc+egmSgty+NCl0watnytmF7N0RhH/9eohuvw9YekH7DVS8ybnIyKsWljKpiMNcYMh0hkVH6UHVaFRxQt7z/BOfQd3Lp8xqO3euXwGk8Zl842nq3q9XAdCXZuPT/18C+t3nuS+lfO5LslAg4Hi9QifWTGb57/wLj52+Qzef0kZv/7U8rg55PqLiPC5ays42dzJE5uPhctDm9yFNr1btagUY+C5PSPT9abio/QsMlX1GTUEg4a1Lx5g1sTcQXUpgZVL7Vu3XEjVyRa+9ofd/cqhFuJcV4D/fu0w1699jTcO1fOdWy/iMytmD+o44zFrYh5fu3kB//qBhcyc0P8tJhLxrooJXDZzPP/nzwfDkW/7z7QyJtPLtCJr/dLcSflUlOTxp10jM+RaxUcJo/EGo4enttWw/0wrn7+2Ykj2mFl5USmfu2YOT26p4ZM/30JNY0fCa9p9fl7dX8tXntrJsgdf4sENb7OgdBx/vPeqQbfOUo2I8MBNF1DX1sUjrx0GrPxxcyfl9doPadWiUt462kDtMG1VMZzoZnKKWjyjjLo2H/+y4W0qZxTxvjiboJ0vX3zvXMbnZvFvz+3j3f/7Va6umEDljCImF4xhTKaX9i4/je1dHDrbxv4zbVSdaMYfNORmeVl5USl3LJtGZfn4IRtfqlk8rZBVi0r579cO88ElU9l+rImPXDa9V51VC0v5j5cO8uye0zEXy6YzKj5KGLV8Rj7BoOH+3+6m3efnXz+wcEh3HRUR/v7KmVy3YBK/3HiMF6pO8+r+vtsKTMjLZk5JLp961yyumF1M5YzxgxLWnA7cd8N8XtlXy3u//xpdgSArotYvVUzKZ96kfP6065SKjzLy6JnzUUY6D792mJfePsM3bl5AxaT8YemzrGgs9984n/tvnE+bz099m49z3QFyszIoHJtJfk7msIzDjUwvHsv3P7SYLz+5g+WzJ3L1nAl96qxaVMralw5wurlzwIlN3YiKjxJGF5mObJ7ecYJ/f34fqxaV8vEry1MyhrzsjCELlU5XVl40mfcuuAGPEHOh700LS/n+iwd4ds8pPn7lzBSMcGjQgANFLZ5RwLO7T/GlJ3dy2czxfO/2iwclm4EyeHg94vh/Mqckj/mT81mfxHYM6YSKj6JutxHO/7xxhH/49TYWlRXw07suJSdzdMynjCQ+uKSM7cea2He6JdVDGTSSEh8RWSki+0WkWkTuj3E+W0SesM9vEpHyiHMP2OX7ReSGRG2KyEy7jYN2m1nx+hCR94rIVhHZbf+8ZqAPY7SjXreRxbmuAA/8bhf/9Me9XHfBJH79yeXq8kpTPri0jKwMD7/edCxx5TQhofiIiBf4EXAjsAC4Q0QWRFW7G2g0xswB1gIP2dcuANYAFwIrgR+LiDdBmw8Ba40xFUCj3bZjH0Ad8D5jzELgLuAX/XsEito8I4+qk83c8sPXWbf5OJ9ZMZuf3Ll01ESQjUTG52axamEpv9t2gnaH7RjSjWQsn2VAtTHmsDGmC1gHrI6qsxp4zD5+CrhWLAfmamCdMcZnjDkCVNvtxWzTvuYauw3sNm+N14cxZrsxJuQMrQJyRCQ72QegRFo8KkLpTpvPzz8/s5dbfvgGjR3d/PwTy7hv5Xy8QxhSrQwPH71sOm0+f1JbcacDydjgU4HI9Ks1wGVOdYwxfhFpBort8o1R14bS58ZqsxhoMsb4Y9R36iNy56oPAtuNMX02wRCRe4B7AKZPnx59elQTkhx1u6Uv3YEgv91aw9qXDlDb6uOOZdP5yg3zKBybleqhKYPE0hlFzJuUz682HWPNsvR/hyUjPrG+MkW/ppzqOJXHsrji1U84DhG5EMsVd32MehhjHgEeAaisrNTXbASBoG6pkK4Eg4Y/7jrJ2hcPcLS+g4unFfJfdy5lyfTz23VTcR8iwkeXT+cbT1exq6aJRYOwuV0qSUZ8aoDIfOJlQLTdF6pTIyIZQAHQkODaWOV1QKGIZNjWT2R9pz4QkTLg98DHjDGHkrgnJYKBJH5UUsu5rgBPbT3Oz944ypG6duZNyue/P1bJdReUaBj1CObWS6byrxv28auNx1h0W3qLTzJzPpuBCjsKLQsrgGB9VJ31WJP9ALcBLxtrxeJ6YI0dqTYTqADecmrTvuYVuw3sNp+O14eIFAJ/Ah4wxrzRn5tXLAKhzeRUg1zP6eZOvvv8fi7/tz/z9aerGJeTwQ8/cgnPfv5q3rtgkgrPCGdcTiarF09h/c6TtHR2p3o450VCy8eeX7kXeB7wAj8zxlSJyLeBLcaY9cCjwC9EpBrLGlljX1slIk8CewE/8FljTAAgVpt2l/cB60TkO8B2u22c+gDuBeYAXxeRr9tl1xtjBn8f3xFKMLyTqaqPG+kOBHllXy1PbD7OK/trMcB1F0ziU1fP4tLyIhWcUcYdy6azbvNxnt5xkr9L42zfMhpTqlRWVpotW7akehiuYe5Xn6UrEOQ/77iE9108dFmOlf5x+GwbT22t4Tdbazjb6mNifja3Ly3jw5dOY0bx4O8xo6QHxhhu+sHrCPCnz101rF8+RGSrMaZyMNrSFWdK2O2mpJ5Tzed4Zucp1u88ye4TzXgErplfwocvnc575k0ckr13lPRCRPjIsml8/ekqdp9oTtvAAxUfRaPdUkxDexcbdluCs/loA8bAorICvrbqAt538RQmjRs5mYyVwWH1JVN5cMPbPP7WcRUfJf0ZjS7YVHG6uZMX957m+aozbDxcjz9omD0xly9cO5dbFk8Zkq2blZHDuJxMbl40hfU7TvC1VReQm4Zpk9JvxIqSplTXtvF81Wle2HuGncebAJg1IZdPXj2L911cyoLScRo8oCTNHcum8dTWGv6482RaLjpV8VGUISIQNOw43sRLb5/h+arTHD7bDsDFZQX8/zfM44YLJzF7Yp4KjjIglkwvYu6kPB7ffFzFR0lv1Ot2/tS3+Xjt4Fle3X+Wvxw4S1NHNxkeYfmsYv7+inLeu2ASpQVjUj1MZQQgIqy5dDrffmYve0+2sGDKuFQPqV+o+ChhdJ1P/wkGDbtPNPPq/rO8sr+WnTVNGAPFuVlcM7+EFfNKeHfFRArGjt6topWh4wNLpvJvz+7jN1uP880pF6Z6OP1CxUcJo5ZPctS3+XjjUD2v7q/lL/vPUt/ehQhcXFbIF66dy4p5E1k4tQCPZpJWhpjCsdaXnD/uPMlXb7ogrULxVXwUJQHnugJsPtrAG9V1/PVgHXtPWbtJFo7N5N1zJ/KeeSVcXTGB4jzdyUMZft6/ZCrPVZ3mr9V1vGdeSaqHkzQqPkoYtXwsAkHDnhPNvF5dx+sH69j6TiNdgSCZXmHpjCL+8fq5XDlnAovKCnWfHCXlvGdeCYVjM/n9thMqPsrwUtvSyUPP7ec7t17U790quwPB8PFo1R5jDO/Ud4TF5m+H6mjptLaUuqB0HHddMYMr50xg2czxjM3SPxnFXWRleLh5USlPba2hzedPm63S02OUSlzWvnSQ326rYfms8dxeOS3xBRH4/MHElUYYxhiO1new8XA9mw7Xs/FwA6dbOgGYUpDDyosmc1XFRK6YXcwEdaUpacD7L5nKLzce49ndp/r9DkgVKj4jgDZ7T/esjP5PNvq6A+HjkZrhwBjD4bp2W2wa2Hi4ntpWa7PbCXnZLJ81nstmFXPVnAmUF4/VdTdK2rFkehEzisfy++0nVHyU4eNclyUgOZn9c7kBdPpHntvNGMOhs228ebghbNnUtVliU5KfzfJZxVw2azzLZxUza0Kuio2S9ogIty6eyg9ePsjp5k4mF7g/H6CKzwigM8J66S8dttUEpK36+ANB9p1uZfPRBrYcbWTTkXrq2roAmDwuh6vmFHPZrGKWzypWy0YZsbzv4in8nz8f5Nk9p/j4lTNTPZyEqPiMANq7LAHpGsD8TfO59NsNsc3nZ/uxRrYcbWTrO41sP9ZIu239TSnI4eqKiSy3LZvp41VslNHBnJI85k/OZ8NuFR9lmOjwWS/eyMi1ZGnq6BEft2Y4ONV8ji1HG9lytIEt7zTy9qkWggZEYP7kcXxwaRlLZxRRWT6eqYWaukYZvdy0sJTvv3ggLVxvKj4jgJDlMyDxibB83BBvEAgaDpxpDQvNlqONnGg6B8CYTC+XTC/k3vfMobJ8PJdMLyQ/R9PWKEqIkPikg+tNxWcE0G7P23QF+q8eZ+wQY0jNlE9Hl58dx5vYerSRze80sv2dRlrt+ynJz+bS8vHcfdVMKsuLuKB0HJlplD5EUYabkOvtT7tUfJQhxhhDo+066x7AnM+x+o7BHlJcTjadY+s71lzNtmONVJ1sIRA0iMDcknzet3gKl5YXUTljPGVFY3S+RlH6yaqFpXwvDVxvKj5pTiiqC6BrAG63g7WtTBs/huMN5wbd7dYdCPL2qZaw2Gx9p5FTzZalNSbTy+JphXzm3bNZOqOIJdOLNPOzogwCNy2yxMftrjcVnzRnh70jJvTf8jnXFWDPiRZuWjjZEp/zdLw1dXSx7ZglMluONrKzponObmtMUwpyqCwfz9LphSydMZ75pfnqQlOUIWD2xPRwvan4pDkbdp8iN8tLe1eg3wEHz+45RVcgyHvml/CHHSf7bfmcaelko72I860j9Ryyd+r0eoQLp4zjjmXTw1bNFI1CU5RhIx1cbyo+acyLe8/whx0nuPvKmfx84zv9Cjho6ujiey8cYO6kPK6cMwGwIs3i0dDexV8PnuXNQ/VsOtLAkTpLbPKzM6gsL+IDS6yQ50VlBZqAU1FSSDq43vQNkYa0+fz88OVqHnntEBdNKeBL189l3ebjSVs+B8+08v89vp2zbT5+cMdyxthpeXz+vpkSjta188yuk7y8r5btx61dOsflZLBs5ng+smw6y2cVs2DKON1aQFFcRMj15uYFpyo+aURLZzePbzrGT18/wtlWH7ctLeNbt1zI2KwMMr2SMMPBwTOtPPzaYf6w/QTjxmTy049VsnRGEX5btHz2/EwwaHhh7xkeff0wm482ArCorIDPXVPBe+aXsHBqgYqNoricmxaWsvalA5xp6WTSOPe53lR8XE4waNh6rJHfbz/B+h0nafP5uXxWMQ//3VKWTC8K18v0emJaPrUtnazfeZI/7DjBnhMt5GR6uHP5DD77njlMzLe2C8jwevB6BJ8/yPGGDr7wxA62vtNIefFYvrJyHu+/ZCqlBTpnoyjpRHjB6e5T/L0LrR8VHxdyrivAxsP1/OXAWV7ce4YTTefIyfSw8sLJ3H3VLBaWFfS5JivDEw61Pt7QwYbdp3i+6nTYVbZwagFfv3kBqxdPiblHTXaGh4O1rXzwv/7Gue4A//7BRXxgydS02hNeUZQe5pTkMW9SPht2n1bxUWLT1NHFzppmdh5vYvPRBjYdaaDLHyQn08Pls4r5xxvmcv2CyeTG2aEwy+uhsb2Lz6/bzvqdVuTaRVPH8cXr5nLTwlLmlOTFHUNWhofnq84wLieD337mCuZOyh/s21QUZZi5aWEp//HnA9S2dFLiMtebis8w4/MH2HuyhZ3Hm9hxvImdNc3hqLHQKv+PLZ/Bu+dN5NLy8Unv0ZPp9fDK/rNkeIRPv3s2H1k2nWnjxyY9rhY7x9v9N16gwqMoI4SbFk5m7UsHeHbPae66ojzVw+mFis8QEgwajtS39wjN8Sb2nmqh2w6JLsnPZvG0Qm6vLGNxWSELywoGnCiz6ZyV6eAzK2bz5evn9fv6L1w3l6qTzay5ND12QVQUJTEVk/KpKMljw+5TKj4jmeZz3ew43sQ2O2/ZzuNNtHRaSTLHZnlZVFbAJ66aySXTCrl4WuGgTuLfcvEU/vx2Lf+wYs6Arv/ctRWDNhZFUdzDTQtL+cHLB6lt7aQk3z2uNzFuyKM/zFRWVpotW7YMSltVJ5t5fs9pXtl/lt0nmgHwCMydlM8l0wtZPK2QxdOKmFOSN+ThycGgwaMh0IqiRHDgTCvXr32Nf159IX93efl5tSUiW40xlYMxLrV8Bsgr+2tZ++IBdtU04xFYMr2IL143l8ryIi6eVkhenOCAoUKFR1GUaOZOymdOSR5/2n3qvMVnMFHx6SeBoOGb6/fwy43HKC8ey7duuZBbLp5CUW5WqoemKIoSk5sWlvJDl7neklrEISIrRWS/iFSLyP0xzmeLyBP2+U0iUh5x7gG7fL+I3JCoTRGZabdx0G4za6B9DAU/+cshfrnxGPe8axYvfPHd3HVFuQqPoiiu5paLSwka+M2WmlQPJUxC8RERL/Aj4EZgAXCHiCyIqnY30GiMmQOsBR6yr10ArAEuBFYCPxYRb4I2HwLWGmMqgEa77X730d8HkQyd3QF+8pdDXL9gEg/cOJ+sDF2AqSiK+5lTks+75k7kf/52lM7uvjkcU0Eyb89lQLUx5rAxpgtYB6yOqrMaeMw+fgq4VqwtKFcD64wxPmPMEaDabi9mm/Y119htYLd56wD7GHSO1rfT2unn5oun6A6biqKkFZ9dMZuzrT6+/JudqR4KkNycz1TgeMTnGuAypzrGGL+INAPFdvnGqGun2sex2iwGmowx/hj1B9JHGBG5B7jH/tgmIvVAneNdx2H1QwO5ytVMYIDPYgSiz8JCn0MPI+pZ/Bj48UcHdOkEYMZgjSMZ8Yn1FT86PtupjlN5LIsrXv2B9NG7wJhHgEdCn0Vky2CFDKY7+ix60Gdhoc+hB30WFvZzKB+s9pJxu9UAkcvey4CTTnVEJAMoABriXOtUXgcU2m1E99XfPhRFURSXkoz4bAYq7Ci0LKzJ/fVRddYDd9nHtwEvG2v16npgjR2pNhOoAN5yatO+5hW7Dew2nx5gH4qiKIpLSeh2s+dX7gWeB7zAz4wxVSLybWCLMWY98CjwCxGpxrJG1tjXVonIk8BewA981hgTAIjVpt3lfcA6EfkOsN1um4H0kYBHElcZNeiz6EGfhYU+hx70WVgM6nMYlel1FEVRlNSiC1UURVGUYUfFR1EURRl2RqX4JEoXNBIQkZ+JSK2I7IkoGy8iL9qpi14UkSK7XETkB/bz2CUiSyKuucuuf1BE7orVl5sRkWki8oqIvC0iVSLyebt8VD0LEckRkbdEZKf9HL5ll7s6ndVQYmdb2S4iz9ifR+WzEJGjIrJbRHaIyBa7bOj/Powxo+ofVoDDIWAWkAXsBBakelxDcJ/vApYAeyLK/h243z6+H3jIPr4JeBZrzdRyYJNdPh44bP8sso+LUn1v/XwOpcAS+zgfOICV0mlUPQv7fvLs40xgk31/TwJr7PKfAJ+xj/8B+Il9vAZ4wj5eYP/NZAMz7b8lb6rvb4DP5EvAr4Fn7M+j8lkAR4EJUWVD/vcxGi2fZNIFpT3GmNewogIjiUxRFJ266OfGYiPWWqtS4AbgRWNMgzGmEXgRK39e2mCMOWWM2WYftwJvY2XAGFXPwr6fNvtjpv3P4OJ0VkOJiJQBq4Cf2p9dndorBQz538doFJ9Y6YL6pOMZoUwyxpwC66UMlNjlTs9kRD0r211yCda3/lH3LGw30w6gFuvlcIgk01kBkems0vo52PwH8BUgaH9OOrUXI+9ZGOAFEdkqVhoyGIa/j9G4n09S6XhGGeeVuigdEJE84LfAF4wxLeKcGHbEPgtjrX9bLCKFwO+BC2JVs3+O2OcgIjcDtcaYrSKyIlQco+qIfxY2VxpjTopICfCiiOyLU3fQnsVotHxGczqeM7aJjP2z1i7vbxqktEJEMrGE51fGmN/ZxaPyWQAYY5qAV7F89qMxndWVwC0ichTL7X4NliU0Gp8FxpiT9s9arC8lyxiGv4/RKD7JpAsaqUSmKIpOXfQxO5JlOdBsm9rPA9eLSJEd7XK9XZY22L75R4G3jTHfjzg1qp6FiEy0LR5EZAxwHdb816hLZ2WMecAYU2asJJlrsO7to4zCZyEiuSKSHzrG+r3ew3D8faQ60iIV/7AiNg5g+by/murxDNE9Pg6cArqxvpXcjeWn/jNw0P453q4rWJv7HQJ2A5UR7XwCayK1Gvh4qu9rAM/hKizzfxeww/5302h7FsAirHRVu+yXyzfs8llYL8xq4DdAtl2eY3+uts/Pimjrq/bz2Q/cmOp7O8/nsoKeaLdR9yzse95p/6sKvQ+H4+9D0+soiqIow85odLspiqIoKUbFR1EURRl2VHwURVGUYUfFR1EURRl2VHwURVGUYUfFR1EURRl2VHwURVGUYef/AXsV0R4QSSFAAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], "source": [ "plt.clf()\n", "# plt.plot(x_part, calcs, '.')\n", "plt.plot(test_q, calcs_test, label = 'pdf')\n", "# plt.plot(test_q, res_y, label = 'res')\n", "plt.legend()\n", - "plt.ylim(0.0, 5e-6)\n", + "plt.ylim(0.0, 10e-6)\n", "# plt.yscale('log')\n", "# plt.xlim(3080, 3110)\n", "plt.savefig('test3.png')\n", @@ -1272,7 +1240,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 38, "metadata": {}, "outputs": [], "source": [ @@ -1286,7 +1254,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 39, "metadata": {}, "outputs": [], "source": [ @@ -1307,7 +1275,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 40, "metadata": {}, "outputs": [], "source": [ @@ -1317,10 +1285,60 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 41, "metadata": {}, "outputs": [], - "source": [] + "source": [ + " def jpsi_res(q):\n", + " return resonance(q, _mass = jpsi_mass, scale = jpsi_scale,\n", + " phase = jpsi_phase, width = jpsi_width)\n", + "\n", + " def psi2s_res(q):\n", + " return resonance(q, _mass = psi2s_mass, scale = psi2s_scale,\n", + " phase = psi2s_phase, width = psi2s_width)\n", + " \n", + " def p3770_res(q):\n", + " return resonance(q, _mass = p3770_mass, scale = p3770_scale,\n", + " phase = p3770_phase, width = p3770_width)\n", + " \n", + " def p4040_res(q):\n", + " return resonance(q, _mass = p4040_mass, scale = p4040_scale,\n", + " phase = p4040_phase, width = p4040_width)\n", + " \n", + " def p4160_res(q):\n", + " return resonance(q, _mass = p4160_mass, scale = p4160_scale,\n", + " phase = p4160_phase, width = p4160_width)\n", + " \n", + " def p4415_res(q):\n", + " return resonance(q, _mass = p4415_mass, scale = p4415_scale,\n", + " phase = p4415_phase, width = p4415_width)" + ] + }, + { + "cell_type": "code", + "execution_count": 59, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "2.404870995999285" + ] + }, + "execution_count": 59, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "_0 = jpsi_scale*np.cos(jpsi_phase)*jpsi_width/jpsi_mass**3 + psi2s_scale*np.cos(psi2s_phase)*psi2s_width/psi2s_mass**3\n", + "_1 = p3770_scale*np.cos(p3770_phase)*p3770_width/p3770_mass**3 + p4040_scale*np.cos(p4040_phase)*p4040_width/p4040_mass**3\n", + "_2 = p4160_scale*np.cos(p4160_phase)*p4160_width/p4160_mass**3 + p4415_scale*np.cos(p4415_phase)*p4415_width/p4415_mass**3\n", + "\n", + "R = (0.1/(1300**2) - ((_0 + _1 + _2)))\n", + "\n", + "R*10*2010**2" + ] }, { "cell_type": "code", diff --git a/raremodel-nb.py b/raremodel-nb.py index 1c50f0d..886b976 100644 --- a/raremodel-nb.py +++ b/raremodel-nb.py @@ -36,7 +36,7 @@ # In[2]: -# chunksize = 1000000 +# chunksize = 10000 # zfit.run.chunking.active = True # zfit.run.chunking.max_n_points = chunksize @@ -292,22 +292,19 @@ # In[5]: - - - # r = rho_scale * rho_width/rho_mass * np.cos(rho_phase)*(1-np.tan(rho_phase)*rho_width/rho_mass) # o = omega_scale*np.cos(omega_phase)*omega_width/omega_mass # p = phi_scale*np.cos(phi_phase)*phi_width/phi_mass -# phi_s = np.linspace(-500, 5000, 100000) +# # phi_s = np.linspace(-500, 5000, 100000) -# p_ = phi_s*np.cos(phi_phase)*phi_width/phi_mass +# # p_ = phi_s*np.cos(phi_phase)*phi_width/phi_mass -# p_y = r+o+p_ +# # p_y = r+o+p_ -# plt.plot(phi_s, p_y) +# # plt.plot(phi_s, p_y) -# # print(r + o + p) +# print(r + o + p) # ## Build pdf @@ -408,9 +405,10 @@ rho_mass, rho_width, rho_phase, rho_scale = pdg["rho"] -rho_m = zfit.Parameter("rho_m", ztf.constant(rho_mass), floating = False) +rho_m = zfit.Parameter("rho_m", ztf.constant(rho_mass), floating = False) #lower_limit = rho_mass - rho_width, +# upper_limit = rho_mass + rho_width) rho_w = zfit.Parameter("rho_w", ztf.constant(rho_width), floating = False) -rho_p = zfit.Parameter("rho_p", ztf.constant(rho_phase)) +rho_p = zfit.Parameter("rho_p", ztf.constant(rho_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi) rho_s = zfit.Parameter("rho_s", ztf.constant(rho_scale), floating = False) #omega @@ -419,7 +417,7 @@ omega_m = zfit.Parameter("omega_m", ztf.constant(omega_mass), floating = False) omega_w = zfit.Parameter("omega_w", ztf.constant(omega_width), floating = False) -omega_p = zfit.Parameter("omega_p", ztf.constant(omega_phase)) +omega_p = zfit.Parameter("omega_p", ztf.constant(omega_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi) omega_s = zfit.Parameter("omega_s", ztf.constant(omega_scale), floating = False) @@ -429,7 +427,7 @@ phi_m = zfit.Parameter("phi_m", ztf.constant(phi_mass), floating = False) phi_w = zfit.Parameter("phi_w", ztf.constant(phi_width), floating = False) -phi_p = zfit.Parameter("phi_p", ztf.constant(phi_phase)) +phi_p = zfit.Parameter("phi_p", ztf.constant(phi_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi) phi_s = zfit.Parameter("phi_s", ztf.constant(phi_scale), floating = False) #jpsi @@ -439,7 +437,7 @@ jpsi_m = zfit.Parameter("jpsi_m", ztf.constant(jpsi_mass), floating = False) jpsi_w = zfit.Parameter("jpsi_w", ztf.constant(jpsi_width), floating = False) -jpsi_p = zfit.Parameter("jpsi_p", ztf.constant(jpsi_phase)) +jpsi_p = zfit.Parameter("jpsi_p", ztf.constant(jpsi_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi) jpsi_s = zfit.Parameter("jpsi_s", ztf.constant(jpsi_scale), floating = False) #psi2s @@ -448,7 +446,7 @@ psi2s_m = zfit.Parameter("psi2s_m", ztf.constant(psi2s_mass), floating = False) psi2s_w = zfit.Parameter("psi2s_w", ztf.constant(psi2s_width), floating = False) -psi2s_p = zfit.Parameter("psi2s_p", ztf.constant(psi2s_phase)) +psi2s_p = zfit.Parameter("psi2s_p", ztf.constant(psi2s_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi) psi2s_s = zfit.Parameter("psi2s_s", ztf.constant(psi2s_scale), floating = False) #psi(3770) @@ -457,7 +455,7 @@ p3770_m = zfit.Parameter("p3770_m", ztf.constant(p3770_mass), floating = False) p3770_w = zfit.Parameter("p3770_w", ztf.constant(p3770_width), floating = False) -p3770_p = zfit.Parameter("p3770_p", ztf.constant(p3770_phase)) +p3770_p = zfit.Parameter("p3770_p", ztf.constant(p3770_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi) p3770_s = zfit.Parameter("p3770_s", ztf.constant(p3770_scale), floating = False) #psi(4040) @@ -466,7 +464,7 @@ p4040_m = zfit.Parameter("p4040_m", ztf.constant(p4040_mass), floating = False) p4040_w = zfit.Parameter("p4040_w", ztf.constant(p4040_width), floating = False) -p4040_p = zfit.Parameter("p4040_p", ztf.constant(p4040_phase)) +p4040_p = zfit.Parameter("p4040_p", ztf.constant(p4040_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi) p4040_s = zfit.Parameter("p4040_s", ztf.constant(p4040_scale), floating = False) #psi(4160) @@ -475,7 +473,7 @@ p4160_m = zfit.Parameter("p4160_m", ztf.constant(p4160_mass), floating = False) p4160_w = zfit.Parameter("p4160_w", ztf.constant(p4160_width), floating = False) -p4160_p = zfit.Parameter("p4160_p", ztf.constant(p4160_phase)) +p4160_p = zfit.Parameter("p4160_p", ztf.constant(p4160_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi) p4160_s = zfit.Parameter("p4160_s", ztf.constant(p4160_scale), floating = False) #psi(4415) @@ -484,7 +482,7 @@ p4415_m = zfit.Parameter("p4415_m", ztf.constant(p4415_mass), floating = False) p4415_w = zfit.Parameter("p4415_w", ztf.constant(p4415_width), floating = False) -p4415_p = zfit.Parameter("p4415_p", ztf.constant(p4415_phase)) +p4415_p = zfit.Parameter("p4415_p", ztf.constant(p4415_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi) p4415_s = zfit.Parameter("p4415_s", ztf.constant(p4415_scale), floating = False) @@ -507,6 +505,9 @@ # print(calcs_test) +# for param in total_f.get_dependents(): +# print(zfit.run(param)) + # ## Test if graphs actually work and compute values @@ -592,7 +593,7 @@ # In[14]: -# total_f.update_integration_options(draws_per_dim=2000000, mc_sampler=None) +total_f.update_integration_options(draws_per_dim=200000, mc_sampler=None) # inte = total_f.integrate(limits = (1000, 1040), norm_range=False) # inte_fl = zfit.run(inte) # print(inte_fl) @@ -840,16 +841,19 @@ tf.constant(0.93, dtype=dtype), tf.constant(0.05, dtype=dtype), tf.constant(0.065, dtype=dtype), + tf.constant(0.04, dtype=dtype), tf.constant(0.05, dtype=dtype)]), components_distribution=tfd.Uniform(low=[tf.constant(x_min, dtype=dtype), tf.constant(3090, dtype=dtype), tf.constant(3681, dtype=dtype), tf.constant(3070, dtype=dtype), + tf.constant(1000, dtype=dtype), tf.constant(3660, dtype=dtype)], high=[tf.constant(x_max, dtype=dtype), tf.constant(3102, dtype=dtype), tf.constant(3691, dtype=dtype), - tf.constant(3110, dtype=dtype), + tf.constant(3110, dtype=dtype), + tf.constant(1040, dtype=dtype), tf.constant(3710, dtype=dtype)])) # dtype = tf.float64 # mixture = tfd.MixtureSameFamily(mixture_distribution=tfd.Categorical(probs=[tf.constant(0.04, dtype=dtype), @@ -891,7 +895,7 @@ # In[26]: -0.00133/(0.00133+0.213+0.015)*(x_max-3750)/(x_max-x_min) +# 0.00133/(0.00133+0.213+0.015)*(x_max-3750)/(x_max-x_min) # In[27]: @@ -1036,26 +1040,34 @@ # # Fitting -# In[ ]: +# In[35]: start = time.time() +for param in total_f.get_dependents(): + param.randomize() + +# for param in total_f.get_dependents(): +# print(zfit.run(param)) + nll = zfit.loss.UnbinnedNLL(model=total_f, data=data2, fit_range = (x_min, x_max)) minimizer = zfit.minimize.MinuitMinimizer(verbosity = 5) # minimizer._use_tfgrad = False result = minimizer.minimize(nll) -param_errors = result.error() +# param_errors = result.error() -for var, errors in param_errors.items(): - print('{}: ^{{+{}}}_{{{}}}'.format(var.name, errors['upper'], errors['lower'])) +# for var, errors in param_errors.items(): +# print('{}: ^{{+{}}}_{{{}}}'.format(var.name, errors['upper'], errors['lower'])) print("Function minimum:", result.fmin) +# print("Results:", result.params) +print("Hesse errors:", result.hesse()) -# In[ ]: +# In[36]: print("Time taken for fitting: {}".format(display_time(int(time.time()-start)))) @@ -1066,7 +1078,7 @@ res_y = zfit.run(jpsi_res(test_q)) -# In[ ]: +# In[60]: plt.clf() @@ -1074,7 +1086,7 @@ plt.plot(test_q, calcs_test, label = 'pdf') # plt.plot(test_q, res_y, label = 'res') plt.legend() -plt.ylim(0.0, 5e-6) +plt.ylim(0.0, 10e-6) # plt.yscale('log') # plt.xlim(3080, 3110) plt.savefig('test3.png') @@ -1113,14 +1125,48 @@ # In[40]: -print((zfit.run(jpsi_p)%(2*np.pi))/np.pi) -print((zfit.run(psi2s_p)%(2*np.pi))/np.pi) +# print((zfit.run(jpsi_p)%(2*np.pi))/np.pi) +# print((zfit.run(psi2s_p)%(2*np.pi))/np.pi) -# In[ ]: +# In[41]: +def jpsi_res(q): + return resonance(q, _mass = jpsi_mass, scale = jpsi_scale, + phase = jpsi_phase, width = jpsi_width) +def psi2s_res(q): + return resonance(q, _mass = psi2s_mass, scale = psi2s_scale, + phase = psi2s_phase, width = psi2s_width) + +def p3770_res(q): + return resonance(q, _mass = p3770_mass, scale = p3770_scale, + phase = p3770_phase, width = p3770_width) + +def p4040_res(q): + return resonance(q, _mass = p4040_mass, scale = p4040_scale, + phase = p4040_phase, width = p4040_width) + +def p4160_res(q): + return resonance(q, _mass = p4160_mass, scale = p4160_scale, + phase = p4160_phase, width = p4160_width) + +def p4415_res(q): + return resonance(q, _mass = p4415_mass, scale = p4415_scale, + phase = p4415_phase, width = p4415_width) + + +# In[59]: + + +_0 = jpsi_scale*np.cos(jpsi_phase)*jpsi_width/jpsi_mass**3 + psi2s_scale*np.cos(psi2s_phase)*psi2s_width/psi2s_mass**3 +_1 = p3770_scale*np.cos(p3770_phase)*p3770_width/p3770_mass**3 + p4040_scale*np.cos(p4040_phase)*p4040_width/p4040_mass**3 +_2 = p4160_scale*np.cos(p4160_phase)*p4160_width/p4160_mass**3 + p4415_scale*np.cos(p4415_phase)*p4415_width/p4415_mass**3 + +R = (0.1/(1300**2) - ((_0 + _1 + _2))) + +R*10*2010**2 # In[ ]: diff --git a/test2.png b/test2.png index 38efb57..ef7d7d5 100644 --- a/test2.png +++ b/test2.png Binary files differ diff --git a/test3.png b/test3.png index 6a37a23..43cb0f4 100644 --- a/test3.png +++ b/test3.png Binary files differ