diff --git a/.ipynb_checkpoints/raremodel-nb-checkpoint.ipynb b/.ipynb_checkpoints/raremodel-nb-checkpoint.ipynb index 3980d05..4047084 100644 --- a/.ipynb_checkpoints/raremodel-nb-checkpoint.ipynb +++ b/.ipynb_checkpoints/raremodel-nb-checkpoint.ipynb @@ -16,7 +16,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\util\\execution.py:53: UserWarning: Not running on Linux. Determining available cpus for thread can failand be overestimated. Workaround (only if too many cpus are used):`zfit.run.set_n_cpu(your_cpu_number)`\n", + "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" ] }, @@ -151,16 +151,13 @@ "\n", " #Rotate by the phase\n", "\n", - " r = tf.abs(com)\n", + " r = ztf.to_complex(scale*tf.abs(com))\n", "\n", " _phase = tf.angle(com)\n", "\n", " _phase += phase\n", "\n", - " x = tf.cos(phase)*r\n", - " y = tf.sin(phase)*r\n", - "\n", - " com = tf.complex(scale* x, scale * y)\n", + " com = r * tf.exp(tf.complex(ztf.constant(0.0), _phase))\n", "\n", " return com\n", "\n", @@ -269,7 +266,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ @@ -324,15 +321,15 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "class total_pdf(zfit.pdf.ZPDF):\n", " _N_OBS = 1 # dimension, can be omitted\n", " _PARAMS = ['jpsi_mass', 'jpsi_scale', 'jpsi_phase', 'jpsi_width',\n", - " 'psi2s_mass', 'psi2s_scale', 'psi2s_phase', 'psi2s_width',\n", - " 'cusp_mass', 'sigma_L', 'sigma_R', 'cusp_scale'\n", + " 'psi2s_mass', 'psi2s_scale', 'psi2s_phase', 'psi2s_width'#,\n", + " #'cusp_mass', 'sigma_L', 'sigma_R', 'cusp_scale'\n", " ] # the name of the parameters\n", "\n", " def _unnormalized_pdf(self, x):\n", @@ -348,7 +345,7 @@ " def cusp(q):\n", " return bifur_gauss(q, mean = self.params['cusp_mass'], sigma_L = self.params['sigma_L'], sigma_R = self.params['sigma_R'], scale = self.params['cusp_scale'])\n", "\n", - " funcs = jpsi_res(x) + psi2s_res(x) + cusp(x)\n", + " funcs = jpsi_res(x) + psi2s_res(x) #+ cusp(x)\n", "\n", " vec_f = vec(x, funcs)\n", "\n", @@ -368,7 +365,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 5, "metadata": {}, "outputs": [], "source": [ @@ -396,7 +393,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 6, "metadata": {}, "outputs": [ { @@ -413,7 +410,7 @@ "#jpsi\n", "\n", "jpsi_mass, jpsi_width, jpsi_phase, jpsi_scale = pdg[\"jpsi\"]\n", - "jpsi_scale *= pdg[\"factor_jpsi\"]\n", + "# jpsi_scale *= pdg[\"factor_jpsi\"]\n", "\n", "jpsi_m = zfit.Parameter(\"jpsi_m\", ztf.constant(jpsi_mass), floating = False)\n", "jpsi_w = zfit.Parameter(\"jpsi_w\", ztf.constant(jpsi_width), floating = False)\n", @@ -423,7 +420,6 @@ "#psi2s\n", "\n", "psi2s_mass, psi2s_width, psi2s_phase, psi2s_scale = pdg[\"psi2s\"]\n", - "psi2s_scale *= pdg[\"factor_psi2s\"]\n", "\n", "psi2s_m = zfit.Parameter(\"psi2s_m\", ztf.constant(psi2s_mass), floating = False)\n", "psi2s_w = zfit.Parameter(\"psi2s_w\", ztf.constant(psi2s_width), floating = False)\n", @@ -432,12 +428,12 @@ "\n", "#cusp\n", "\n", - "cusp_mass, sigma_R, sigma_L, cusp_scale = 3550, 3e-7, 200, 0\n", + "# cusp_mass, sigma_R, sigma_L, cusp_scale = 3550, 3e-7, 200, 0\n", "\n", - "cusp_m = zfit.Parameter(\"cusp_m\", ztf.constant(cusp_mass), floating = False)\n", - "sig_L = zfit.Parameter(\"sig_L\", ztf.constant(sigma_L), floating = False)\n", - "sig_R = zfit.Parameter(\"sig_R\", ztf.constant(sigma_R), floating = False)\n", - "cusp_s = zfit.Parameter(\"cusp_s\", ztf.constant(cusp_scale), floating = False)" + "# cusp_m = zfit.Parameter(\"cusp_m\", ztf.constant(cusp_mass), floating = False)\n", + "# sig_L = zfit.Parameter(\"sig_L\", ztf.constant(sigma_L), floating = False)\n", + "# sig_R = zfit.Parameter(\"sig_R\", ztf.constant(sigma_R), floating = False)\n", + "# cusp_s = zfit.Parameter(\"cusp_s\", ztf.constant(cusp_scale), floating = False)" ] }, { @@ -449,13 +445,13 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "total_f = total_pdf(obs=obs, jpsi_mass = jpsi_m, jpsi_scale = jpsi_s, jpsi_phase = jpsi_p, jpsi_width = jpsi_w,\n", - " psi2s_mass = psi2s_m, psi2s_scale = psi2s_s, psi2s_phase = psi2s_p, psi2s_width = psi2s_w,\n", - " cusp_mass = cusp_m, sigma_L = sig_L, sigma_R = sig_R, cusp_scale = cusp_s)\n", + " psi2s_mass = psi2s_m, psi2s_scale = psi2s_s, psi2s_phase = psi2s_p, psi2s_width = psi2s_w)#,\n", + " #cusp_mass = cusp_m, sigma_L = sig_L, sigma_R = sig_R, cusp_scale = cusp_s)\n", "\n", "# print(total_pdf.obs)" ] @@ -469,7 +465,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 8, "metadata": {}, "outputs": [], "source": [ @@ -509,7 +505,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 9, "metadata": {}, "outputs": [ { @@ -521,7 +517,7 @@ }, { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAD8CAYAAACl69mTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8FdX5+PHPk4WEkBAghLAkkABhCbIIEXEtdSm44oISrYoLoha/Xfy1Vtpv1Vrtt7YWq9YNBcUVKHWJFhcUtypbkEXCGvZAgEBIAtmX5/fHHTDEhFwgYe7yvF+vvJg798yZZ4bkPvfMmXNGVBVjjDGmJYS4HYAxxpjAZUnGGGNMi7EkY4wxpsVYkjHGGNNiLMkYY4xpMZZkjDHGtBivkoyIjBaRdSKSIyL3NfB+hIjMct5fJCLJdd6b7KxfJyKjjqHOp0TkoDf7MMYY45uaTDIiEgo8DVwEpAHXiUhavWK3AftVtTfwOPCos20akAEMAEYDz4hIaFN1ikg60M6bfRhjjPFd3rRkhgM5qrpJVSuBmcCYemXGADOc5TnA+SIizvqZqlqhqpuBHKe+Rut0EtDfgHu93IcxxhgfFeZFmW7A9jqvc4HTGyujqtUiUgTEOesX1tu2m7PcWJ13A5mqmlcvhzS2j711C4nIRGAiQJs2bYb169fPi0M0xhyy92AFeUXlxLVpRdd2rd0O5wjVtcqavGIABnaLdTmawLV06dK9qhrfHHV5k2Qaai3Un4umsTKNrW+oBaUi0hW4Bhh5nHGgqlOBqQDp6emalZXVwGbGmMZM++9m/vT+am4+M5kHLx/gdjhH2HuwgvSHPwEg6y+XuBxN4BKRrc1VlzeXy3KBpDqvE4GdjZURkTAgFig4yraNrT8V6A3kiMgWIEpEcprYhzGmGfnyNeham2vR73iTZJYAqSKSIiKt8HTkZ9YrkwmMd5bHAvPVM/NmJpDh3BmWAqQCixurU1X/o6qdVTVZVZOBUqej/2j7MMY0I5/u6bS/eL/T5OUyp//jbuAjIBSYrqrZIvIQkKWqmcA04FWn1VGAJ2nglJsNrAaqgUmqWgPQUJ1NhNLgPowxLcMXv8PV+l5Ipgne9MmgqnOBufXW3V9nuRxPX0pD2z4CPOJNnQ2UifZmH8eiqqqK3NxcysvLT7Qq10VGRpKYmEh4eLjboZgAcqgh44uf5+qTUZmj8SrJBJLc3FxiYmJITk7Gn++AVlX27dtHbm4uKSkpbodjzElhLRn/E3TTypSXlxMXF+fXCQZARIiLiwuIFpkx3qq1LON3gi7JAH6fYA4JlOMwvuXQ75UPdskYPxSUScYY0zhf/u5itzD7H0syPurzzz/n0ksvBaCiooILLriAIUOGMGvWLJcjM8HCFzvZLcf4n6Dr+PdHy5Yto6qqiuXLl7sdigkCh+8u88EPdGvJ+B9rybhgy5Yt9OvXj/HjxzNo0CDGjh1LaWkpH374If369ePss8/mrbfeAmDPnj3ccMMNLF++nCFDhrBx40aXozfGPdbv73+CuiXzx/eyWb2zuFnrTOvalgcua3q+p3Xr1jFt2jTOOussbr31VqZMmcLzzz/P/Pnz6d27N+PGjQOgU6dOvPjiizz22GO8//77zRqrMQ059Dnui30zzdGSKausYc2uYtbmHWDrvhL2l1ZSWllD6/BQYluHkxLfhgFdYxnYLZbQEB88CX4mqJOMm5KSkjjrrLMAuOGGG3jyySdJSUkhNTX18LqpU6e6GaIJUr58Raqyuva4tss/UMF7K3by2bo9LNpccLieVmEhdIhqRetWoVRU1VBQWkl5lee99lHhjBrQmRtG9OAUm/H5uAV1kvGmxdFS6t9+XFRUZLckG5/gi9PJHFJV432SUVW+ztnHy99s4bN1e6ipVXp3iubGET0Y0TOO/l1i6Nau9RF/d7W1ys6iMr7dVshna/fw7vKdzFyynZF94/ndxf3pkxDTEocV0II6ybhp27ZtLFiwgDPOOIM333yTCy64gOeff56NGzfSq1cv3nzzTbdDNEHKd1OM9y2ZFdsLefTDtXyzcR/xMRFMOCeFa4Yl0rvT0ZNESIiQ2D6KxPZRXD64Kw9eXsWbi7fx9Gc5XPTEV9z94978z3m9CQu17mxvWZJxSf/+/ZkxYwZ33HEHqampPPHEEwwbNoxLLrmEjh07cvbZZ7Nq1Sq3wzRByIcbMlTVHD24jfkH+fvH65j73S7i2rTigcvSuP707kSEhR7X/mJbh3Pnj3oxLj2Jh95fzROfbuDrnL08f+Mw4qIjjqvOYGNJxiUhISE899xzR6wbPXo0a9eu/UHZkSNHMnLkyJMUmQl2PpxjqKiuObysqocvde0qKueJT9czOyuXyLAQfnlBKhPO6Ul0RPN8xLVv04rHxw1hZN947p2zkque/YaXbxlOSsc2zVJ/ILMkY4w5wqE+GfHBx5cVllYdXj5YUU1tLTzzRQ4vf72FWlVuHNGDu8/rTccWamWMGdKNxPZR3P5KFtdNXci/7jyDpA5RLbKvQGFJxgXJycl2Kcz4rEOXy3zxPpT9pZWHl//f7BUs3LSPAxXVXDmkG7+6sM9J+cAf1qM9r084nYypC7n+xYX8+84z6dQ2ssX366+CsvfKl++eORaBchzGt/jidDKHbNh9kIiwEHp3imbemt2c2asjc39+DlPGDTmpLYr+Xdryyq3D2XugkjtfW3rEZTxzpKBryURGRrJv3z6/n+7/0PNkIiPtG5RpXr763aWkoppP1+5mZN94nr5+KJU1tUS1cu8jbHBSOx67ZjCT3viWBzOz+b+rBrkWiy8LuiSTmJhIbm4u+fn5bodywg49GdOY5uSLOaayupZ756xkX0klt5/Tk7DQEJ+4jfiSQV1YtbMXz36+kXNT47loYBe3Q/I5XiUZERkNPAGEAi+q6l/qvR8BvAIMA/YB41R1i/PeZOA2oAb4uap+dLQ6RWQakI5nnr71wM2qelBEbgb+BuxwdvtPVX3xWA84PDzcniRpzFH4WkvmYEU1d722lK827OV3F/cjPbmD2yEd4Z4L+/B1zl5+9/Z3DEtuT6cYu7pQV5NfBUQkFHgauAhIA64TkbR6xW4D9qtqb+Bx4FFn2zQgAxgAjAaeEZHQJur8laoOVtVBwDbg7jr7maWqQ5yfY04wxpimHeqT8YWLyXsPVnDd1IV8s3Effx07iInn9nI7pB8IDw1hyrWDKa2s4Q/v2A099XnT3hwO5KjqJlWtBGYCY+qVGQPMcJbnAOeLp8NjDDBTVStUdTOQ49TXaJ2qWgzgbN8a32y9GxOwfKUls3VfCVc/+w0b9hzghZuGcW16ktshNap3pxh+fn4qH2Xv5vN1e9wOx6d4k2S6AdvrvM511jVYRlWrgSIg7ijbHrVOEXkJ2AX0A56qU+5qEVkpInNExHd/44wxJ2TVjiKufvYbisqqeOP2EZzXL8HtkJo04ZwUUjq24Y/vrba7zerwJsk01Gqu/12nsTLHut6zoHoL0BVYA4xzVr8HJDuX0T7h+5bTkYGITBSRLBHJCoTOfWNONrdvjf9qQz7jnl9ARFgoc+48k6Hd27saj7ciwkJ54LI0Nu8t4dUFW90Ox2d4k2RygbqthkRgZ2NlRCQMiAUKjrJtk3Wqag0wC7jaeb1PVSuct1/Ac5PBD6jqVFVNV9X0+Ph4Lw7PGFOXmznmnWU7uOWlJSR1iOKtn51J707R7gVzHEb27cQ5qR155vONlFRUux2OT/AmySwBUkUkRURa4enIz6xXJhMY7yyPBear5+tQJpAhIhEikgKkAosbq1M8esPhPpnLgLXO67r3Bl6Op5VjjGlmbuWYF77cxC9nLSc9uT2z7zyDBD8dRX/PhX0oKKnk5W+2uB2KT2jyFmZVrRaRu4GP8NxuPF1Vs0XkISBLVTOBacCrIpKDpwWT4WybLSKzgdVANTDJaaHQSJ0hwAwRaYvnktoK4C4nlJ+LyOVOPQXAzc1yBowxR/h+WpmTc39Zba3yyNw1TPvvZi4Z2IUp4wYf96zJvuDU7u05v18nnv9iIzee0YO2keFuh+Qqr8bJqOpcYG69dffXWS4Hrmlk20eAR7yssxY4q5F6JgOTvYnXGHP8Tua0MhXVNfz6Xyt5b8VObj4zmfsvTSMkAB55/MsL+nDZP//LrMXbuf3cnm6H4yr3h8waY3zKyeqTKa+qYeIrS3lvxU7uu6gfD1wWGAkGYGBiLCN6dmD615uP6WmegciSjDHmCIdyTEveZVZWWcOEGVl8uSGfv1w1kDt/1Muv5xJsyMRze5JXVM5/Vua5HYqrLMkYY47kJJeWSjElFdXc8vJivtm4l8fGDiZjePcW2pO7RvbpRK/4Nrz4302u3xbuJksyxpgjtOTH4YHyKsZPX8zizQU8Pm4IVw8L3AleQ0KEm89MZtWOYr7bUeR2OK6xJGOMOcKhL93N/eW7qKyKm6YvZtn2Qp66bihjhtSfOCTwjDm1G5HhIby5eHvThQOUJRljzBFqD18ua74sU1hayY3TFrFqRxFPXz+USwYFx5T4bSPDuWRgVzKX7wjawZmWZIwxR/i+47956isoqeT6FxaxNu8Az90wjNGndG6eiv3E9acnUVJZw/sr60+UEhwsyRhjjlDbjB3/ew9WcP0LC9mYf5AXxqdzfn/fn+iyuQ3t3p7enaKZszTX7VBcYUnGGHOEmhonyZxgltlTXE7G1IVs2VfC9JtP40d9gnMuQRHhiiFdWbJlPzsKy9wO56SzJGOMOUKNk11qa48/y+wq8iSYnYVlvHzLcM7q3bG5wvNLlw3uCsD7K4LvkpklGWPMEWqc5HK8I9V3FJYxbuoC9hyo4JVbhzOiZ1xzhueXesS1YXBSOzItyRhjgt3hJHMcLZntBaWMe34BBSWVvHrbcNKTOzR3eH7r8sFdyd5ZzMb8g26HclJZkjHGHOFQkqk+xpbMlr0ljHt+AQfKq3l9wumc6icPGztZLh3UBRF4f0VwTTNjScYYc4TjuVy2Mf8g46YuoKyqhjduP51Bie1aKjy/ldA2klOT2vHJmt1uh3JSWZIxxhzh+yTj3eWyDbsPkDF1ITW1ysyJZzCga2xLhufXLkzrzHc7isgrCp67zCzJGGOOcOjuMm9aMmt3FZMxdSEAMyeOoG/nmBaNzd9dmNYJgE9WB09rxpKMMeYI1Yf7ZI7eklm1o4jrpi4kPDSEWRNH0LuTJZim9IqPJqVjGz62JGOMCVaHOvwrj9KSWbG9kOtfWEhUqzBm3TGCnvHRJys8vyYiXJiWwMJN+ygur3I7nJPCkowx5ghlVZ7kUl3bcJJZunU/N7y4iNiocGZOHEGPuDYnMzy/d0H/BKpqlP9u2Ot2KCeFV0lGREaLyDoRyRGR+xp4P0JEZjnvLxKR5DrvTXbWrxORUU3VKSLTRGSFiKwUkTkiEt3UPowxzae8sgaA0oqaH7y3aNM+bpq2iA7RrZg18QySOkSd7PD83qnd2xEdEcZXlmQ8RCQUeBq4CEgDrhORtHrFbgP2q2pv4HHgUWfbNCADGACMBp4RkdAm6vyVqg5W1UHANuDuo+3DGNO8yqs9yWV/aeUR679cn8/4lxbTOTaSWRPPoGu71m6E5/fCQ0M4o1ccX23ID4onZnrTkhkO5KjqJlWtBGYCY+qVGQPMcJbnAOeL54HdY4CZqlqhqpuBHKe+RutU1WIAZ/vWfD8ZbGP7MMY0ozKnJVNYVnX4duZ5q3czYUYWKR2jmXXHGXSOjXQzRL93TmpHcveXsXVfqduhtDhvkkw3oO5j3XKddQ2WUdVqoAiIO8q2R61TRF4CdgH9gKea2McRRGSiiGSJSFZ+fr4Xh2eMqauozNMhrQp5RWW8smALd762lP5d2zLz9hF0jI5wN8AAcE6qZ0bqr3IC/5KZN0mmodZC/TZeY2WOdb1nQfUWoCuwBhh3DHGgqlNVNV1V0+Pjg3NqcWOOV22tsq+kkvP6ecZzjHt+Ife/m83IPvG8PuF0YqPCXY4wMCTHRdGtXWu+Wh/4X4S9STK5QFKd14lA/alED5cRkTAgFig4yrZN1qmqNcAs4Oom9mGMaSZ7SyqoqVV+1Cee64Z7/kR/O7ofU29KJzoizOXoAoeIcG6fjizYuO+Y54jzN9781iwBUkUkBdiBpyP/+nplMoHxwAJgLDBfVVVEMoE3RGQKnpZJKrAYT6vkB3U6fSy9VDXHWb4MWHu0fRzncRtjGrA27wAAqQnRjD8z2d1gAtyZvTry5uLtZO8sZnBS4M711mSSUdVqEbkb+AgIBaararaIPARkqWomMA14VURy8LQuMpxts0VkNrAaqAYmOS0UGqkzBJghIm3xJKIVwF1OKA3uwxjTfL7euJewELH5x06C01M8j0FYvLkgoJOMBHJjID09XbOystwOwxi/UFpZzY/+9jmndG3LS7cMdzucoDDyb5+RmhDDCzelux3KEURkqao2S1A24t8YA8CTn+aQf6CCu8/r7XYoQWN4SgeWbCk4oUdd+zpLMsYYPlu3h+e+2Mi49CSG9bCnWZ4sw1PiKCytYsOewH1apiUZY4Lc+t0H+OXM5fTrHMMfxwxwO5yg8n2/zD6XI2k5lmSMCWI7Csu4adpiWoWF8MJN6USGh7odUlBJbN+azm0jWbxlv9uhtBhLMsYEqYKSSm6atoiSimpm3DLcJrt0gYgwPKUDizfvC9h5zCzJGBOESiurufXlJWzfX8aL49NJ69rW7ZCC1rAe7dldXEFeUbnbobQISzLGBJnyqhrufO1bVuYW8tR1p3J6zx9MAWhOokNjZJZvL3Q5kpZhScaYIFJZXcvdb3zLl+vz+ctVgxg1oLPbIQW9/l1iaBUaYknGGOPfqmtq+cXMZXyyZg9/uuIUrj0tqemNTIuLCAslrWtbSzLGGP9VU6vcM3sFH6zaxR8uTePGET3cDsnUMSSpHd/lFgXkZJmWZIwJcLW1ym//vZLMFTv57eh+3HZ2itshmXqGJLWjrKqG9bsDb1CmJRljApiq8r/vrmLO0lx+dUEf7hrZy+2QTAOGBHDnvyUZYwKUqvLH91bzxqJt/GxkL35+vs1J5qt6xEXRPiqc5dsDb1CmJRljApCq8pcP1vLyN1uYcHYKvxnVF88jmowvEhEGJbZjxfYit0NpdpZkjAlAj89bz/NfbuKmM3rw+0v6W4LxA6d0a0tO/kHKq2rcDqVZWZIxJsD8c/4GnpyfQ8ZpSTx42QBLMH5iQNdYamqV9bsPuB1Ks7IkY0wAefqzHB77eD1XndqNP185kJAQSzD+YoAztU/2zmKXI2leTT5+2RjjH/45fwOPfbyeK0/txt+uGWwJxs8ktY8iJiKM7J2B1S/jVUtGREaLyDoRyRGR+xp4P0JEZjnvLxKR5DrvTXbWrxORUU3VKSKvO+tXich0EQl31o8UkSIRWe783H8iB25MIKmbYB67ZjChlmD8TkiI0L9r24BryTSZZEQkFHgauAhIA64TkbR6xW4D9qtqb+Bx4FFn2zQgAxgAjAaeEZHQJup8HegHDARaAxPq7OcrVR3i/Dx0PAdsTKB56tMNhy+RWYLxbwO6tmVt3gFqAuhxzN60ZIYDOaq6SVUrgZnAmHplxgAznOU5wPni6W0cA8xU1QpV3QzkOPU1WqeqzlUHsBhIPLFDNCZwPfnpBv4+z5Ng/mYJxu+ldWlLWVUNm/eWuB1Ks/EmyXQDttd5neusa7CMqlYDRUDcUbZtsk7nMtmNwId1Vp8hIitE5AMRafA5sSIyUUSyRCQrPz/fi8Mzxj89+ekGpsxbz1VDLcEEigFdYwECql/GmyTT0G9u/bZcY2WOdX1dzwBfqupXzutvgR6qOhh4CninoWBVdaqqpqtqenx8fENFjPF7T3xSJ8GMtQQTKFITomkVGsLqAOqX8SbJ5AJ15wRPBHY2VkZEwoBYoOAo2x61ThF5AIgH7jm0TlWLVfWgszwXCBeRjl7Eb0xAeeKTDTz+iSWYQBQeGkKfztGszguuJLMESBWRFBFphacjP7NemUxgvLM8Fpjv9KlkAhnO3WcpQCqefpZG6xSRCcAo4DpVPTzvtYh0dvp5EJHhTuz7juegjfFX//hkPY9/sp6rhyZagglQfRJiAmpAZpPjZFS1WkTuBj4CQoHpqpotIg8BWaqaCUwDXhWRHDwtmAxn22wRmQ2sBqqBSapaA9BQnc4unwO2AgucnPKWcyfZWOAuEakGyoAMJ5EZE/BUlSc+3cA/PtnA1UMT+evYQZZgAlTfhBje+nYHRaVVxEaFux3OCZNA/pxOT0/XrKwst8Mw5oSoKo99vI6nP9toCSYIfLZ2D7e8vIR/3XkGpyV3cCUGEVmqqunNUZdNK2OMD1NVHv7PGp7+bCPXDU/ib5ZgAl6fzjEArNsVGJfMbFoZY3xUba1yf+YqXlu4jZvPTOaBy9Jssssg0DU2kuiIMDYESL+MJRljfFBNrTL5rZXMzsrljh/15L7R/SzBBAkRITUhmnUBkmTscpkxPqa6ppZ7Zi9ndlYuvzg/1RJMEOqbEMOG3QfdDqNZWJIxxodUVtfyP28u493lO/nNqL786sI+lmCCUGpCDPtKKtl7sMLtUE6YJRljfER5VQ13vbaUD1bt4g+XpjHpx73dDsm4pG+Cp/N/fQB0/luSMcYHlFXWcPsrWXy6dg9/uuIUbjs7xe2QjIv6JEQDBMSgTOv4N8ZlJRXV3DZjCYs2F/DXsYO4Nj2p6Y1MQIuPiaBdVDjrAqBfxpKMMS4qLK3k5peW8N2OIv4xbghjhtSf4NwEIxGhV3w0m/L9P8nY5TJjXLKnuJxxzy9k9c5inv3pUEsw5gg9O7ZhUwA8V8aSjDEu2F5QytjnFrB9fykv33IaPxnQ2e2QjI/pGR9N/oEKDpRXuR3KCbEkY8xJtn73Aa5+9huKy6t44/YRnNnbnlhhfqhnfBsANuX7d2vGkowxJ9GK7YVc+/wCAGZNPIMhSe1cjsj4ql6Hksxe/+6XsY5/Y06Sbzbu5fYZWXSIbsXrt42ge1yU2yEZH9a9QxtCQ4SNe/y7JWNJxpiTYN7q3Ux641uS46J49bbTSWgb6XZIxse1CgshqX1ra8kYY47u7WW5/PpfKzmla1tevmU47du0cjsk4yd6xkdbn4wxpnEvf72ZX81awfDkDrx++whLMOaY9OzYhs17S6it9d+HS1qSMaYFqCp//XAtD763mgvTEnjpltOIjrALB+bY9IyPpqK6lh2FZW6Hcty8SjIiMlpE1olIjojc18D7ESIyy3l/kYgk13lvsrN+nYiMaqpOEXndWb9KRKaLSLizXkTkSaf8ShEZeiIHbkxLqaqp5TdzVvLM5xu5bnh3nv3pUCLDQ90Oy/ihw7cx+/GgzCaTjIiEAk8DFwFpwHUiklav2G3AflXtDTwOPOpsmwZkAAOA0cAzIhLaRJ2vA/2AgUBrYIKz/iIg1fmZCDx7PAdsTEsqraxm4itZzFnqeRbMn688hbBQu2Bgjs/3Y2X8t/Pfm9/+4UCOqm5S1UpgJjCmXpkxwAxneQ5wvngegjEGmKmqFaq6Gchx6mu0TlWdqw5gMZBYZx+vOG8tBNqJSJfjPG5jml1BSSXXv7CIL9bn8/AVp9izYMwJi4+OICYijM2B3JIBugHb67zOddY1WEZVq4EiIO4o2zZZp3OZ7Ebgw2OIAxGZKCJZIpKVn5/vxeEZc+Jy95cy9rlvWJ1XzDM/HcYNI3q4HZIJACJC97gotu4rdTuU4+ZNkmnoq1j9Wx0aK3Os6+t6BvhSVb86hjhQ1amqmq6q6fHx8Q1sYkzzWpNXzFXPfMPeAxW8dtvpjD7F5iEzzadHXBTbCgI7yeQCdR9wkQjsbKyMiIQBsUDBUbY9ap0i8gAQD9xzjHEYc1It3LSPa59fQIgI/7rzTIandHA7JBNgundoQ+7+Umr89DZmb5LMEiBVRFJEpBWejvzMemUygfHO8lhgvtOnkglkOHefpeDptF98tDpFZAIwCrhOVWvr7eMm5y6zEUCRquYdxzEb0yzmfpfHTdMX0ykmgn//7Ez6do5xOyQTgHrERVFVo+z009uYm7xxX1WrReRu4CMgFJiuqtki8hCQpaqZwDTgVRHJwdOCyXC2zRaR2cBqoBqYpKo1AA3V6ezyOWArsMDpNH1LVR8C5gIX47l5oBS4pTlOgDHHSlV54atN/HnuWob1aM+LN6XbIEvTYnp08Mxxt62glKQO/jffnVejw1R1Lp4P+brr7q+zXA5c08i2jwCPeFOns77BmJyW0SRv4jWmpVTX1PLge9m8tnAblwzswt+vHWxjYEyLOjSR6tZ9pZzV2+VgjoMNQTbGSyUV1dz9xrd8ti6fO37Uk9+O6kdIiN2ibFpWl9jWhIcKWwv88zZmSzLGeGF3cTm3vryENXnFPHzFKXaLsjlpQkOEpPZRbPPT25gtyRjThLW7irn1pSUUllUxbfxp/LhfJ7dDMkHGn8fK2HwXxhzFfzfs5ZpnF1Bdq8y+4wxLMMYVPTp4xsp4uqb9iyUZYxoxO2s7N7+0mK7tWvPOpLM4pVus2yGZINU9rg0HK6opKKl0O5RjZpfLjKmntlb5+7x1PP3ZRs5J7cjTPx1K28hwt8MyQezQbcxbC0qJi45wOZpjYy0ZY+ooqajmrteX8vRnG8k4LYnpN59mCca4rodzG7M/dv5bS8YYx87CMibMyGLtrmL+cGkat56VbLMoG59waBCmP3b+W5IxBli2bT+3v7KU8qoapt18Gj/uax38xndEhoeS0DaC7fstyRjjd95dvoPfzFlJQtsI3rj9dPok2Bxkxvckto9ix37/m7/MkowJWrW1yj8+Wc+T83MYntyB524cRgebg8z4qG7tWrNs+363wzhm1vFvglJpZTWT3viWJ+fncM2wRF6bcLolGOPTEtu3Jq+w3O+m/LeWjAk6OwvLmPhqFtk7i/n9xf2ZcE6KdfAbn9etfWuqa5XdxeV0bdfa7XC8ZknGBJVFm/Yx6Y1vKa+q5cWb0jm/f4LbIRnjlcT2njvMdhSW+VWSsctlJiioKq8u2MJPX1xE28hw3pl0piUY41e6OYkl18/uMLOWjAl4FdU13P9ONrOytnNev078I2OIDbA0ficovkUsAAAX20lEQVSxvSfJ+NsdZpZkTEDbVVTOna8tZfn2Qv7nvN786oI+9gwY45ciw0PpGN2KXEsyxviGpVsLuPO1bymtqOa5G4Yy+pQubodkzAnp1j6KHYX+lWS86pMRkdEisk5EckTkvgbejxCRWc77i0Qkuc57k53160RkVFN1isjdzjoVkY511o8UkSIRWe78HH78szH1vbFoGxlTF9KmVShvTzrLEowJCIntWgdeS0ZEQoGngQuBXGCJiGSq6uo6xW4D9qtqbxHJAB4FxolIGpABDAC6Ap+ISB9nm8bq/Bp4H/i8gXC+UtVLj+M4TZCoqK7hj++t5o1F2xjZN54nxp1KbJT1v5jAkNi+NfNW76a2Vv3msq83l8uGAzmquglARGYCY4C6SWYM8KCzPAf4p3gGHowBZqpqBbBZRHKc+misTlVd5qw7keMyQWhHYRk/e20pK3KLuGtkL379k76E+skfojHe6Na+NZU1tew9WEGntpFuh+MVby6XdQO213md66xrsIyqVgNFQNxRtvWmzoacISIrROQDERnQUAERmSgiWSKSlZ+f70WVJhB8sT6fS5/8ik35JTx3wzB+O7qfJRgTcA7dYbbdjy6ZeZNkGvpLrT+vQWNljnX90XwL9FDVwcBTwDsNFVLVqaqarqrp8fHxTVRp/N2h+cdufmkxCW0jyfyfsxl9Sme3wzKmRXRr9/2ATH/hzeWyXCCpzutEYGcjZXJFJAyIBQqa2LapOo+gqsV1lueKyDMi0lFV93pxDCYAFZRU8stZy/lyfT5XDe3GI1cMpHWrULfDMqbFdGvvfwMyvWnJLAFSRSRFRFrh6cjPrFcmExjvLI8F5quqOusznLvPUoBUYLGXdR5BRDo7/TyIyHAn9n3eHKQJPMu3F3LZU/9l4cZ9/PnKgfz9msGWYEzAi44Io11UuF8NyGyyJaOq1SJyN/AREApMV9VsEXkIyFLVTGAa8KrTsV+AJ2nglJuN5yaBamCSqtaA51bl+nU6638O3At0BlaKyFxVnYAned0lItVAGZDhJDITRFSV1xZt46H3skloG8m/7zqTgYmxbodlzEnTrV1rv7pcJoH8OZ2enq5ZWVluh2GaSUlFNf/7zireXraDH/eN5/FxQ2gXZdPzm+AyYUYWuftL+fCX57bYPkRkqaqmN0ddNuLf+IU1ecVMeuNbNu8t4f9d2IdJP+7tN+MEjGlOXWIjWbzZf3oKLMkYn6aqvLF4G398bzXtWofz+oTTObNXx6Y3NCZAdWkXSXF5NSUV1bSJ8P2PcN+P0ASt4vIqJr/1Hf9Zmce5feKZcu1gOkZHuB2WMa7qEusZhJlXVE7vTtEuR9M0SzLGJ63YXsj/vLmMHYVl/HZ0P+44t6ddHjMG6BLruY05r6jMkowxx0pVmfbfzTz64Vo6xUQy+44RDOvRwe2wjPEZXQ8nmXKXI/GOJRnjM/aXVPKbOSv4ZM0eLkxL4G9jB9ndY8bUkxDruWScV2hJxhivLd5cwC9mLmPvwQoeuCyNm89MtklSjWlARJjn4WW7iv1jrIwlGeOqqppanvhkA898nkNShyj+fdeZDEps53ZYxvi0LrGt2WktGWOObsveEn4xazkrthdyzbBEHrh8ANF+cEumMW7rHBvJ1n0lbofhFfuLNiedqvKvpbk8mJlNWIjw9PVDuWSQPbnSGG91jY1k4Sb/GJBpScacVIWllfzu7e+Y+90uRvTswJRrh9C1XWu3wzLGr3SObc2B8moOVlT7fOvft6MzAeWbjXv5f7NXkH+ggvsu6sft5/S0B4sZcxy6tvMMyNxVVEbvTjEuR3N0lmRMi6usrmXKvPU8/+VGUuLa8PbPzrKZk405AYcGZO4sLLckY4Lb+t0HuGf2clbtKOa64d35w6X9iWplv3bGnIjvp5bx/duY7a/dtIiaWmXafzfx2MfriYkI4/kbhzFqgD0W2ZjmkNA2EhH/GPVvScY0u637Svj1v1awZMt+Rg1I4JErB9rElsY0o1ZhIXSMjvCLUf+WZEyzUVVeX7SNP89dQ2iIMOXawVx5ajcbuW9MC+gSG0lesSUZEyR2FZVz779X8uX6fM5J7cijVw+yW5ONaUFdYiPZlO/7AzJDvCkkIqNFZJ2I5IjIfQ28HyEis5z3F4lIcp33Jjvr14nIqKbqFJG7nXUqIh3rrBcRedJ5b6WIDD3egzbNR1V5e1kuP3n8C5ZsLuBPYwbwyq3DLcEY08K6xLYOjD4ZEQkFngYuBHKBJSKSqaqr6xS7Ddivqr1FJAN4FBgnImlABjAA6Ap8IiJ9nG0aq/Nr4H3g83qhXASkOj+nA886/xqX7DtYwe/fXsWH2bsY1qM9f79mMMkd27gdljFBIaFtJAcrfH9ApjeRDQdyVHUTgIjMBMYAdZPMGOBBZ3kO8E/xXIgfA8xU1Qpgs4jkOPXRWJ2qusxZVz+OMcArqqrAQhFpJyJdVDXvWA7YnDhV5b2VeTyYmc3B8mobWGmMCzo7U/7vLi4nOt53H17mTZLpBmyv8zqXH7YgDpdR1WoRKQLinPUL623bzVluqk5v4ugGHJFkRGQiMBGge/fuTVRpjtWe4nJ+/84q5q3ezeDEWP46djB9O/v2YDBjAlFCW89Ymd3F5fTy8yTT0NdT9bJMY+sb6guqX+fxxIGqTgWmAqSnpzdVp/GSqjJnaS5/en81FdW1/O7iftx6VgphoV516xljmlndJOPLvEkyuUBSndeJwM5GyuSKSBgQCxQ0sW1TdR5PHKYF7CgsY/Jb3/Hl+nyGJ3fg0bGDSLG+F2NcdSjJ7CqqcDmSo/Pma+gSIFVEUkSkFZ6O/Mx6ZTKB8c7yWGC+03eSCWQ4d5+l4Om0X+xlnfVlAjc5d5mNAIqsP6Zl1dYqry3cyk+mfEHWlgIeGjOAmRNHWIIxxgdER4QRHRHm/y0Zp4/lbuAjIBSYrqrZIvIQkKWqmcA04FWnY78AT9LAKTcbz00C1cAkVa0Bz63K9et01v8cuBfoDKwUkbmqOgGYC1wM5AClwC3NdRLMD23dV8Jv/72ShZsKOLt3R/7vqoEkdYhyOyxjTB0JbSN8PsmIp8ERmNLT0zUrK8vtMPxKdU0t07/ezJR56wkPCeF/L+3PtelJNmrfGB90/QsLKa+q4a2fndWs9YrIUlVNb466fPfmanPSrdheyOS3vmN1XjEX9E/g4StOobMz26sxxvd0bhvJos0FbodxVJZkDAcrqnnso3W8smAL8TERPHfDUEYN6GytF2N8XEJsJHsOlFNbq4T46Dg1SzJB7uPsXTyQmc2u4nJuHNGDX4/qS9vIcLfDMsZ4ISEmgqoapaC00mdnOrckE6Tyisp44N1sPl69m36dY3j6p0MZ2r2922EZY47BocvZu4rKLckY31BTq7y6YAuPfbye6tpafju6HxPOSSHcBlUa43c6OWNl9hwoxzM80fdYkgkiq3YU8ft3VrFieyHnpHbkkSsG0j3Obks2xl919oMBmZZkgkBRWRVTPl7Hqwu30qFNK57IGMLlg7tax74xfi4+JgIR355axpJMAFNV3vp2B//3wRoKSiq5cUQP7vlJX2JbW8e+MYEgPDSEuDa+PSDTkkyAWpNXzP3vrmLJlv2c2r0dL98ynFO6+eY1W2PM8esca0nGnEQHyqt4fN4GZizYQmzrcP569SDGDkv02XvojTEnJiEmkp0+/IRMSzIBQlXJXLGTh/+zhr0HK7h+eHd+M6ov7aJauR2aMaYFJcRGsmx7odthNMqSTABYu6uYBzOzWbipgEGJsbx4UzqDk9q5HZYx5iRIiImkoKSSiuoaIsJC3Q7nByzJ+LGCkkoen7ee1xdtJSYynEeuPIWM07rbY5CNCSKHHsO8p7jCJ2dKtyTjh6pqanlt4VYen7eeksoabhzRg19e0If2bezSmDHBpu6ATEsy5oR9sT6fP72/mpw9BzkntSN/uDSNPgkxbodljHGJrw/ItCTjJzbvLeHh91fz6do9JMdF8cJN6VzQv5MNqDQmyB1KMr56G7MlGR9XXF7FP+fn8NLXm4kIC2XyRf24+axkn+zgM8acfO2iwmkVFmJJxhybqppa3ly8jSc/3cC+kkquGZbIr0f1pVOMPUTMGPM9ESGhbQS7fDTJeDX1roiMFpF1IpIjIvc18H6EiMxy3l8kIsl13pvsrF8nIqOaqlNEUpw6Njh1tnLW3ywi+SKy3PmZcCIH7qtUlfdX7uTCKV9w/7vZ9IqPJnPS2fx17GBLMMaYBiXERPpvS0ZEQoGngQuBXGCJiGSq6uo6xW4D9qtqbxHJAB4FxolIGpABDAC6Ap+ISB9nm8bqfBR4XFVnishzTt3POtvMUtW7T/CYfdY3G/fylw/WsjK3iL4JMbx082mM7Btv/S7GmKNKiI1k9c5it8NokDeXy4YDOaq6CUBEZgJjgLpJZgzwoLM8B/ineD4ZxwAzVbUC2CwiOU59NFSniKwBzgOud8rMcOo9lGQC0pq8Yv7ywVq+WJ9P19hIHrtmMFee2s3GuxhjvJIQE8lnxXtQVZ/7UupNkukGbK/zOhc4vbEyqlotIkVAnLN+Yb1tuznLDdUZBxSqanUD5QGuFpFzgfXAr1S1bh1+J3d/KVM+Xs/by3fQNjKc313cj5vOSCYy3Dr1jTHe6xwbQWllDQcqqn3u8eneJJmG0qJ6Waax9Q31BR2tPMB7wJuqWiEid+Jp5Zz3g2BFJgITAbp3795Ade7bX1LJM5/nMOObrSAw8dye/OxHvYmN8q1fDmOMf0g4NCCzuNwvk0wukFTndSKws5EyuSIShuc5oAVNbNvQ+r1AOxEJc1ozh8ur6r465V/A03fzA6o6FZgKkJ6eXj8Zuqq8qoaXvt7CM5/ncLCimrFDE/nVhX3o2q6126EZY/xYQp0Bmb07+dbgbG+SzBIgVURSgB14OvKvr1cmExgPLADGAvNVVUUkE3hDRKbg6fhPBRbjabH8oE5nm8+cOmY6db4LICJdVDXP2d/lwJrjPOaTrqZW+ffSXKbMW8+u4nLO79eJe0f3o29n3/plMMb4p8Oj/n3wDrMmk4zTx3I38BEQCkxX1WwReQjIUtVMYBrwqtOxX4AnaeCUm43nJoFqYJKq1gA0VKezy98CM0XkYWCZUzfAz0XkcqeeAuDmEz76FqaqfLpmD49+uJYNew4yJKkd/8gYwoiecW6HZowJIAk+POpfVH3qilKzSk9P16ysLFf2vXTrfh79YC2LtxSQ0rEN947qy+hTOvvcnR/GmMAw5KGPuXRQFx6+YuAJ1yUiS1U1vRnCshH/zW1j/kH+9uE6PszeRcfoCB6+4hTGnZZEeKhX416NMea4dIltTV6h77VkLMk0k70HK3jikw28sXgbkWEh3HNhH247O4U2EXaKjTEtr0tsJHk++Bhm+wQ8QWWVNUz/ejPPfr6Rsqoarh/enV9ckErH6Ai3QzPGBJEusZEs27bf7TB+wJLMCfgoexd/zMxmZ1E5F6YlcN9F/egVH+12WMaYINQlNpL9pVWUV9X41IBuSzLHYWdhGfe/m80na3bTr3MMU8bZHWPGGHd1ifWMt8srKielYxuXo/meJZlj9FH2Lu6ds5LK6lomX9SPW89OsU59Y4zrusR6bmPOKyqzJOOPVJUp89bz1PwcBnaL5anrTiXZh/4jjTHBrYszc4iv3WFmScYLqsr/vrOK1xdtY1x6Eg9dMcCeTGmM8Sm+OurfkowXnpqfw+uLtnHHj3py3+h+NqDSGONzWrcKpX1UODsLy9wO5QjWmdCEhZv2MWXeeq46tZslGGOMT+sc25pdPjZWxpLMUVTV1DL5re/o3iGKh688xRKMMcandY2NZKclGf/xzrIdbN5bwv2XphHVyq4sGmN8W5d2keQV2eUyvzH96y3079KW8/t3cjsUY4xpUpfY1hSWVlFWWeN2KIdZkmlEzp4DrMkr5tr0RLtMZozxC3XHyvgKSzKN+HTNHgAuHtjF5UiMMcY7h56ym7vfkozPW7p1P8lxUYcfBmSMMb7u0Ej/LftKXI7ke5ZkGrEit5BTu7d3OwxjjPFap5gIolqFsnmvJRmfVl5Vw+7iCp+a/8cYY5oiIiTHtbEk4+sOXc9M6tDa5UiMMebYpHRswxZ/SzIiMlpE1olIjojc18D7ESIyy3l/kYgk13lvsrN+nYiMaqpOEUlx6tjg1NmqqX00t4KSSgDio60/xhjjX5I7RrF9fxlVNbVuhwJ4kWREJBR4GrgISAOuE5G0esVuA/aram/gceBRZ9s0IAMYAIwGnhGR0CbqfBR4XFVTgf1O3Y3uoyWUVlYDnrmAjDHGn/Tv0paaWiV7Z7HboQDetWSGAzmquklVK4GZwJh6ZcYAM5zlOcD54hlcMgaYqaoVqroZyHHqa7BOZ5vznDpw6ryiiX00u/Iqz0CmKEsyxhg/c3qK5wGKH67a5XIkHt7MldIN2F7ndS5wemNlVLVaRIqAOGf9wnrbdnOWG6ozDihU1eoGyje2j711AxGRicBE5+VBEdlXv4y30lqsreSajhznuQhAdi487Dx8L6DOxeRHYfLxbdoR6NFccXiTZBpqLaiXZRpb31AL6mjlvY0DVZ0KTD0cmEiWqqY3sG3QsXPxPTsXHnYevmfnwsM5D8nNVZ83l8tygaQ6rxOBnY2VEZEwIBYoOMq2ja3fC7Rz6qi/r8b2YYwxxkd5k2SWAKnOXV+t8HTkZ9YrkwmMd5bHAvNVVZ31Gc6dYSlAKrC4sTqdbT5z6sCp890m9mGMMcZHNXm5zOn/uBv4CAgFpqtqtog8BGSpaiYwDXhVRHLwtC4ynG2zRWQ2sBqoBiapag1AQ3U6u/wtMFNEHgaWOXXT2D68MLXpIkHDzsX37Fx42Hn4np0Lj2Y9D2KNAWOMMS3FRvwbY4xpMZZkjDHGtJiATjJNTYfj70RkuojsEZFVddZ1EJF5zrQ880SkvbNeRORJ51ysFJGhdbYZ75TfICLjG9qXrxORJBH5TETWiEi2iPzCWR9U50NEIkVksYiscM7DH531xzxdU2NTQvkbZ5aRZSLyvvM6KM+FiGwRke9EZLmIZDnrWv7vQ1UD8gfPDQUbgZ5AK2AFkOZ2XM18jOcCQ4FVddb9FbjPWb4PeNRZvhj4AM94oxHAImd9B2CT8297Z7m928d2HOeiCzDUWY4B1uOZsiiozodzPNHOcjiwyDm+2UCGs/454C5n+WfAc85yBjDLWU5z/mYigBTnbynU7eM7znNyD/AG8L7zOijPBbAF6FhvXYv/fQRyS8ab6XD8mqp+yQ/HCtWdfqf+tDyvqMdCPOORugCjgHmqWqCq+4F5eOaZ8yuqmqeq3zrLB4A1eGaJCKrz4RzPQedluPOjHPt0TY1NCeVXRCQRuAR40Xl9PFNXBcS5aESL/30EcpJpaDqcbo2UDSQJqpoHng9eoJOzvrHzEXDnybnMcSqeb/FBdz6cy0PLgT14PgQ24uV0TUDdKaH8+jw4/gHcCxyaktjrqasIvHOhwMcislQ802/BSfj78GZaGX/l1TQ0QeRYp/7xSyISDfwb+KWqFkvjc6gG7PlQz1i0ISLSDngb6N9QMeffgD0PInIpsEdVl4rIyEOrGyga8OfCcZaq7hSRTsA8EVl7lLLNdi4CuSXjzXQ4gWi306zF+XePs/5Yp/jxOyISjifBvK6qbzmrg/Z8qGoh8Dmea+rHOl1TIJyHs4DLRWQLnsvl5+Fp2QTjuUBVdzr/7sHz5WM4J+HvI5CTjDfT4QSiutPv1J+W5ybnrpERQJHTPP4I+ImItHfuLPmJs86vONfOpwFrVHVKnbeC6nyISLzTgkFEWgMX4OmfOtbpmhqbEspvqOpkVU1Uz2SPGXiO7acE4bkQkTYiEnNoGc/v9SpOxt+H23c8tOQPnjsk1uO5Jv17t+NpgeN7E8gDqvB8w7gNzzXkT4ENzr8dnLKC50FxG4HvgPQ69dyKpzMzB7jF7eM6znNxNp5m+0pgufNzcbCdD2AQnumYVjofIvc763vi+WDMAf4FRDjrI53XOc77PevU9Xvn/KwDLnL72E7wvIzk+7vLgu5cOMe8wvnJPvR5eDL+PmxaGWOMMS0mkC+XGWOMcZklGWOMMS3GkowxxpgWY0nGGGNMi7EkY4wxpsVYkjHGGNNiLMkYY4xpMf8fCbWrKjzpHYAAAAAASUVORK5CYII=\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAD8CAYAAACl69mTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl4VdW5+PHvS+YQSCAJJJAAAQIYBBHCJNaLVQuOtBYrjsggXivXtvZXh3tvbeutvdp6tdWqiKKiBYGitqkTdawTUxgljGEOUyCEhBAynOT9/XE2EGJCDpCwz/B+nieP+6yz9trv2ZLzZu219tqiqhhjjDEtoZXbARhjjAlelmSMMca0GEsyxhhjWowlGWOMMS3GkowxxpgWY0nGGGNMi/EpyYjIaBHZICL5IvJgA+9Hichc5/3FItKtznsPOeUbRGTUabT5jIiU+XIMY4wx/qnJJCMiYcCzwJVAFnCTiGTVqzYJKFbVnsBTwOPOvlnAOKAvMBp4TkTCmmpTRLKBBF+OYYwxxn/50pMZAuSr6hZVrQLmAGPq1RkDzHS25wOXiYg45XNUtVJVtwL5TnuNtukkoD8A9/t4DGOMMX4q3Ic6nYGddV4XAEMbq6OqHhEpARKd8kX19u3sbDfW5lQgR1X31MshjR3jQN1KIjIFmALQunXrQX369PHhIxpjjtl/uJK9pRUkxUWRGh/tdjgn8dQq6/aUAtCvc7zL0QSvZcuWHVDV5OZoy5ck01Bvof5aNI3Vaay8oR6Uikgn4AZg5BnGgapOB6YDZGdna25ubgO7GWMaM+1fm3ns/fVMuaQ7/3nVeW6Hc5L9hysZ/OhHAOQ+drXL0QQvEdneXG35crmsAEiv8zoN2N1YHREJB+KBg6fYt7HyC4GeQL6IbANiRSS/iWMYY5qRP1+D1m//XWn8nC9JZimQKSIZIhKJdyA/p16dHGC8sz0W+ES9K2/mAOOcmWEZQCawpLE2VfVdVU1R1W6q2g0odwb6T3UMY0yosN/4gNPk5TJn/GMqsAAIA15W1TwReQTIVdUcYAbwutPrOIg3aeDUmwesBTzAPapaA9BQm02E0uAxjDGhw3JM4PFlTAZVfQ94r17Zw3W2K/COpTS076PAo7602UCdOF+OcTqqq6spKCigoqLibJtyXXR0NGlpaURERLgdiglC/nihwA9DMk3wKckEk4KCAtq0aUO3bt0I5BnQqkpRUREFBQVkZGS4HY4x54SNyQSekFtWpqKigsTExIBOMAAiQmJiYlD0yIzxlfVkAk/IJRkg4BPMMcHyOYx/8ed/VpZjAk9IJhljTGDyx3Eic2qWZPzUZ599xjXXXANAZWUll19+OQMGDGDu3LkuR2ZChT9+n/tjTObUQm7gPxCtWLGC6upqVq5c6XYoJgSIH9+OaUkm8FhPxgXbtm2jT58+jB8/nv79+zN27FjKy8v54IMP6NOnDxdffDFvvfUWAIWFhdx6662sXLmSAQMGsHnzZpejN6HCH7/PPbW1bodgTlNI92R+84881u4ubdY2szq15VfX9m2y3oYNG5gxYwYjRoxg4sSJPPnkk7zwwgt88skn9OzZkxtvvBGADh068NJLL/HEE0/wzjvvNGusxgSamlp/TH3mVEI6ybgpPT2dESNGAHDrrbfy9NNPk5GRQWZm5vGy6dOnuxmiMX7H0wxJ5nBFNSt2HGLjvsNs3n+EkqNVlFXWEBXeirbREXRPbs15qW0YmpFI6yj7ijxbIX0GfelxtJT6049LSkpsSrLxC8f+Gfrj+MeZ9mQKSyt4e8UuPl5XyPIdxceTVfvWkbSLjSAuKpxKTy2Hyqt5c3kBAJFhrfhOZhK3De/KJZnJtGplv59nIqSTjJt27NjBwoULGT58OG+88QaXX345L7zwAps3b6ZHjx688cYbbodoQtSx5OKPf/OcTk9GVfkqv4gZX27hXxv3U6twfue2TLmkOxf1SOK81DYkxkV9a7/DFdWsLijhk/WF5KzazR2vLCUrtS2/vCaL4T0Sm/PjhARLMi4577zzmDlzJnfddReZmZn86U9/YtCgQVx99dUkJSVx8cUXs2bNGrfDNCHMD3MMNT4M/KsqC/L28dxn+awuKCG5TRR3j+zBDwem0T05rsn920RHMKJnEiN6JvHA6D68s3o3//fPjdz04iJ+lJ3Gr6/rS2ykfXX6ys6US1q1asW0adNOKhs9ejTr16//Vt2RI0cycuTIcxSZMf6r0nPqJLNy5yF++85acrcXk5HUmseu78cPBnYmKjzsjI4XGd6K6wemcVW/VJ7+eBPP/2szuduLeXn8YLoltT6jNkONJRljzEmOLULpj5fLyio8ANQfHik+UsX/vLuWt5bvIikuiseu78cN2emENdM4SnREGPeP7sN3MpO5Z/Zyrn/+a2aMz+bCLu2apf1gZvfJuKBbt252Kcz4rRNjMv6XZcoqvUmmdZ3LVe99s4crnvoXOSt3c/fIHnz2i5GMG9Kl2RJMXcN7JPLm3RcRFxXObTOW8E1BSbMfI9iEZJIJlvWPguVzGOOrA2WVAERFhFFe5eG+uSv58azlpMRHkzP1Yh4Y3Ye4Fp52nJHUmrl3DSM+JoLbX15MfuHhFj1eoAu5JBMdHU1RUVHAf0Efe55MdHS026GYIOPPvxkb95UB3mTz/We/4u2Vu/jJZZn87ccjyOrU9pzFkRofw+w7hxLWqhWTZ+ZScrT6nB070ITcmExaWhoFBQXs37/f7VDO2rEnYxrTnPz1769D5VV8tG7f8df7D1fy2sQhfCcz2ZV4uia25vlbB3LT9EX8bO5KXro92+6laYBPSUZERgN/AsKAl1T1sXrvRwGvAYOAIuBGVd3mvPcQMAmoAe5V1QWnalNEZgDZeGdQbgTuUNUyEbkD+AOwyznsn1X1pdP9wBEREfYkSWNO4fjAv8tx1HXwSBUTXl3K4QoPsyYPZfHWg9w4OJ3OCTGuxjW4W3sevjaLh/+ex6tfb2PixfbdUl+Tl8tEJAx4FrgSyAJuEpGsetUmAcWq2hN4Cnjc2TcLGAf0BUYDz4lIWBNt/kxVL1DV/sAOYGqd48xV1QHOz2knGGNM0/ytJ7PzYDljn/+a9XtKef6WgYzomcR9V/RyPcEcc9uwrlzWpwOPf7CezfvL3A7H7/gyJjMEyFfVLapaBcwBxtSrMwaY6WzPBy4T79SUMcAcVa1U1a1AvtNeo22qaimAs38M/n2J2BjTgtbuLuX657/mQFklsyYP5Xt9U9wO6VtEhP+9vh8xkWH8fN4qam0Rz5P4kmQ6AzvrvC5wyhqso6oeoARIPMW+p2xTRF4B9gJ9gGfq1PuhiKwWkfkiku5D7MaY03R8UozL18sWbi7ixhcWEt5KmH/3RWR3a+9uQKfQoW00D1+Txcqdh5i/rMDtcPyKL0mmoX9q9VN1Y3VOt9y7oToB6ASsA250iv8BdHMuo33EiZ7TyYGITBGRXBHJDYbBfWPOtdrjOca9LPPeN3sY//ISOsZH8+bdF9GrYxvXYvHVDy7szKCu7fj9gvWUVthss2N8STIFQN1eQxqwu7E6IhIOxAMHT7Fvk22qag0wF/ih87pIVSudt1/EO8ngW1R1uqpmq2p2crI7s06MCWRuj8m8vnAb98xeTr+0eOb/+3A6+cnYS1NEhF9f25eiI1X8+ZN8t8PxG74kmaVApohkiEgk3oH8nHp1coDxzvZY4BP19rlzgHEiEiUiGUAmsKSxNsWrJxwfk7kWWO+8Tq1zvOvw9nKMMc3MrWVlVJUn/7mBX/49j8v6dOAvk4aSEBt5boM4S/3S4rn+wjRmfr2NwtIKt8PxC00mGWeMZSqwAO8X+zxVzRORR0TkOqfaDCBRRPKB+4AHnX3zgHnAWuAD4B5VrWmsTbyX0WaKyDfAN0Aq8IhzjHtFJE9EVgH3Anec9ac3xnyLGz2ZmlrlP99ew9Of5POj7DSm3TqImMgzW9TSbT+5LBNPrfLcZ/aodPDxPhlVfQ94r17Zw3W2K4AbGtn3UeBRH9usBUY00s5DwEO+xGuMOXPHcsy5SjYV1TX8ZM4KFuTt455Le/D/vtfbL9dN81WXxFhuGJTG7MU7uOvfupMaHxiX+1pKyC0rY4xpwjnsypQcreb2l5ewIG8fv7o2i1+M6hPQCeaYey7tSY0qM77Y6nYorrMkY4w5yfGeTAvfolZYWsGNLyxkxY5i/jRuABNGBM/d8untY7m6Xypzlu7kcIjPNLMkY4w5Se2xnkwL5pidB8sZO20hOw6WM2P8YMYMqH/rXeCb/J0Myio9zF26s+nKQcySjDHmJC19tWzrgSPc+MJCDpVXMWvyUC7pFZy3GvRPS2BIt/a88tU2PDVNPzY6WFmSMcacROv9tzlt3HeYH72wkApPLW9MGRb0T5aceHE3dh06ymcbQvfGcEsyxpiTHL9a1sxdmjW7Shg3fREAc6cMo2+n+GZt3x9ddl5HkuKimJsbupfMLMkYY07SEgP+K3YUc/OLi4gOb8W8u4aTGQDLxDSHiLBWjB2UxifrC0P25kxLMsaYkx3vyTRPc4u3FHHrS4tJiI1k3r8PJyOpdfM0HCBuHJxOTa0yf3loLpxpScYYc5Jjs8uaI8d8uekA419ZQkp8NPPuGk5au9hmaDWwZCS1ZmhGe/6aWxDwj30/E5ZkjDEn8TjLMJ/t9+GXmw4waeZSuiW2Zu5dw0mJj26G6ALTDwemsfXAEdbsKnU7lHPOkowx5iQ1zfDQra/yvQkmI6k1s+8cRlJcVDNEFrhG9U0hIkz4x+r6C9gHP0syxpiTHO/JnOEFs6/rJJhZk4fSvnVgraTcEuJjI7gkM5l3Vu0OuSdnWpIxxpykpubML5d9nX+AiTOX0rW9N8EkhngPpq5rL+jE7pIKlu8odjuUc8qSjDHmJDVnOBizcHMRE2cupUv7WGbdaQmmvsuzOhIV3oqcVaF1ycySjDHmJGcyJrNoSxETX11KertYG4NpRFxUOCN7J/PPvH0hNcvMkowx5iTHxmRqffwiXLSliAmvLCWtXYwlmCZcfl5H9pZWkLc7dGaZWZIxxpykpta7mKPHhx7NYifBdHYSTHIbSzCncmmfDojAR+v2uR3KOWNJxhhzEo8z8F/tOfXKwct3FDPh1aV0Sohm9p1DLcH4ICkuioFd2lmSMcaErqqapnsyebtLuOPlJSS3ieKNO4fRoU3o3mh5ui4/ryNrdpWyp+So26GcEz4lGREZLSIbRCRfRB5s4P0oEZnrvL9YRLrVee8hp3yDiIxqqk0RmSEiq0RktYjMF5G4po5hjGk+R6tqgBPJpr78wsPcPmMJcVHhzJo8lA5tLcGcjsvP6wDAx+sKXY7k3GgyyYhIGPAscCWQBdwkIln1qk0CilW1J/AU8LizbxYwDugLjAaeE5GwJtr8mapeoKr9gR3A1FMdwxjTvI5We5NMQ5fLdhSVc8tLixERZt05LCTXIjtbPTvEkdYuhs83hsYzZnzpyQwB8lV1i6pWAXOAMfXqjAFmOtvzgctERJzyOapaqapbgXynvUbbVNVSAGf/GE6s09fYMYwxzehYT6a6Xk9m96Gj3PzSIio9tcyaPDTkVlNuLiLCdzKTWLilKCSemOlLkukM1H3iToFT1mAdVfUAJUDiKfY9ZZsi8gqwF+gDPNPEMU4iIlNEJFdEcvfvD42/FIxpTuXHk8yJMZn9hyu59aXFlJRX8/rEofROCY3nwbSUET2TOFzhYfWuErdDaXG+JJmGegv1RwQbq3O65d4N1QlAJ2AdcONpxIGqTlfVbFXNTk4OzmeHG9NSVJWiI5XAiTGZQ+VV3DZjMXtKKnh5wmD6pQX/Ey1b2ogeSYh4V6oOdr4kmQIgvc7rNKD+ugjH64hIOBAPHDzFvk22qao1wFzgh00cwxjTTMoqPVRUe5NLeZWHwxXVjH95CVv2H+HF27MZ3K29yxEGh3atIzm/U7wlGcdSIFNEMkQkEu9Afk69OjnAeGd7LPCJetdNyAHGOTPDMoBMYEljbYpXTzg+JnMtsL6JYxhjmsmOg+WAdwmULfuPMOnVXPJ2l/LcLQO5ODPJ5eiCy4ieSSzfUcyRSo/bobSoJpOMM/4xFViA9/LVPFXNE5FHROQ6p9oMIFFE8oH7gAedffOAecBa4APgHlWtaaxNvJfEZorIN8A3QCrwyKmOYYxpPmucMYKr+6VSXlXDkm0H+eO4AVye1dHlyILPdzKT8NQqi7cWuR1Kiwr3pZKqvge8V6/s4TrbFcANjez7KPCoj23WAiMaaafRYxhjmseCvH2kxkfz0FV9iIkMY2TvZEb27uB2WEFpUNd2RIQJi7ce5Lt9gjeJ+5RkjDHBb82uEj7dUMg9I3uSEBvJr6/r63ZIQS06Ioz+aQks3RrcQ8u2rIwxhorqGh54czUJMRHceUl3t8MJGYO7teebXSVUODfABiNLMsaEOFXlN/9YS97uUp644QLiYyLcDilkDO7WjuoaZcWOQ26H0mIsyRgT4l74fAtvLNnB3SN7cNl5wTs24I+yu7ZHBJZuC95LZpZkjAlh83J38tj767n2gk784nu93Q4n5MTHRtC7YxtLMsaY4PPh2n089NY3fCczif+74QJatbKlAN0wuFt7lm8vDtp1zCzJGBOClmw9yNTZyzm/U1um3TqIyHD7KnDL4Iz2HKmqYd2ew26H0iLsX5YxIWblzkNMfHUpnRNiePmOwbSOsjsZ3HRhegIAK3cWuxxJy7AkY0wIWbOrhNtnLKZd6whm3TmUxDh7ZLLb0trFkBQXycqdwbkisyUZY0LE+r2l3DpjMW2iI5g9eRip8TFuh2TwPl/mgrQEVhUE5zRmSzLGhID8wsPc8uJiosJbMfvOoaS3tyda+pML0hPYvL+M0opqt0NpdpZkjAlyWw8c4eYXvY9Mnn3nMLom2hMt/c0F6QmowjcFwXfJzJKMMUFsR1E5N7+4CE+tMvvOofRIjnM7JNOAC5wHwa3cGXyXzCzJGBOkCorLuenFRRytruEvk4bSq6M9MtlfJcRGkpHUmlWWZIwxgWBvSQU3v7iY0opq/jJpKFmd2rodkmnCBWnxQTn4b0nGmCBTeLiCm19cxMEjVbw2cQjnd453OyTjgwvSE9hXWsm+0gq3Q2lWlmSMCSL7D1dyy4uL2VtawSsTBnNhl3Zuh2R81LeT94+BtbtLXY6keVmSMSZIFB6u4KYXF1FQfJQZ4wczuFt7t0Myp+G8VO+YWd7u4Jph5lOSEZHRIrJBRPJF5MEG3o8SkbnO+4tFpFud9x5yyjeIyKim2hSRWU75GhF5WUQinPKRIlIiIiudn4cxxgBQWFrBTdMXsav4KK9MGMzwHoluh2ROU5voCLomxrJ2T4j1ZEQkDHgWuBLIAm4Skax61SYBxaraE3gKeNzZNwsYB/QFRgPPiUhYE23OAvoA/YAYYHKd43yhqgOcn0fO5AMbE2wKSysY9+Ii9pRU8OqEwQzrbgkmUGWltg3Jy2VDgHxV3aKqVcAcYEy9OmOAmc72fOAyERGnfI6qVqrqViDfaa/RNlX1PXUAS4C0s/uIxgSvfaUVjJu+iL0lFbw6YQhDLcEEtKzUtmwrKqes0uN2KM3GlyTTGdhZ53WBU9ZgHVX1ACVA4in2bbJN5zLZbcAHdYqHi8gqEXlfRPo2FKyITBGRXBHJ3b9/vw8fz5jAtM+5RLavtIKZE4cwJMPGYALdsanm64LokpkvSaahJxmpj3VOt7yu54DPVfUL5/VyoKuqXgA8A/ytoWBVdbqqZqtqdnJyckNVjAl4e0u8PZhjCcYG+YNDMM4w8yXJFADpdV6nAbsbqyMi4UA8cPAU+56yTRH5FZAM3HesTFVLVbXM2X4PiBCRJB/iNyao7Ck5yrjpC9l/uJLXJg0h2xJM0OjYNor2rSNDLsksBTJFJENEIvEO5OfUq5MDjHe2xwKfOGMqOcA4Z/ZZBpCJd5yl0TZFZDIwCrhJVY8/j1REUpxxHkRkiBN70Zl8aGMC1e5DRxk3fREHyqqYOXEIg7paggkmIkJWalvy9gTPNOYmH4mnqh4RmQosAMKAl1U1T0QeAXJVNQeYAbwuIvl4ezDjnH3zRGQesBbwAPeoag1AQ206h5wGbAcWOjnlLWcm2VjgbhHxAEeBcU4iMyYkHEswxUeqeG3SEAbajZZBKatTW179ahuemlrCwwL/VkYJ5u/p7Oxszc3NdTsMY87azoPl3PzSIg4dqea1SUPsTv4gNn9ZAf/vr6v4+Of/5tqq2SKyTFWzm6OtwE+TxgS5rQeOcOMLCykpr+b1yUMtwQS5Xh29iWXTvsMuR9I8LMkY48c27TvMj15YSIWnljlThjMgPcHtkEwL69nBm2Q27itzOZLm0eSYjDHGHXm7S7htxhLCWglzpwwj054HExJiI8NJbx/DRuvJGGNaysqdh7hp+iKiw1sx767hlmBCTK8ObdgUJD0ZSzLG+Jml2w5y60uLSYiNZO5dw8lIau12SOYcy+zYhi0HyqiuqW26sp+zJGOMH/kq/wC3z1hCh7ZRzLtrOOntY90OybigV8c4qmuUbQeOuB3KWbMkY4yf+HR9IRNeXUrXxFjmThlOSny02yEZl/RyLo8Gw+C/JRlj/MAHa/Yy5fVcenWM4407h5HcJsrtkIyLenaIo5UQFIP/NrvMGJf9bcUufv7XVVyQFs8rE4YQHxPhdkjGZdERYXRpH8umwsBPMtaTMcZFry3cxk/nrmRIt/a8NmmoJRhzXGbHNna5zBhzZlSVZz7exMN/z+OKrI68MmEwcVF2YcGc0LNDHNuLjuAJ8BlmlmSMOcdqa5XfvruO//twI9cP7MzztwwkOiLM7bCMn8lIak11jVJQfNTtUM6KJRljziFPTS33v7maGV9u5Y6LuvHE2AuCYqVd0/x6JHvvj9pyILAvmdm/bmPOkYrqGu6ZvZz5ywr42eW9+NW1WbRq1dBDYo2BjCTvGmZb9gf2vTJ2EdiYc6Cs0sNdr+fyVX4Rv7o2iwkjMtwOyfi59q0jSYiNYEuA35BpScaYFlZ8pIo7Xl3Kml0lPPmjC7h+YJrbIZkAkZHUmq0B3pOxy2XGtKB9pRXcOH0h6/aUMu3WQZZgzGnpnhRnYzLGmIblF5Zx/XNfs6v4KK9OGMwVWR3dDskEmO7JrdlXWklZpcftUM6YT0lGREaLyAYRyReRBxt4P0pE5jrvLxaRbnXee8gp3yAio5pqU0RmOeVrRORlEYlwykVEnnbqrxaRgWfzwY1pScu2FzN22tdUemqZe9dwLuqR5HZIJgB1d1bgDuSFMptMMiISBjwLXAlkATeJSFa9apOAYlXtCTwFPO7smwWMA/oCo4HnRCSsiTZnAX2AfkAMMNkpvxLIdH6mAM+fyQc2pqV9vG4ft7y0iISYCN66+yLO7xzvdkgmQGU405g37w/cS2a+9GSGAPmqukVVq4A5wJh6dcYAM53t+cBlIiJO+RxVrVTVrUC+016jbarqe+oAlgBpdY7xmvPWIiBBRFLP8HMb0yLmLt3BlNeX0btjG+bffRFdEm2pfnPmuiW2RgS2BnNPBugM7KzzusApa7COqnqAEiDxFPs22aZzmew24IPTiAMRmSIiuSKSu3//fh8+njFnT1V5+uNNPPDmN1zcM4nZdw4jKc5WUjZnJzoijE7xMQF9r4wvSaahu8XUxzqnW17Xc8DnqvrFacSBqk5X1WxVzU5OTm5gF2OaV02t8t9/W8OTzjIxL43PprWtQ2aaSffk1gE9w8yX34QCIL3O6zRgdyN1CkQkHIgHDjaxb6NtisivgGTgrtOMw5hzqqK6hnvfWME/1+7j7pE9uH9Ub7xXio1pHt0SW7Ny5yFUNSD/bfnSk1kKZIpIhohE4h3Iz6lXJwcY72yPBT5xxlRygHHO7LMMvIP2S07VpohMBkYBN6lqbb1j3O7MMhsGlKjqnjP4zMY0i5Lyam6bsZgP1+3jV9dm8cDoPgH5JWD8W9fEWA5XeDhUXu12KGekyZ6MqnpEZCqwAAgDXlbVPBF5BMhV1RxgBvC6iOTj7cGMc/bNE5F5wFrAA9yjqjUADbXpHHIasB1Y6PzCvqWqjwDvAVfhnTxQDkxojhNgzJnYUVTOHa8uoeDgUZ656UKu6d/J7ZBMkEpv7508suNgOe1aR7oczenz6cKxqr6H90u+btnDdbYrgBsa2fdR4FFf2nTKG4zJ6Rnd40u8xrSkFTuKmTwzF0+t8vqkIQztnuh2SCaIdXVmKG4/WM4F6QkuR3P6bHTSmNPw/jd7+OnclXRsG80rEwbTIznO7ZBMkOvi9GR2Hix3OZIzY0nGGB+oKi99sZXfvb+OAekJvHR7Nok2RdmcA7GR4STFRbG9KDCnMVuSMaYJnppafv2PPP6yaAdX9UvhyR8NsCdZmnOqa2Is24usJ2NM0DlS6WHq7OV8umE/d/1bdx4Y1cceNGbOua7tY1m0pcjtMM6IrcJsTCP2llRww7SFfL7pAL/7QT8euvI8SzDGFentY9lTWkGlp8btUE6b9WSMacDa3aVMmrmU0qPVzBifzcjeHdwOyYSwromxqMLOg0fp2SGwJptYT8aYehbk7WXstK8B+Ou/X2QJxrju2DTmQJxhZj0ZYxyqynOfbeYPCzYwID2B6bcNokPbaLfDMub4DZmBOMPMkowxeNcge/DN1fxt5W7GDOjE4z/sbzPIjN9IjosiNjKMHQePuh3KabMkY0Je4eEK7np9GSt2HOIXo3rz45E9bA0y41dEhC7tY9lx0HoyxgSUNbtKmPJaLsXl1Uy7dRCjz09xOyRjGpTePjYgH8NsA/8mZH2wZg83TFsIwPy7h1uCMX4tvV0suw4dxbuMY+CwnowJOarKs5/m88Q/N3oH+G8fRIc2NsBv/FvndjGUV9VQXF5N+wBajdmSjAkpRyo93D9/Ne9+s4fvD+jEYzbAbwJE54QYAHYVH7UkY4w/2l50hCmvLWNT4WEeurIPUy7pbgP8JmCktXOSzKFy+qXFuxyN7yzJmJDw2YZC7n1jBSLCzIlD+E5mstshGXNajvVkCooDaxqzJRkT1FSV5//lvcGyd8c2TL8tmy7O3dPGBJKE2AhiI8PYdciSjDF+oe74yzVaZifVAAAXXElEQVT9U/n92P7ERto/eROYRIS0djEB15PxaQqziIwWkQ0iki8iDzbwfpSIzHXeXywi3eq895BTvkFERjXVpohMdcpURJLqlI8UkRIRWen8HH/8szH1bS86wvXPfc37a/bw0JV9eOamCy3BmIDXOSGGXQGWZJr8rRORMOBZ4AqgAFgqIjmqurZOtUlAsar2FJFxwOPAjSKSBYwD+gKdgI9EpJezT2NtfgW8A3zWQDhfqOo1Z/A5TQix8RcTrDq3i2H5jkNuh3FafOnJDAHyVXWLqlYBc4Ax9eqMAWY62/OBy8Q7bWcMMEdVK1V1K5DvtNdom6q6QlW3neXnMiGotlZ5+uNNTHh1KZ0SYvjH1IstwZig0jkhlpKj1ZRVetwOxWe+JJnOwM46rwucsgbrqKoHKAEST7GvL202ZLiIrBKR90Wkb0MVRGSKiOSKSO7+/ft9aNIEg+IjVUx4dSlPfriRMRd04q0fX2QD/CbodG534l6ZQOHLReqGbiSov65BY3UaK28ouTW1VsJyoKuqlonIVcDfgMxvNaI6HZgOkJ2dHVjrL5gzsnLnIe6ZtZz9hyv57ffP55ahXez+FxOUTkxjLqd3ShuXo/GNLz2ZAiC9zus0YHdjdUQkHIgHDp5iX1/aPImqlqpqmbP9HhBRd2KACT2qyusLt3GD84Cx+XcP59ZhXS3BmKCVfvyGzMDpyfiSZJYCmSKSISKReAfyc+rVyQHGO9tjgU/Uu4pbDjDOmX2WgbfnscTHNk8iIinOOA8iMsSJvciXD2mCz5FKDz+du5Jf/j2Pi3sm8e69F9M/LcHtsIxpUUlxUUSGtQquy2Wq6hGRqcACIAx4WVXzROQRIFdVc4AZwOsiko+3BzPO2TdPROYBawEPcI+q1oB3qnL9Np3ye4H7gRRgtYi8p6qT8Savu0XEAxwFxmmgLUdqmkV+4WHu/styNu8v4/99rxc/HtmTVq2s92KCX6tWQqeEaAoCqCcjwfw9nZ2drbm5uW6HYZrRP1bt5oE3VxMTEcbTN13IiJ52xdSEllteWsSRyhr+ds+IFjuGiCxT1ezmaMvuTjMBoaK6hv95Zy2zFu9gUNd2/PnmC0mNj3E7LGPOuc4JMXyyPnBmzlqSMX4vv7CMqbOXs37vYaZc0p1fjOpNRJg9b8+Eps4JsRwoq6TSU0NUuP8/psKSjPFbqsr8ZQU8/Pc8YiLDeGXCYC7t3cHtsIxxVWq89wF7haWVpLf3/3vBLMkYv1RW6eGXf1vD2yt2Max7e/5444WkxNvTK41JTfD+Huw+dNSSjDFnIm93CVNnr2B70RF+enkm//HdTMJs9pgxwImezN7SCpcj8Y0lGeM3VJXXFm7n0XfX0a51BLPvHMaw7oluh2WMX0lxJrzsKbEkY4zPDpVX8cCbq1mQt49LeyfzxA0XkBgX5XZYxviduKhw2kSHsydA7pWxJGNc9/XmA9w3dxUHyir576vPY+KIDLu50phTSI2Ptp6MMU2p8tTyf//cwPQvtpCR2Jq3fnyRLQ1jjA9S42NsTMaYU8kvPMxP5qwkb3cpNw/twn9ffZ49udIYH6XGR5O3u9TtMHxiv9XmnFJV/rJ4B4++u5bYyHCm3zaI7/VNcTssYwJKSnw0B8oqqfLUEhnu3zcmW5Ix58yBskoemL+aj9cXckmvZJ4Y258Obe3eF2NOVydnhtm+0gq/v1fGkow5Jz7dUMgv/rqK0goPv7o2i/HDu9ngvjFn6NiNyXtKLMmYEFde5eHx99czc+F2+qS0YdbkYQHzRD9j/FWnhGNJxv+nMVuSMS1m2faD/HzeKrYVlTNxRAb3j+5NdIT/L+hnjL8LpBsyLcmYZlfpqeGpDzcx/fPNpMbH8Madwxjew+7cN6a5HLshc68lGRNq1uwq4efzVrFh32FuGpLOf12dRVyU/TMzprl5b8i0y2UmRFTX1PLcp5t55pNNtG8dySt3DObSPrYsvzEtJSU+JiAul/k0wVpERovIBhHJF5EHG3g/SkTmOu8vFpFudd57yCnfICKjmmpTRKY6ZSoiSXXKRUSedt5bLSIDz/RDm+a1ad9hrn/ua576aCNX90/lnz+7xBKMMS2sU4AsLdNkT0ZEwoBngSuAAmCpiOSo6to61SYBxaraU0TGAY8DN4pIFjAO6At0Aj4SkV7OPo21+RXwDvBZvVCuBDKdn6HA885/jUtqapUZX27hiX9uJC4qnOduGchV/VLdDsuYkBAoN2T6crlsCJCvqlsARGQOMAaom2TGAL92tucDfxYRccrnqGolsFVE8p32aKxNVV3hlNWPYwzwmqoqsEhEEkQkVVX3nM4HNs0jv/AwD7z5Dcu2F3NFVkd+94N+JLexVZONOVc6xceg6v83ZPqSZDoDO+u8LuDbPYjjdVTVIyIlQKJTvqjevp2d7aba9CWOzsBJSUZEpgBTALp06dJEk+Z0VdfUMv3zLfzpo03ERoXx5I8u4AcXdm7ojwJjTAtKqfPwskBPMg19e6iPdRorb6hvV7/NM4kDVZ0OTAfIzs5uqk1zGvJ2l3D//NXk7S7lqn4p/Oa68633YoxLjj0hc7efP1fGlyRTAKTXeZ0G7G6kToGIhAPxwMEm9m2qzTOJw7SASk8Nz3ycz7R/bSYhNpJptw5k9Pk29mKMm471ZPb5+ZL/vowWLQUyRSRDRCLxDuTn1KuTA4x3tscCnzhjJznAOGf2WQbeQfslPrZZXw5wuzPLbBhQYuMxLW/Z9mKufvpL/vxpPtcN6MRH911iCcYYP9AmOoK4qHC/n2HWZE/GGWOZCiwAwoCXVTVPRB4BclU1B5gBvO4M7B/EmzRw6s3DO0nAA9yjqjXgnapcv02n/F7gfiAFWC0i76nqZOA94CogHygHJjTXSTDfVl7l4YkFG3nl662kto3mlQmDubS3TUs2xp+kxEf7/V3/4u1wBKfs7GzNzc11O4yA8+WmA/zn29+w42A5tw7rwgOj+9AmOsLtsIwx9dw2YzFllR7e/vGIZm1XRJapanZztGV3/JvjDpRV8ui763h7xS66JcYyZ8owhnW3NceM8VcpbaP5Mv+A22GckiUZg6ry19wCfvf+Oo5Uerj3uz358aU9bcVkY/xcSnw0hYcr8dTUEh7mnzdkWpIJcfmFZfzn29+wZOtBBndrx+9+0I/Mjva8F2MCQUp8NDW1yoGyquOzzfyNJZkQVVFdw3Ofbeb5z/KJiQjjsev78aPsdHtapTEBJDX+xMPLLMkYv/H15gP899tr2HLgCN8f0In/ujrLbqo0JgCltPU+vMyfZ5hZkgkhRWWV/O699by5vICuibG8PmkI38lMdjssY8wZOtGTsSRjXFRTq8xesoMnFmzgSKWHey7twX98N9MG9o0JcAmxEUSFt/Lru/4tyQS5FTuK+eXf17BmVynDuyfyyJi+NrBvTJAQEVL8/LkylmSCVFFZJb//YANzc3fSsW0Uz9x0Idf0T7XVko0JMilt/fuuf0syQab+pbG7LunOf1yWSVyU/a82JhilxkezbEex22E0yr55gohdGjMm9KTEx7CvZC+1teqXtyBYkgkCRWWV/GHBBuYs9V4ae/qmC7nWLo0ZExJS46OpqqnlYHkVSXH+dyuCJZkAVuWp5bWF2/jTx5s4WlXDlEu6c69dGjMmpBx/QmZJhSUZ0zxUlU/WF/Lou+vYcuAI/9YrmV9ecx49O9ilMWNCTUrbE0nm/M7xLkfzbZZkAsymfYf5n3fX8fnG/XRPbs0rdwzm0j72nBdjQtXxGzL99F4ZSzIB4lB5FX/8aBOvL9pObGQYv7wmi9uHdyXCT1deNcacG4lxUYS3EvaWHHU7lAZZkvFznppaZi3ewVMfbaT0aDU3D+3CfVf0pn3rSLdDM8b4gbBWQse2/ntDpiUZP6WqfLqhkP99bz2bCsu4qEciD1+bRZ+Utm6HZozxM/78GGafrrWIyGgR2SAi+SLyYAPvR4nIXOf9xSLSrc57DznlG0RkVFNtikiG08Ymp81Ip/wOEdkvIiudn8ln88H92fIdxdw4fRETX82luqaWF24bxKzJQy3BGGMalBIfzd5AHZMRkTDgWeAKoABYKiI5qrq2TrVJQLGq9hSRccDjwI0ikgWMA/oCnYCPRKSXs09jbT4OPKWqc0RkmtP2884+c1V16ll+Zr+VX1jGHxasZ0HePpLiovif75/PuMHpNu5ijDmllLbRfLq+EFX1u/vjfLlcNgTIV9UtACIyBxgD1E0yY4BfO9vzgT+L95OOAeaoaiWwVUTynfZoqE0RWQd8F7jZqTPTafdYkglK+0or+ONHG5mXW0B0eCvuu6IXky7OoLXd72KM8UFqfDTlVTWUVniIj4lwO5yT+PIt1hnYWed1ATC0sTqq6hGREiDRKV9Ub9/OznZDbSYCh1TV00B9gB+KyCXARuBnqlq3jYBTcrSaF/61mZe/2kpNrXL78K5MvbQniX54Q5Uxxn/VvSEzEJNMQ30v9bFOY+UNXf85VX2AfwBvqGqliPw73l7Od78VrMgUYApAly5dGmjOfRXVNfxl0Xb+/Gk+h8qr+f6ATtx3RW+6JMa6HZoxJgDVfQxz7xT/uinblyRTAKTXeZ0G7G6kToGIhAPxwMEm9m2o/ACQICLhTm/meH1VLapT/0W8YzffoqrTgekA2dnZ9ZOhq2pqlb+t2MWTH25k16GjXNIrmftH9fbLu3SNMYEjJd5/H8PsS5JZCmSKSAawC+9A/s316uQA44GFwFjgE1VVEckBZovIk3gH/jOBJXh7LN9q09nnU6eNOU6bfwcQkVRV3eMc7zpg3Rl+5nNOVflsw34e/2A96/cepl/neH4/tj8jeia5HZoxJgh0aBOFCH45w6zJJOOMsUwFFgBhwMuqmicijwC5qpoDzABedwb2D+JNGjj15uGdJOAB7lHVGoCG2nQO+QAwR0R+C6xw2ga4V0Suc9o5CNxx1p/+HFixo5jH3l/P4q0H6ZoYy59vvpCrzk/1yyW5jTGBKSKsFUlxUX7ZkxFVv7qi1Kyys7M1NzfXlWNv3l/GEws28P6avSTFRfKTyzK5cXAXIsNtOrIxpvld9+cvaRcbycyJQ5qu3AQRWaaq2c0Qlt3x39wKSyv408ebmLN0J9HhrfjZ5b2Y/B2bjmyMaVkpbaPZXlTudhjfYt98zeRIpYcXPt/Ci59vwVNby23DujL1uz398vkOxpjgkxofzaItRU1XPMcsyZwlT00tc3N38tSHmzhQVsnV/VO5f1Rvuia2djs0Y0wISYmPobTCw5FKj19dOfGfSALQws1FPPz3NWwqLGNwt3a8ePsgLuzSzu2wjDEh6Ni9MntLK+iRHOdyNCdYkjkDB49U8ei763hzeQHp7WOYdusgRvXt6HdrBhljQkdH5wmZ+0osyQS0r/MP8NO5Kykur+KeS3sw9dJMYiLD3A7LGBPiTtz171/TmC3J+EhVefGLLfzv++vpntSaVycMIauTLb1vjPEPKXUul/kTSzI+UFUee389L3y+hav7p/KHsf2JjbRTZ4zxH9ERYbSLjWCPnz2G2b4pffDSF1t54fMt3DasK7+5rq/drW+M8Usp8TF+d9e/3X7ehOU7ivnf99dxVb8USzDGGL+WGh/td2MylmROoaZWeWD+alLaRvP4D/tbgjHG+LWObaPZ52djMpZkTuGd1bvZVFjGL6/Jok20fz0IyBhj6kuNj+ZAWRWVnhq3QznOkswpzPhyK706xjGqb4rboRhjTJOOzTArLK10OZITLMk0Ysv+MlYXlPCj7HS7TGaMCQj+eK+MJZlGfLyuEICr+6e6HIkxxvim7mOY/YUlmUYs215Ml/axpDqPNTXGGH/nj49htiTTiFUFhxjYJcHtMIwxxmdxUeHERYX71V3/lmQaUFFdw56SCrr70SJzxhjji5T4aHYV2+Uyv1bg/A9Kb2+XyowxgaVnchz5hWVuh3GcT0lGREaLyAYRyReRBxt4P0pE5jrvLxaRbnXee8gp3yAio5pqU0QynDY2OW1GNnWM5nbwSBUAyXHRLXUIY4xpEb1T2rCt6AhHq/zjXpkmk4yIhAHPAlcCWcBNIpJVr9okoFhVewJPAY87+2YB44C+wGjgOREJa6LNx4GnVDUTKHbabvQYLaG8ygNgS/gbYwLOealtqFXYVHjY7VAA33oyQ4B8Vd2iqlXAHGBMvTpjgJnO9nzgMvE+wWsMMEdVK1V1K5DvtNdgm84+33XawGnz+00co9lVVHv/Aoi1JGOMCTD907wTlr7YdMDlSLx8WYW5M7CzzusCYGhjdVTVIyIlQKJTvqjevp2d7YbaTAQOqaqngfqNHeOkMykiU4ApzssyESmqX8dXWS3WV3JNEmd4LoKQnQsvOw8nBNW5mPo4TD2zXZOArs0Vhy9JpqHegvpYp7HyhnpQp6rvaxyo6nRg+vHARHJVNbuBfUOOnYsT7Fx42Xk4wc6Fl3MeujVXe75cLisA0uu8TgN2N1ZHRMKBeODgKfZtrPwAkOC0Uf9YjR3DGGOMn/IlySwFMp1ZX5F4B/Jz6tXJAcY722OBT1RVnfJxzsywDCATWNJYm84+nzpt4LT59yaOYYwxxk81ebnMGf+YCiwAwoCXVTVPRB4BclU1B5gBvC4i+Xh7F+OcffNEZB6wFvAA96hqDUBDbTqHfACYIyK/BVY4bdPYMXwwvekqIcPOxQl2LrzsPJxg58KrWc+DWGfAGGNMS7E7/o0xxrQYSzLGGGNaTFAnmaaWwwl0IvKyiBSKyJo6Ze1F5ENnWZ4PRaSdUy4i8rRzLlaLyMA6+4x36m8SkfENHcvfiUi6iHwqIutEJE9EfuKUh9T5EJFoEVkiIquc8/Abp/y0l2tqbEmoQOOsMrJCRN5xXofkuRCRbSLyjYisFJFcp6zlfz9UNSh/8E4o2Ax0ByKBVUCW23E182e8BBgIrKlT9nvgQWf7QeBxZ/sq4H289xsNAxY75e2BLc5/2znb7dz+bGdwLlKBgc52G2Aj3iWLQup8OJ8nztmOABY7n28eMM4pnwbc7Wz/GJjmbI8D5jrbWc7vTBSQ4fwuhbn9+c7wnNwHzAbecV6H5LkAtgFJ9cpa/PcjmHsyviyHE9BU9XO+fa9Q3eV36i/L85p6LcJ7P1IqMAr4UFUPqmox8CHedeYCiqruUdXlzvZhYB3eVSJC6nw4n+fYErwRzo9y+ss1NbYkVEARkTTgauAl5/WZLF0VFOeiES3++xHMSaah5XA6N1I3mHRU1T3g/eIFOjjljZ2PoDtPzmWOC/H+FR9y58O5PLQSKMT7JbAZH5drAuouCRXQ58HxR+B+oNZ57fPSVQTfuVDgnyKyTLzLb8E5+P3wZVmZQOXTMjQh5HSX/glIIhIHvAn8VFVLpfE1VIP2fKj3XrQBIpIAvA2c11A1579Bex5E5BqgUFWXicjIY8UNVA36c+EYoaq7RaQD8KGIrD9F3WY7F8Hck/FlOZxgtM/p1uL8t9ApP90lfgKOiETgTTCzVPUtpzhkz4eqHgI+w3tN/XSXawqG8zACuE5EtuG9XP5dvD2bUDwXqOpu57+FeP/4GMI5+P0I5iTjy3I4waju8jv1l+W53Zk1MgwocbrHC4DviUg7Z2bJ95yygOJcO58BrFPVJ+u8FVLnQ0SSnR4MIhIDXI53fOp0l2tqbEmogKGqD6lqmnoXexyH97PdQgieCxFpLSJtjm3j/Xe9hnPx++H2jIeW/ME7Q2Ij3mvS/+V2PC3w+d4A9gDVeP/CmIT3GvLHwCbnv+2duoL3QXGbgW+A7DrtTMQ7mJkPTHD7c53hubgYb7d9NbDS+bkq1M4H0B/vckyrnS+Rh53y7ni/GPOBvwJRTnm08zrfeb97nbb+yzk/G4Ar3f5sZ3leRnJidlnInQvnM69yfvKOfR+ei98PW1bGGGNMiwnmy2XGGGNcZknGGGNMi7EkY4wxpsVYkjHGGNNiLMkYY4xpMZZkjDHGtBhLMsYYY1rM/we5gZknDF8RXwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] @@ -554,31 +550,37 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 10, "metadata": {}, "outputs": [], "source": [ - "# total_f.update_integration_options(draws_per_dim=10000000, mc_sampler=None)\n", - "# inte = total_f.integrate(limits = (3000, 3200), norm_range=False)\n", - "# print(zfit.run(inte))\n", - "# print(pdg[\"jpsi_BR\"]/pdg[\"NR_BR\"], zfit.run(inte)/pdg[\"NR_auc\"])" + "# total_f.update_integration_options(draws_per_dim=2000000, mc_sampler=None)\n", + "# inte = total_f.integrate(limits = (3090, 3102), norm_range=False)\n", + "# inte_fl = zfit.run(inte)\n", + "# print(inte_fl)\n", + "# print(pdg[\"jpsi_BR\"]/pdg[\"NR_BR\"], inte_fl/pdg[\"NR_auc\"])" ] }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 11, "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": 13, + "execution_count": 12, "metadata": {}, "outputs": [], "source": [ @@ -610,7 +612,7 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 13, "metadata": {}, "outputs": [ { @@ -627,7 +629,7 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 14, "metadata": {}, "outputs": [], "source": [ @@ -648,58 +650,53 @@ }, { "cell_type": "code", - "execution_count": 16, - "metadata": {}, - "outputs": [], - "source": [ - "nevents = int(pdg[\"number_of_decays\"])\n", - "event_stack = 100000\n", - "\n", - "calls = int(nevents/event_stack + 1)\n", - "\n", - "# total_samp = []\n", - "\n", - "# start = time.time()\n", - "\n", - "# for call in range(calls):\n", - "# samp = total_f.sample(n=event_stack)\n", - "# sam = samp.unstack_x()\n", - "# sam = zfit.run(sam)\n", - "# clear_output(wait=True)\n", - " \n", - "# print(\"{0}/{1}\".format(call + 1, calls))\n", - "# print(\"Time taken: {}\".format(display_time(int(time.time() - start))))\n", - "# c = call + 1\n", - "# print(\"Projected time left: {}\".format(display_time(int((time.time() - start)/c*(calls-c)))))\n", - " \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": 17, - "metadata": {}, - "outputs": [], - "source": [ - "# print(\"Time to generate full toy: {} s\".format(int(time.time()-start)))" - ] - }, - { - "cell_type": "code", - "execution_count": 18, + "execution_count": null, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "(5404696,)\n" + "1/1081\n", + "Time taken: 5 min, 13 s\n", + "Projected time left: 3 d, 21 h\n" ] } ], "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", + "\n", + "for call in range(calls):\n", + " sam = samp.unstack_x()\n", + " sam = zfit.run(sam)\n", + " clear_output(wait=True)\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", @@ -716,22 +713,9 @@ }, { "cell_type": "code", - "execution_count": 19, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAD8CAYAAABDwhLXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAFslJREFUeJzt3X+MXeWd3/H3p3YgJFliAyalNls7jdUuQduGjAhtqlUUtmAgivmDSEZRsbJIVlPSZptWiWmkoiZdibTVskVKWKHgxlRpCGWzwgoQr0WIokqBMOQHP0KIJ4TCLGzsxMCyjZos2W//uM+Ey+TODJ4Znrn2vF/S1T3ne55znnMe+frjc+7DJVWFJEmvtr+10icgSVodDBxJUhcGjiSpCwNHktSFgSNJ6sLAkSR1YeBIkrowcCRJXRg4kqQu1q70CfRy2mmn1ebNm1f6NCTpmPLAAw/8pKo2LMexVk3gbN68mcnJyZU+DUk6piT5P8t1LB+pSZK6MHAkSV0YOJKkLgwcSVIXBo4kqQsDR5LUhYEjSepiwcBJsifJoSQPj9j275JUktPaepJcn2QqyYNJzhlquzPJwfbaOVR/e5KH2j7XJ0mrn5LkQGt/IMn6hfqQJI2vV3KH8zlg2+xikjOBfwY8OVS+CNjaXruAG1rbU4BrgHcA5wLXzARIa7NraL+ZvnYDd1fVVuDutj5nH5Kk8bZg4FTV14EjIzZdB3wUqKHaduDmGrgXWJfkDOBC4EBVHamqZ4EDwLa27eSq+kZVFXAzcOnQsfa25b2z6qP6kPQq2Lz7jpU+BR0nFvUdTpL3An9eVd+dtWkj8NTQ+nSrzVefHlEHeFNVPQPQ3k9foA9J0hg76t9SS/I64OPABaM2j6jVIurznsIr3SfJLgaP3fjN3/zNBQ4rSXo1LeYO5+8BW4DvJnkC2AR8K8nfZnC3ceZQ203A0wvUN42oA/x45lFZez/U6nMd69dU1Y1VNVFVExs2LMuPnUqSFumoA6eqHqqq06tqc1VtZhAA51TVXwD7gCvaTLLzgOfb47D9wAVJ1rfJAhcA+9u2F5Kc12anXQHc3rraB8zMZts5qz6qD0nSGFvwkVqSLwDvAk5LMg1cU1U3zdH8TuBiYAr4GfABgKo6kuSTwP2t3SeqamYiwgcZzIQ7CbirvQCuBW5NciWDmXDvm68PSdJ4WzBwquryBbZvHlou4Ko52u0B9oyoTwJnj6j/FDh/RH3OPiRJ48tfGpAkdWHgSJK6MHAkSV0YOJKkLgwcSVIXBo4kqQsDR5LUhYEjSerCwJEkdWHgSJK6MHAkSV0YOJKkLgwcSVIXBo4kqQsDR5LUhYEjSerCwJEkdWHgSJK6MHAkSV0YOJKkLhYMnCR7khxK8vBQ7b8k+X6SB5P8aZJ1Q9uuTjKV5LEkFw7Vt7XaVJLdQ/UtSe5LcjDJF5Oc0OontvWptn3zQn1IksbXK7nD+RywbVbtAHB2Vf028APgaoAkZwE7gLe2fT6TZE2SNcCngYuAs4DLW1uATwHXVdVW4Fngyla/Eni2qt4CXNfazdnHUV63JKmzBQOnqr4OHJlV+7OqerGt3gtsasvbgVuq6udV9SNgCji3vaaq6vGq+gVwC7A9SYB3A7e1/fcClw4da29bvg04v7Wfqw9J0hhbju9wfg+4qy1vBJ4a2jbdanPVTwWeGwqvmfrLjtW2P9/az3UsSdIYW1LgJPk48CLw+ZnSiGa1iPpijjXq/HYlmUwyefjw4VFNJEmdLDpwkuwE3gO8v6pm/sKfBs4carYJeHqe+k+AdUnWzqq/7Fht+xsZPNqb61i/pqpurKqJqprYsGHDYi5TkrRMFhU4SbYBHwPeW1U/G9q0D9jRZphtAbYC3wTuB7a2GWknMPjSf18LqnuAy9r+O4Hbh461sy1fBny1tZ+rD0nSGFu7UIMkXwDeBZyWZBq4hsGstBOBA4Pv8bm3qv5FVT2S5FbgewwetV1VVb9sx/kQsB9YA+ypqkdaFx8Dbknyn4BvAze1+k3A/0gyxeDOZgfAfH1IksZXXnoadnybmJioycnJlT4N6ZizefcdPHHtJSt9GlohSR6oqonlOJa/NCBJ6sLAkSR1YeBIkrowcCRJXRg4kqQuDBxJUhcGjiSpCwNHktSFgSNJ6sLAkSR1YeBIkrowcCRJXRg4kqQuDBxJUhcGjiSpCwNHktSFgSNJ6sLAkSR1YeBIkrowcCRJXRg4kqQuFgycJHuSHEry8FDtlCQHkhxs7+tbPUmuTzKV5MEk5wzts7O1P5hk51D97UkeavtcnySL7UOSNL5eyR3O54Bts2q7gburaitwd1sHuAjY2l67gBtgEB7ANcA7gHOBa2YCpLXZNbTftsX0IUkabwsGTlV9HTgyq7wd2NuW9wKXDtVvroF7gXVJzgAuBA5U1ZGqehY4AGxr206uqm9UVQE3zzrW0fQhSRpji/0O501V9QxAez+91TcCTw21m261+erTI+qL6ePXJNmVZDLJ5OHDh4/qAiVJy2u5Jw1kRK0WUV9MH79erLqxqiaqamLDhg0LHFaS9GpabOD8eOYxVns/1OrTwJlD7TYBTy9Q3zSivpg+JEljbLGBsw+YmWm2E7h9qH5Fm0l2HvB8exy2H7ggyfo2WeACYH/b9kKS89rstCtmHeto+pAkjbG1CzVI8gXgXcBpSaYZzDa7Frg1yZXAk8D7WvM7gYuBKeBnwAcAqupIkk8C97d2n6iqmYkIH2QwE+4k4K724mj7kCSNtwUDp6oun2PT+SPaFnDVHMfZA+wZUZ8Ezh5R/+nR9iFJGl/+0oAkqQsDR5LUhYEjSerCwJEkdWHgSJK6MHAkSV0YOJKkLgwcSVIXBo4kqQsDR5LUhYEjSerCwJEkdWHgSJK6MHAkSV0YOJKkLgwcSVIXBo4kqQsDR5LUhYEjSerCwJEkdbGkwEnyb5I8kuThJF9I8tokW5Lcl+Rgki8mOaG1PbGtT7Xtm4eOc3WrP5bkwqH6tlabSrJ7qD6yD0nS+Fp04CTZCPxrYKKqzgbWADuATwHXVdVW4FngyrbLlcCzVfUW4LrWjiRntf3eCmwDPpNkTZI1wKeBi4CzgMtbW+bpQ5I0ppb6SG0tcFKStcDrgGeAdwO3te17gUvb8va2Ttt+fpK0+i1V9fOq+hEwBZzbXlNV9XhV/QK4Bdje9pmrD0nSmFp04FTVnwP/FXiSQdA8DzwAPFdVL7Zm08DGtrwReKrt+2Jrf+pwfdY+c9VPnacPSdKYWsojtfUM7k62AH8HeD2Dx1+z1cwuc2xbrvqoc9yVZDLJ5OHDh0c1kSR1spRHar8L/KiqDlfVXwNfAv4JsK49YgPYBDzdlqeBMwHa9jcCR4brs/aZq/6Tefp4maq6saomqmpiw4YNS7hUSdJSLSVwngTOS/K69r3K+cD3gHuAy1qbncDtbXlfW6dt/2pVVavvaLPYtgBbgW8C9wNb24y0ExhMLNjX9pmrD0nSmFrKdzj3Mfji/lvAQ+1YNwIfAz6SZIrB9y03tV1uAk5t9Y8Au9txHgFuZRBWXwGuqqpftu9oPgTsBx4Fbm1tmacPSdKYyuCG4fg3MTFRk5OTK30a0jFn8+47eOLaS1b6NLRCkjxQVRPLcSx/aUCS1IWBI0nqwsCRJHVh4EiSujBwJEldGDiSpC4MHElSFwaOJKkLA0eS1IWBI0nqwsCRJHVh4Eg6Zm3efcdKn4KOgoEjSerCwJEkdWHgSJK6MHAkSV0YOJKkLgwcSVIXBo4kqQsDR5LUxZICJ8m6JLcl+X6SR5P84ySnJDmQ5GB7X9/aJsn1SaaSPJjknKHj7GztDybZOVR/e5KH2j7XJ0mrj+xDkjS+lnqH89+Ar1TVPwD+IfAosBu4u6q2Ane3dYCLgK3ttQu4AQbhAVwDvAM4F7hmKEBuaG1n9tvW6nP1IUkaU4sOnCQnA78D3ARQVb+oqueA7cDe1mwvcGlb3g7cXAP3AuuSnAFcCByoqiNV9SxwANjWtp1cVd+oqgJunnWsUX1IksbUUu5w3gwcBv57km8n+WyS1wNvqqpnANr76a39RuCpof2nW22++vSIOvP0IUkaU0sJnLXAOcANVfU24P8y/6OtjKjVIuqvWJJdSSaTTB4+fPhodpUkLbOlBM40MF1V97X12xgE0I/b4zDa+6Gh9mcO7b8JeHqB+qYRdebp42Wq6saqmqiqiQ0bNizqIiVJy2PRgVNVfwE8leTvt9L5wPeAfcDMTLOdwO1teR9wRZutdh7wfHscth+4IMn6NlngAmB/2/ZCkvPa7LQrZh1rVB+SpDG1don7/yvg80lOAB4HPsAgxG5NciXwJPC+1vZO4GJgCvhZa0tVHUnySeD+1u4TVXWkLX8Q+BxwEnBXewFcO0cfkqQxtaTAqarvABMjNp0/om0BV81xnD3AnhH1SeDsEfWfjupDkjS+/KUBSVIXBo4kqQsDR5LUhYEjSerCwJEkdWHgSJK6MHAkSV0YOJKkLgwcSVIXBo4kqQsDR5LUhYEjSerCwJEkdWHgSJK6MHAkSV0YOJKkLgwcSVIXBo4kqQsDR5LUhYEjSerCwJEkdbHkwEmyJsm3k3y5rW9Jcl+Sg0m+mOSEVj+xrU+17ZuHjnF1qz+W5MKh+rZWm0qye6g+sg9J0vhajjucDwOPDq1/CriuqrYCzwJXtvqVwLNV9RbgutaOJGcBO4C3AtuAz7QQWwN8GrgIOAu4vLWdrw9J0phaUuAk2QRcAny2rQd4N3Bba7IXuLQtb2/rtO3nt/bbgVuq6udV9SNgCji3vaaq6vGq+gVwC7B9gT4kSWNqqXc4fwR8FPibtn4q8FxVvdjWp4GNbXkj8BRA2/58a/+r+qx95qrP18fLJNmVZDLJ5OHDhxd7jZKkZbDowEnyHuBQVT0wXB7RtBbYtlz1Xy9W3VhVE1U1sWHDhlFNJEmdrF3Cvu8E3pvkYuC1wMkM7njWJVnb7kA2AU+39tPAmcB0krXAG4EjQ/UZw/uMqv9knj4kSWNq0Xc4VXV1VW2qqs0MvvT/alW9H7gHuKw12wnc3pb3tXXa9q9WVbX6jjaLbQuwFfgmcD+wtc1IO6H1sa/tM1cfkqQx9Wr8dzgfAz6SZIrB9y03tfpNwKmt/hFgN0BVPQLcCnwP+ApwVVX9st29fAjYz2AW3K2t7Xx9SJLG1FIeqf1KVX0N+FpbfpzBDLPZbf4f8L459v8D4A9G1O8E7hxRH9mHJGl8+UsDkqQuDBxJUhcGjiSpCwNHktSFgSNJ6sLAkSR1YeBIkrowcCRJXRg4kqQuDBxJUhcGjiSpCwNHktSFgSNJ6sLAkSR1YeBIkrowcCRJXRg4kqQuDBxJUhcGjiSpCwNHktTFogMnyZlJ7knyaJJHkny41U9JciDJwfa+vtWT5PokU0keTHLO0LF2tvYHk+wcqr89yUNtn+uTZL4+JEnjayl3OC8C/7aqfgs4D7gqyVnAbuDuqtoK3N3WAS4CtrbXLuAGGIQHcA3wDuBc4JqhALmhtZ3Zb1urz9WHJGlMLTpwquqZqvpWW34BeBTYCGwH9rZme4FL2/J24OYauBdYl+QM4ELgQFUdqapngQPAtrbt5Kr6RlUVcPOsY43qQ5I0ppblO5wkm4G3AfcBb6qqZ2AQSsDprdlG4Kmh3aZbbb769Ig68/Qx+7x2JZlMMnn48OHFXp4kaRksOXCSvAH4E+D3q+ov52s6olaLqL9iVXVjVU1U1cSGDRuOZldJ0jJbUuAkeQ2DsPl8VX2plX/cHofR3g+1+jRw5tDum4CnF6hvGlGfrw9J0phayiy1ADcBj1bVHw5t2gfMzDTbCdw+VL+izVY7D3i+PQ7bD1yQZH2bLHABsL9teyHJea2vK2Yda1QfkqQxtXYJ+74T+OfAQ0m+02r/HrgWuDXJlcCTwPvatjuBi4Ep4GfABwCq6kiSTwL3t3afqKojbfmDwOeAk4C72ot5+pAkjalFB05V/W9Gf88CcP6I9gVcNcex9gB7RtQngbNH1H86qg9J0vjylwYkSV0YOJKkLgwcSVIXBo4kqQsDR5LUhYEjSerCwJEkdWHgSJK6MHAkSV0YOJKkLgwcSVIXBo4kqQsDR5LUhYEjSerCwJEkdWHgSJK6MHAkSV0YOJKkLgwcSVIXBo4kqYtjOnCSbEvyWJKpJLtX+nwkSXM7ZgMnyRrg08BFwFnA5UnOWtmzkiTNZe1Kn8ASnAtMVdXjAEluAbYD31vRs5K0JJt338ET117C5t13vOL24+6Jay9Z6VMYC8dy4GwEnhpanwbesULncsw6Fj6sWnm9/5wcb38ue17POIfbsRw4GVGrlzVIdgG72upfJfkp8JNX+8SOEafhWMxwLAYch5ccs2ORTy3r4U4D/u5yHexYDpxp4Myh9U3A08MNqupG4MaZ9SSTVTXR5/TGm2PxEsdiwHF4iWMx0MZh83Id75idNADcD2xNsiXJCcAOYN8Kn5MkaQ7H7B1OVb2Y5EPAfmANsKeqHlnh05IkzeGYDRyAqroTuPModrlx4SarhmPxEsdiwHF4iWMxsKzjkKpauJUkSUt0LH+HI0k6hqyawDnefwYnyZ4kh5I8PFQ7JcmBJAfb+/pWT5Lr21g8mOScoX12tvYHk+xciWtZqiRnJrknyaNJHkny4VZfVeOR5LVJvpnku20c/mOrb0lyX7umL7ZJNyQ5sa1Pte2bh451das/luTClbmipUuyJsm3k3y5ra/KsUjyRJKHknwnyWSrvfqfj6o67l8MJhX8EHgzcALwXeCslT6vZb7G3wHOAR4eqv1nYHdb3g18qi1fDNzF4L9lOg+4r9VPAR5v7+vb8vqVvrZFjMUZwDlt+TeAHzD4+aNVNR7tet7Qll8D3Neu71ZgR6v/MfDBtvwvgT9uyzuAL7bls9pn5kRgS/ssrVnp61vkmHwE+J/Al9v6qhwL4AngtFm1V/3zsVrucH71MzhV9Qtg5mdwjhtV9XXgyKzydmBvW94LXDpUv7kG7gXWJTkDuBA4UFVHqupZ4ACw7dU/++VVVc9U1bfa8gvAowx+mWJVjUe7nr9qq69prwLeDdzW6rPHYWZ8bgPOT5JWv6Wqfl5VPwKmGHymjilJNgGXAJ9t62GVjsUcXvXPx2oJnFE/g7Nxhc6lpzdV1TMw+EsYOL3V5xqP426c2qOQtzH41/2qG4/2COk7wCEGfyH8EHiuql5sTYav6VfX27Y/D5zKcTAOzR8BHwX+pq2fyuodiwL+LMkDGfwiC3T4fBzT06KPwoI/g7PKzDUex9U4JXkD8CfA71fVXw7+gTq66YjacTEeVfVL4B8lWQf8KfBbo5q19+N2HJK8BzhUVQ8keddMeUTT434smndW1dNJTgcOJPn+PG2XbSxWyx3Ogj+Dc5z6cbv1pb0favW5xuO4Gackr2EQNp+vqi+18qodj6p6Dvgag2fw65LM/GNz+Jp+db1t+xsZPKY9HsbhncB7kzzB4JH6uxnc8azGsaCqnm7vhxj8Q+RcOnw+VkvgrNafwdkHzMwc2QncPlS/os0+OQ94vt1C7wcuSLK+zVC5oNWOKe1Z+03Ao1X1h0ObVtV4JNnQ7mxIchLwuwy+z7oHuKw1mz0OM+NzGfDVGnw7vA/Y0WZubQG2At/scxXLo6qurqpNNfhdsB0Mru39rMKxSPL6JL8xs8zgz/XD9Ph8rPRsiV4vBjMtfsDgGfbHV/p8XoXr+wLwDPDXDP7lcSWDZ853Awfb+ymtbRj8z+t+CDwETAwd5/cYfBE6BXxgpa9rkWPxTxnc2j8IfKe9Ll5t4wH8NvDtNg4PA/+h1d/M4C/JKeB/ASe2+mvb+lTb/uahY328jc9jwEUrfW1LHJd38dIstVU3Fu2av9tej8z8fdjj8+EvDUiSulgtj9QkSSvMwJEkdWHgSJK6MHAkSV0YOJKkLgwcSVIXBo4kqQsDR5LUxf8Htmuht2g0nfcAAAAASUVORK5CYII=\n", - "text/plain": [ - "
" - ] - }, - "metadata": { - "needs_background": "light" - }, - "output_type": "display_data" - } - ], + "outputs": [], "source": [ "bins = int((x_max-x_min)/7)\n", "\n", @@ -741,10 +725,10 @@ "\n", "# plt.plot(sam, calcs, '.')\n", "# plt.plot(test_q, calcs_test)\n", - "# plt.ylim(6000, 10000)\n", - "# plt.xlim(3000, 3750)\n", + "plt.ylim(4000, 12000)\n", + "plt.xlim(3000, 3750)\n", "\n", - "plt.savefig('test.png')" + "plt.savefig('test2.png')" ] }, { @@ -756,7 +740,7 @@ }, { "cell_type": "code", - "execution_count": 20, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -786,316 +770,9 @@ }, { "cell_type": "code", - "execution_count": 21, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": { - "text/html": [ - "
" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/html": [ - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
FCN = 19224270.58100143TOTAL NCALL = 58NCALLS = 58
EDM = 1.5229666301917304e-07GOAL EDM = 5e-06\n", - " UP = 0.5
\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
ValidValid ParamAccurate CovarPosDefMade PosDef
TrueTrueTrueTrueFalse
Hesse FailHasCovAbove EDMReach calllim
FalseTrueFalseFalse
" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/html": [ - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
+NameValueHesse ErrorMinos Error-Minos Error+Limit-Limit+Fixed?
0psi2s_p-9.424780.00651368No
1jpsi_s444.4950.220516No
2jpsi_p-9.424780.00631191No
3psi2s_s74.73490.048852No
\n", - "
\n",
-       "\n",
-       "
" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/html": [ - "
" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/html": [ - "Minos status for psi2s_p: VALID\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
Error-0.0065151409160649580.006517797812598556
ValidTrueTrue
At LimitFalseFalse
Max FCNFalseFalse
New MinFalseFalse
" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/html": [ - "Minos status for jpsi_s: VALID\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
Error-0.220659286900404150.22053488031823795
ValidTrueTrue
At LimitFalseFalse
Max FCNFalseFalse
New MinFalseFalse
" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/html": [ - "Minos status for jpsi_p: VALID\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
Error-0.0063110017718012870.006310920098439091
ValidTrueTrue
At LimitFalseFalse
Max FCNFalseFalse
New MinFalseFalse
" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/html": [ - "Minos status for psi2s_s: VALID\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
Error-0.049227653394393410.04926101726478015
ValidTrueTrue
At LimitFalseFalse
Max FCNFalseFalse
New MinFalseFalse
" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "psi2s_p: ^{+0.006517797812598556}_{-0.006515140916064958}\n", - "jpsi_s: ^{+0.22053488031823795}_{-0.22065928690040415}\n", - "jpsi_p: ^{+0.006310920098439091}_{-0.006311001771801287}\n", - "psi2s_s: ^{+0.04926101726478015}_{-0.04922765339439341}\n", - "Function minimum: 19224270.58100143\n" - ] - } - ], + "outputs": [], "source": [ "nll = zfit.loss.UnbinnedNLL(model=total_f, data=data2, fit_range = (x_min, x_max))\n", "\n", @@ -1113,22 +790,11 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "-3.1420346928204133" - ] - }, - "execution_count": 22, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ - "(-9.42522+2*np.pi)#/np.pi" + "(-3.14+2*np.pi)/np.pi" ] }, { @@ -1140,37 +806,18 @@ }, { "cell_type": "code", - "execution_count": 23, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "'5 h, 55 min'" - ] - }, - "execution_count": 23, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "display_time(int(395*pdg[\"number_of_decays\"]/100000))" ] }, { "cell_type": "code", - "execution_count": 24, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "6 h, 12 min\n" - ] - } - ], + "outputs": [], "source": [ "print(display_time(22376))" ] @@ -1180,6 +827,36 @@ "execution_count": null, "metadata": {}, "outputs": [], + "source": [ + "probs = total_f.pdf(test_q)\n", + "\n", + "calcs_test = zfit.run(probs)\n", + "res_y = zfit.run(jpsi_res(test_q))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.clf()\n", + "# plt.plot(x_part, calcs, '.')\n", + "plt.plot(test_q, calcs_test, label = 'pdf')\n", + "# plt.plot(test_q, res_y, label = 'res')\n", + "plt.legend()\n", + "plt.ylim(0.0, 4e-4)\n", + "# plt.yscale('log')\n", + "# plt.xlim(3080, 3110)\n", + "plt.savefig('test3.png')\n", + "print(jpsi_width)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], "source": [] } ], diff --git a/__pycache__/pdg_const.cpython-37.pyc b/__pycache__/pdg_const.cpython-37.pyc index c499ac2..af58ea9 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_1/0.pkl b/data/zfit_toys/toy_1/0.pkl index e3599d9..e3b4e36 100644 --- a/data/zfit_toys/toy_1/0.pkl +++ b/data/zfit_toys/toy_1/0.pkl Binary files differ diff --git a/pdg_const.py b/pdg_const.py index aedf951..f831ad0 100644 --- a/pdg_const.py +++ b/pdg_const.py @@ -66,15 +66,15 @@ "NR_auc": 0.00133, #Resonances format(mass, width, phase, scale) -"jpsi": (3096.0, 0.09, -1.5, 2e-2), +# "jpsi": (3096.0, 0.09, -1.5, 2e-2), #---> prescaling +"jpsi": (3096.0, 0.09, -1.5, 184.39), #---> after scaling # "jpsi": (3096.0, 0.09, -1.5, 0.0), "jpsi_BR": 6.02e-5, -"jpsi_auc": 2.336e-9, -"factor_jpsi": 8856.2, +"jpsi_auc": 0.2126825758464027, -"psi2s": (3686.0, 0.3, -1.5, 3.14e-3), +# "psi2s": (3686.0, 0.3, -1.5, 3.14e-3), #---> prescaling +"psi2s": (3686.0, 0.3, -1.5, 23.07), #---> after scaling # "psi2s": (3686.0, 0.3, -1.5, 0.0), "psi2s_BR": 4.97e-6, -"psi2s_auc": 2.853e-10, -"factor_psi2s": 7281.4 +"psi2s_auc": 2.802257483178487e-10, } diff --git a/raremodel-nb.ipynb b/raremodel-nb.ipynb index 3980d05..4047084 100644 --- a/raremodel-nb.ipynb +++ b/raremodel-nb.ipynb @@ -16,7 +16,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\zfit\\util\\execution.py:53: UserWarning: Not running on Linux. Determining available cpus for thread can failand be overestimated. Workaround (only if too many cpus are used):`zfit.run.set_n_cpu(your_cpu_number)`\n", + "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" ] }, @@ -151,16 +151,13 @@ "\n", " #Rotate by the phase\n", "\n", - " r = tf.abs(com)\n", + " r = ztf.to_complex(scale*tf.abs(com))\n", "\n", " _phase = tf.angle(com)\n", "\n", " _phase += phase\n", "\n", - " x = tf.cos(phase)*r\n", - " y = tf.sin(phase)*r\n", - "\n", - " com = tf.complex(scale* x, scale * y)\n", + " com = r * tf.exp(tf.complex(ztf.constant(0.0), _phase))\n", "\n", " return com\n", "\n", @@ -269,7 +266,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ @@ -324,15 +321,15 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "class total_pdf(zfit.pdf.ZPDF):\n", " _N_OBS = 1 # dimension, can be omitted\n", " _PARAMS = ['jpsi_mass', 'jpsi_scale', 'jpsi_phase', 'jpsi_width',\n", - " 'psi2s_mass', 'psi2s_scale', 'psi2s_phase', 'psi2s_width',\n", - " 'cusp_mass', 'sigma_L', 'sigma_R', 'cusp_scale'\n", + " 'psi2s_mass', 'psi2s_scale', 'psi2s_phase', 'psi2s_width'#,\n", + " #'cusp_mass', 'sigma_L', 'sigma_R', 'cusp_scale'\n", " ] # the name of the parameters\n", "\n", " def _unnormalized_pdf(self, x):\n", @@ -348,7 +345,7 @@ " def cusp(q):\n", " return bifur_gauss(q, mean = self.params['cusp_mass'], sigma_L = self.params['sigma_L'], sigma_R = self.params['sigma_R'], scale = self.params['cusp_scale'])\n", "\n", - " funcs = jpsi_res(x) + psi2s_res(x) + cusp(x)\n", + " funcs = jpsi_res(x) + psi2s_res(x) #+ cusp(x)\n", "\n", " vec_f = vec(x, funcs)\n", "\n", @@ -368,7 +365,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 5, "metadata": {}, "outputs": [], "source": [ @@ -396,7 +393,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 6, "metadata": {}, "outputs": [ { @@ -413,7 +410,7 @@ "#jpsi\n", "\n", "jpsi_mass, jpsi_width, jpsi_phase, jpsi_scale = pdg[\"jpsi\"]\n", - "jpsi_scale *= pdg[\"factor_jpsi\"]\n", + "# jpsi_scale *= pdg[\"factor_jpsi\"]\n", "\n", "jpsi_m = zfit.Parameter(\"jpsi_m\", ztf.constant(jpsi_mass), floating = False)\n", "jpsi_w = zfit.Parameter(\"jpsi_w\", ztf.constant(jpsi_width), floating = False)\n", @@ -423,7 +420,6 @@ "#psi2s\n", "\n", "psi2s_mass, psi2s_width, psi2s_phase, psi2s_scale = pdg[\"psi2s\"]\n", - "psi2s_scale *= pdg[\"factor_psi2s\"]\n", "\n", "psi2s_m = zfit.Parameter(\"psi2s_m\", ztf.constant(psi2s_mass), floating = False)\n", "psi2s_w = zfit.Parameter(\"psi2s_w\", ztf.constant(psi2s_width), floating = False)\n", @@ -432,12 +428,12 @@ "\n", "#cusp\n", "\n", - "cusp_mass, sigma_R, sigma_L, cusp_scale = 3550, 3e-7, 200, 0\n", + "# cusp_mass, sigma_R, sigma_L, cusp_scale = 3550, 3e-7, 200, 0\n", "\n", - "cusp_m = zfit.Parameter(\"cusp_m\", ztf.constant(cusp_mass), floating = False)\n", - "sig_L = zfit.Parameter(\"sig_L\", ztf.constant(sigma_L), floating = False)\n", - "sig_R = zfit.Parameter(\"sig_R\", ztf.constant(sigma_R), floating = False)\n", - "cusp_s = zfit.Parameter(\"cusp_s\", ztf.constant(cusp_scale), floating = False)" + "# cusp_m = zfit.Parameter(\"cusp_m\", ztf.constant(cusp_mass), floating = False)\n", + "# sig_L = zfit.Parameter(\"sig_L\", ztf.constant(sigma_L), floating = False)\n", + "# sig_R = zfit.Parameter(\"sig_R\", ztf.constant(sigma_R), floating = False)\n", + "# cusp_s = zfit.Parameter(\"cusp_s\", ztf.constant(cusp_scale), floating = False)" ] }, { @@ -449,13 +445,13 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "total_f = total_pdf(obs=obs, jpsi_mass = jpsi_m, jpsi_scale = jpsi_s, jpsi_phase = jpsi_p, jpsi_width = jpsi_w,\n", - " psi2s_mass = psi2s_m, psi2s_scale = psi2s_s, psi2s_phase = psi2s_p, psi2s_width = psi2s_w,\n", - " cusp_mass = cusp_m, sigma_L = sig_L, sigma_R = sig_R, cusp_scale = cusp_s)\n", + " psi2s_mass = psi2s_m, psi2s_scale = psi2s_s, psi2s_phase = psi2s_p, psi2s_width = psi2s_w)#,\n", + " #cusp_mass = cusp_m, sigma_L = sig_L, sigma_R = sig_R, cusp_scale = cusp_s)\n", "\n", "# print(total_pdf.obs)" ] @@ -469,7 +465,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 8, "metadata": {}, "outputs": [], "source": [ @@ -509,7 +505,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 9, "metadata": {}, "outputs": [ { @@ -521,7 +517,7 @@ }, { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAD8CAYAAACl69mTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8FdX5+PHPk4WEkBAghLAkkABhCbIIEXEtdSm44oISrYoLoha/Xfy1Vtpv1Vrtt7YWq9YNBcUVKHWJFhcUtypbkEXCGvZAgEBIAtmX5/fHHTDEhFwgYe7yvF+vvJg798yZZ4bkPvfMmXNGVBVjjDGmJYS4HYAxxpjAZUnGGGNMi7EkY4wxpsVYkjHGGNNiLMkYY4xpMZZkjDHGtBivkoyIjBaRdSKSIyL3NfB+hIjMct5fJCLJdd6b7KxfJyKjjqHOp0TkoDf7MMYY45uaTDIiEgo8DVwEpAHXiUhavWK3AftVtTfwOPCos20akAEMAEYDz4hIaFN1ikg60M6bfRhjjPFd3rRkhgM5qrpJVSuBmcCYemXGADOc5TnA+SIizvqZqlqhqpuBHKe+Rut0EtDfgHu93IcxxhgfFeZFmW7A9jqvc4HTGyujqtUiUgTEOesX1tu2m7PcWJ13A5mqmlcvhzS2j711C4nIRGAiQJs2bYb169fPi0M0xhyy92AFeUXlxLVpRdd2rd0O5wjVtcqavGIABnaLdTmawLV06dK9qhrfHHV5k2Qaai3Un4umsTKNrW+oBaUi0hW4Bhh5nHGgqlOBqQDp6emalZXVwGbGmMZM++9m/vT+am4+M5kHLx/gdjhH2HuwgvSHPwEg6y+XuBxN4BKRrc1VlzeXy3KBpDqvE4GdjZURkTAgFig4yraNrT8V6A3kiMgWIEpEcprYhzGmGfnyNeham2vR73iTZJYAqSKSIiKt8HTkZ9YrkwmMd5bHAvPVM/NmJpDh3BmWAqQCixurU1X/o6qdVTVZVZOBUqej/2j7MMY0I5/u6bS/eL/T5OUyp//jbuAjIBSYrqrZIvIQkKWqmcA04FWn1VGAJ2nglJsNrAaqgUmqWgPQUJ1NhNLgPowxLcMXv8PV+l5Ipgne9MmgqnOBufXW3V9nuRxPX0pD2z4CPOJNnQ2UifZmH8eiqqqK3NxcysvLT7Qq10VGRpKYmEh4eLjboZgAcqgh44uf5+qTUZmj8SrJBJLc3FxiYmJITk7Gn++AVlX27dtHbm4uKSkpbodjzElhLRn/E3TTypSXlxMXF+fXCQZARIiLiwuIFpkx3qq1LON3gi7JAH6fYA4JlOMwvuXQ75UPdskYPxSUScYY0zhf/u5itzD7H0syPurzzz/n0ksvBaCiooILLriAIUOGMGvWLJcjM8HCFzvZLcf4n6Dr+PdHy5Yto6qqiuXLl7sdigkCh+8u88EPdGvJ+B9rybhgy5Yt9OvXj/HjxzNo0CDGjh1LaWkpH374If369ePss8/mrbfeAmDPnj3ccMMNLF++nCFDhrBx40aXozfGPdbv73+CuiXzx/eyWb2zuFnrTOvalgcua3q+p3Xr1jFt2jTOOussbr31VqZMmcLzzz/P/Pnz6d27N+PGjQOgU6dOvPjiizz22GO8//77zRqrMQ059Dnui30zzdGSKausYc2uYtbmHWDrvhL2l1ZSWllD6/BQYluHkxLfhgFdYxnYLZbQEB88CX4mqJOMm5KSkjjrrLMAuOGGG3jyySdJSUkhNTX18LqpU6e6GaIJUr58Raqyuva4tss/UMF7K3by2bo9LNpccLieVmEhdIhqRetWoVRU1VBQWkl5lee99lHhjBrQmRtG9OAUm/H5uAV1kvGmxdFS6t9+XFRUZLckG5/gi9PJHFJV432SUVW+ztnHy99s4bN1e6ipVXp3iubGET0Y0TOO/l1i6Nau9RF/d7W1ys6iMr7dVshna/fw7vKdzFyynZF94/ndxf3pkxDTEocV0II6ybhp27ZtLFiwgDPOOIM333yTCy64gOeff56NGzfSq1cv3nzzTbdDNEHKd1OM9y2ZFdsLefTDtXyzcR/xMRFMOCeFa4Yl0rvT0ZNESIiQ2D6KxPZRXD64Kw9eXsWbi7fx9Gc5XPTEV9z94978z3m9CQu17mxvWZJxSf/+/ZkxYwZ33HEHqampPPHEEwwbNoxLLrmEjh07cvbZZ7Nq1Sq3wzRByIcbMlTVHD24jfkH+fvH65j73S7i2rTigcvSuP707kSEhR7X/mJbh3Pnj3oxLj2Jh95fzROfbuDrnL08f+Mw4qIjjqvOYGNJxiUhISE899xzR6wbPXo0a9eu/UHZkSNHMnLkyJMUmQl2PpxjqKiuObysqocvde0qKueJT9czOyuXyLAQfnlBKhPO6Ul0RPN8xLVv04rHxw1hZN947p2zkque/YaXbxlOSsc2zVJ/ILMkY4w5wqE+GfHBx5cVllYdXj5YUU1tLTzzRQ4vf72FWlVuHNGDu8/rTccWamWMGdKNxPZR3P5KFtdNXci/7jyDpA5RLbKvQGFJxgXJycl2Kcz4rEOXy3zxPpT9pZWHl//f7BUs3LSPAxXVXDmkG7+6sM9J+cAf1qM9r084nYypC7n+xYX8+84z6dQ2ssX366+CsvfKl++eORaBchzGt/jidDKHbNh9kIiwEHp3imbemt2c2asjc39+DlPGDTmpLYr+Xdryyq3D2XugkjtfW3rEZTxzpKBryURGRrJv3z6/n+7/0PNkIiPtG5RpXr763aWkoppP1+5mZN94nr5+KJU1tUS1cu8jbHBSOx67ZjCT3viWBzOz+b+rBrkWiy8LuiSTmJhIbm4u+fn5bodywg49GdOY5uSLOaayupZ756xkX0klt5/Tk7DQEJ+4jfiSQV1YtbMXz36+kXNT47loYBe3Q/I5XiUZERkNPAGEAi+q6l/qvR8BvAIMA/YB41R1i/PeZOA2oAb4uap+dLQ6RWQakI5nnr71wM2qelBEbgb+BuxwdvtPVX3xWA84PDzcniRpzFH4WkvmYEU1d722lK827OV3F/cjPbmD2yEd4Z4L+/B1zl5+9/Z3DEtuT6cYu7pQV5NfBUQkFHgauAhIA64TkbR6xW4D9qtqb+Bx4FFn2zQgAxgAjAaeEZHQJur8laoOVtVBwDbg7jr7maWqQ5yfY04wxpimHeqT8YWLyXsPVnDd1IV8s3Effx07iInn9nI7pB8IDw1hyrWDKa2s4Q/v2A099XnT3hwO5KjqJlWtBGYCY+qVGQPMcJbnAOeLp8NjDDBTVStUdTOQ49TXaJ2qWgzgbN8a32y9GxOwfKUls3VfCVc/+w0b9hzghZuGcW16ktshNap3pxh+fn4qH2Xv5vN1e9wOx6d4k2S6AdvrvM511jVYRlWrgSIg7ijbHrVOEXkJ2AX0A56qU+5qEVkpInNExHd/44wxJ2TVjiKufvYbisqqeOP2EZzXL8HtkJo04ZwUUjq24Y/vrba7zerwJsk01Gqu/12nsTLHut6zoHoL0BVYA4xzVr8HJDuX0T7h+5bTkYGITBSRLBHJCoTOfWNONrdvjf9qQz7jnl9ARFgoc+48k6Hd27saj7ciwkJ54LI0Nu8t4dUFW90Ox2d4k2RygbqthkRgZ2NlRCQMiAUKjrJtk3Wqag0wC7jaeb1PVSuct1/Ac5PBD6jqVFVNV9X0+Ph4Lw7PGFOXmznmnWU7uOWlJSR1iOKtn51J707R7gVzHEb27cQ5qR155vONlFRUux2OT/AmySwBUkUkRURa4enIz6xXJhMY7yyPBear5+tQJpAhIhEikgKkAosbq1M8esPhPpnLgLXO67r3Bl6Op5VjjGlmbuWYF77cxC9nLSc9uT2z7zyDBD8dRX/PhX0oKKnk5W+2uB2KT2jyFmZVrRaRu4GP8NxuPF1Vs0XkISBLVTOBacCrIpKDpwWT4WybLSKzgdVANTDJaaHQSJ0hwAwRaYvnktoK4C4nlJ+LyOVOPQXAzc1yBowxR/h+WpmTc39Zba3yyNw1TPvvZi4Z2IUp4wYf96zJvuDU7u05v18nnv9iIzee0YO2keFuh+Qqr8bJqOpcYG69dffXWS4Hrmlk20eAR7yssxY4q5F6JgOTvYnXGHP8Tua0MhXVNfz6Xyt5b8VObj4zmfsvTSMkAB55/MsL+nDZP//LrMXbuf3cnm6H4yr3h8waY3zKyeqTKa+qYeIrS3lvxU7uu6gfD1wWGAkGYGBiLCN6dmD615uP6WmegciSjDHmCIdyTEveZVZWWcOEGVl8uSGfv1w1kDt/1Muv5xJsyMRze5JXVM5/Vua5HYqrLMkYY47kJJeWSjElFdXc8vJivtm4l8fGDiZjePcW2pO7RvbpRK/4Nrz4302u3xbuJksyxpgjtOTH4YHyKsZPX8zizQU8Pm4IVw8L3AleQ0KEm89MZtWOYr7bUeR2OK6xJGOMOcKhL93N/eW7qKyKm6YvZtn2Qp66bihjhtSfOCTwjDm1G5HhIby5eHvThQOUJRljzBFqD18ua74sU1hayY3TFrFqRxFPXz+USwYFx5T4bSPDuWRgVzKX7wjawZmWZIwxR/i+47956isoqeT6FxaxNu8Az90wjNGndG6eiv3E9acnUVJZw/sr60+UEhwsyRhjjlDbjB3/ew9WcP0LC9mYf5AXxqdzfn/fn+iyuQ3t3p7enaKZszTX7VBcYUnGGHOEmhonyZxgltlTXE7G1IVs2VfC9JtP40d9gnMuQRHhiiFdWbJlPzsKy9wO56SzJGOMOUKNk11qa48/y+wq8iSYnYVlvHzLcM7q3bG5wvNLlw3uCsD7K4LvkpklGWPMEWqc5HK8I9V3FJYxbuoC9hyo4JVbhzOiZ1xzhueXesS1YXBSOzItyRhjgt3hJHMcLZntBaWMe34BBSWVvHrbcNKTOzR3eH7r8sFdyd5ZzMb8g26HclJZkjHGHOFQkqk+xpbMlr0ljHt+AQfKq3l9wumc6icPGztZLh3UBRF4f0VwTTNjScYYc4TjuVy2Mf8g46YuoKyqhjduP51Bie1aKjy/ldA2klOT2vHJmt1uh3JSWZIxxhzh+yTj3eWyDbsPkDF1ITW1ysyJZzCga2xLhufXLkzrzHc7isgrCp67zCzJGGOOcOjuMm9aMmt3FZMxdSEAMyeOoG/nmBaNzd9dmNYJgE9WB09rxpKMMeYI1Yf7ZI7eklm1o4jrpi4kPDSEWRNH0LuTJZim9IqPJqVjGz62JGOMCVaHOvwrj9KSWbG9kOtfWEhUqzBm3TGCnvHRJys8vyYiXJiWwMJN+ygur3I7nJPCkowx5ghlVZ7kUl3bcJJZunU/N7y4iNiocGZOHEGPuDYnMzy/d0H/BKpqlP9u2Ot2KCeFV0lGREaLyDoRyRGR+xp4P0JEZjnvLxKR5DrvTXbWrxORUU3VKSLTRGSFiKwUkTkiEt3UPowxzae8sgaA0oqaH7y3aNM+bpq2iA7RrZg18QySOkSd7PD83qnd2xEdEcZXlmQ8RCQUeBq4CEgDrhORtHrFbgP2q2pv4HHgUWfbNCADGACMBp4RkdAm6vyVqg5W1UHANuDuo+3DGNO8yqs9yWV/aeUR679cn8/4lxbTOTaSWRPPoGu71m6E5/fCQ0M4o1ccX23ID4onZnrTkhkO5KjqJlWtBGYCY+qVGQPMcJbnAOeL54HdY4CZqlqhqpuBHKe+RutU1WIAZ/vWfD8ZbGP7MMY0ozKnJVNYVnX4duZ5q3czYUYWKR2jmXXHGXSOjXQzRL93TmpHcveXsXVfqduhtDhvkkw3oO5j3XKddQ2WUdVqoAiIO8q2R61TRF4CdgH9gKea2McRRGSiiGSJSFZ+fr4Xh2eMqauozNMhrQp5RWW8smALd762lP5d2zLz9hF0jI5wN8AAcE6qZ0bqr3IC/5KZN0mmodZC/TZeY2WOdb1nQfUWoCuwBhh3DHGgqlNVNV1V0+Pjg3NqcWOOV22tsq+kkvP6ecZzjHt+Ife/m83IPvG8PuF0YqPCXY4wMCTHRdGtXWu+Wh/4X4S9STK5QFKd14lA/alED5cRkTAgFig4yrZN1qmqNcAs4Oom9mGMaSZ7SyqoqVV+1Cee64Z7/kR/O7ofU29KJzoizOXoAoeIcG6fjizYuO+Y54jzN9781iwBUkUkBdiBpyP/+nplMoHxwAJgLDBfVVVEMoE3RGQKnpZJKrAYT6vkB3U6fSy9VDXHWb4MWHu0fRzncRtjGrA27wAAqQnRjD8z2d1gAtyZvTry5uLtZO8sZnBS4M711mSSUdVqEbkb+AgIBaararaIPARkqWomMA14VURy8LQuMpxts0VkNrAaqAYmOS0UGqkzBJghIm3xJKIVwF1OKA3uwxjTfL7euJewELH5x06C01M8j0FYvLkgoJOMBHJjID09XbOystwOwxi/UFpZzY/+9jmndG3LS7cMdzucoDDyb5+RmhDDCzelux3KEURkqao2S1A24t8YA8CTn+aQf6CCu8/r7XYoQWN4SgeWbCk4oUdd+zpLMsYYPlu3h+e+2Mi49CSG9bCnWZ4sw1PiKCytYsOewH1apiUZY4Lc+t0H+OXM5fTrHMMfxwxwO5yg8n2/zD6XI2k5lmSMCWI7Csu4adpiWoWF8MJN6USGh7odUlBJbN+azm0jWbxlv9uhtBhLMsYEqYKSSm6atoiSimpm3DLcJrt0gYgwPKUDizfvC9h5zCzJGBOESiurufXlJWzfX8aL49NJ69rW7ZCC1rAe7dldXEFeUbnbobQISzLGBJnyqhrufO1bVuYW8tR1p3J6zx9MAWhOokNjZJZvL3Q5kpZhScaYIFJZXcvdb3zLl+vz+ctVgxg1oLPbIQW9/l1iaBUaYknGGOPfqmtq+cXMZXyyZg9/uuIUrj0tqemNTIuLCAslrWtbSzLGGP9VU6vcM3sFH6zaxR8uTePGET3cDsnUMSSpHd/lFgXkZJmWZIwJcLW1ym//vZLMFTv57eh+3HZ2itshmXqGJLWjrKqG9bsDb1CmJRljApiq8r/vrmLO0lx+dUEf7hrZy+2QTAOGBHDnvyUZYwKUqvLH91bzxqJt/GxkL35+vs1J5qt6xEXRPiqc5dsDb1CmJRljApCq8pcP1vLyN1uYcHYKvxnVF88jmowvEhEGJbZjxfYit0NpdpZkjAlAj89bz/NfbuKmM3rw+0v6W4LxA6d0a0tO/kHKq2rcDqVZWZIxJsD8c/4GnpyfQ8ZpSTx42QBLMH5iQNdYamqV9bsPuB1Ks7IkY0wAefqzHB77eD1XndqNP185kJAQSzD+YoAztU/2zmKXI2leTT5+2RjjH/45fwOPfbyeK0/txt+uGWwJxs8ktY8iJiKM7J2B1S/jVUtGREaLyDoRyRGR+xp4P0JEZjnvLxKR5DrvTXbWrxORUU3VKSKvO+tXich0EQl31o8UkSIRWe783H8iB25MIKmbYB67ZjChlmD8TkiI0L9r24BryTSZZEQkFHgauAhIA64TkbR6xW4D9qtqb+Bx4FFn2zQgAxgAjAaeEZHQJup8HegHDARaAxPq7OcrVR3i/Dx0PAdsTKB56tMNhy+RWYLxbwO6tmVt3gFqAuhxzN60ZIYDOaq6SVUrgZnAmHplxgAznOU5wPni6W0cA8xU1QpV3QzkOPU1WqeqzlUHsBhIPLFDNCZwPfnpBv4+z5Ng/mYJxu+ldWlLWVUNm/eWuB1Ks/EmyXQDttd5neusa7CMqlYDRUDcUbZtsk7nMtmNwId1Vp8hIitE5AMRafA5sSIyUUSyRCQrPz/fi8Mzxj89+ekGpsxbz1VDLcEEigFdYwECql/GmyTT0G9u/bZcY2WOdX1dzwBfqupXzutvgR6qOhh4CninoWBVdaqqpqtqenx8fENFjPF7T3xSJ8GMtQQTKFITomkVGsLqAOqX8SbJ5AJ15wRPBHY2VkZEwoBYoOAo2x61ThF5AIgH7jm0TlWLVfWgszwXCBeRjl7Eb0xAeeKTDTz+iSWYQBQeGkKfztGszguuJLMESBWRFBFphacjP7NemUxgvLM8Fpjv9KlkAhnO3WcpQCqefpZG6xSRCcAo4DpVPTzvtYh0dvp5EJHhTuz7juegjfFX//hkPY9/sp6rhyZagglQfRJiAmpAZpPjZFS1WkTuBj4CQoHpqpotIg8BWaqaCUwDXhWRHDwtmAxn22wRmQ2sBqqBSapaA9BQnc4unwO2AgucnPKWcyfZWOAuEakGyoAMJ5EZE/BUlSc+3cA/PtnA1UMT+evYQZZgAlTfhBje+nYHRaVVxEaFux3OCZNA/pxOT0/XrKwst8Mw5oSoKo99vI6nP9toCSYIfLZ2D7e8vIR/3XkGpyV3cCUGEVmqqunNUZdNK2OMD1NVHv7PGp7+bCPXDU/ib5ZgAl6fzjEArNsVGJfMbFoZY3xUba1yf+YqXlu4jZvPTOaBy9Jssssg0DU2kuiIMDYESL+MJRljfFBNrTL5rZXMzsrljh/15L7R/SzBBAkRITUhmnUBkmTscpkxPqa6ppZ7Zi9ndlYuvzg/1RJMEOqbEMOG3QfdDqNZWJIxxodUVtfyP28u493lO/nNqL786sI+lmCCUGpCDPtKKtl7sMLtUE6YJRljfER5VQ13vbaUD1bt4g+XpjHpx73dDsm4pG+Cp/N/fQB0/luSMcYHlFXWcPsrWXy6dg9/uuIUbjs7xe2QjIv6JEQDBMSgTOv4N8ZlJRXV3DZjCYs2F/DXsYO4Nj2p6Y1MQIuPiaBdVDjrAqBfxpKMMS4qLK3k5peW8N2OIv4xbghjhtSf4NwEIxGhV3w0m/L9P8nY5TJjXLKnuJxxzy9k9c5inv3pUEsw5gg9O7ZhUwA8V8aSjDEu2F5QytjnFrB9fykv33IaPxnQ2e2QjI/pGR9N/oEKDpRXuR3KCbEkY8xJtn73Aa5+9huKy6t44/YRnNnbnlhhfqhnfBsANuX7d2vGkowxJ9GK7YVc+/wCAGZNPIMhSe1cjsj4ql6Hksxe/+6XsY5/Y06Sbzbu5fYZWXSIbsXrt42ge1yU2yEZH9a9QxtCQ4SNe/y7JWNJxpiTYN7q3Ux641uS46J49bbTSWgb6XZIxse1CgshqX1ra8kYY47u7WW5/PpfKzmla1tevmU47du0cjsk4yd6xkdbn4wxpnEvf72ZX81awfDkDrx++whLMOaY9OzYhs17S6it9d+HS1qSMaYFqCp//XAtD763mgvTEnjpltOIjrALB+bY9IyPpqK6lh2FZW6Hcty8SjIiMlpE1olIjojc18D7ESIyy3l/kYgk13lvsrN+nYiMaqpOEXndWb9KRKaLSLizXkTkSaf8ShEZeiIHbkxLqaqp5TdzVvLM5xu5bnh3nv3pUCLDQ90Oy/ihw7cx+/GgzCaTjIiEAk8DFwFpwHUiklav2G3AflXtDTwOPOpsmwZkAAOA0cAzIhLaRJ2vA/2AgUBrYIKz/iIg1fmZCDx7PAdsTEsqraxm4itZzFnqeRbMn688hbBQu2Bgjs/3Y2X8t/Pfm9/+4UCOqm5S1UpgJjCmXpkxwAxneQ5wvngegjEGmKmqFaq6Gchx6mu0TlWdqw5gMZBYZx+vOG8tBNqJSJfjPG5jml1BSSXXv7CIL9bn8/AVp9izYMwJi4+OICYijM2B3JIBugHb67zOddY1WEZVq4EiIO4o2zZZp3OZ7Ebgw2OIAxGZKCJZIpKVn5/vxeEZc+Jy95cy9rlvWJ1XzDM/HcYNI3q4HZIJACJC97gotu4rdTuU4+ZNkmnoq1j9Wx0aK3Os6+t6BvhSVb86hjhQ1amqmq6q6fHx8Q1sYkzzWpNXzFXPfMPeAxW8dtvpjD7F5iEzzadHXBTbCgI7yeQCdR9wkQjsbKyMiIQBsUDBUbY9ap0i8gAQD9xzjHEYc1It3LSPa59fQIgI/7rzTIandHA7JBNgundoQ+7+Umr89DZmb5LMEiBVRFJEpBWejvzMemUygfHO8lhgvtOnkglkOHefpeDptF98tDpFZAIwCrhOVWvr7eMm5y6zEUCRquYdxzEb0yzmfpfHTdMX0ykmgn//7Ez6do5xOyQTgHrERVFVo+z009uYm7xxX1WrReRu4CMgFJiuqtki8hCQpaqZwDTgVRHJwdOCyXC2zRaR2cBqoBqYpKo1AA3V6ezyOWArsMDpNH1LVR8C5gIX47l5oBS4pTlOgDHHSlV54atN/HnuWob1aM+LN6XbIEvTYnp08Mxxt62glKQO/jffnVejw1R1Lp4P+brr7q+zXA5c08i2jwCPeFOns77BmJyW0SRv4jWmpVTX1PLge9m8tnAblwzswt+vHWxjYEyLOjSR6tZ9pZzV2+VgjoMNQTbGSyUV1dz9xrd8ti6fO37Uk9+O6kdIiN2ibFpWl9jWhIcKWwv88zZmSzLGeGF3cTm3vryENXnFPHzFKXaLsjlpQkOEpPZRbPPT25gtyRjThLW7irn1pSUUllUxbfxp/LhfJ7dDMkHGn8fK2HwXxhzFfzfs5ZpnF1Bdq8y+4wxLMMYVPTp4xsp4uqb9iyUZYxoxO2s7N7+0mK7tWvPOpLM4pVus2yGZINU9rg0HK6opKKl0O5RjZpfLjKmntlb5+7x1PP3ZRs5J7cjTPx1K28hwt8MyQezQbcxbC0qJi45wOZpjYy0ZY+ooqajmrteX8vRnG8k4LYnpN59mCca4rodzG7M/dv5bS8YYx87CMibMyGLtrmL+cGkat56VbLMoG59waBCmP3b+W5IxBli2bT+3v7KU8qoapt18Gj/uax38xndEhoeS0DaC7fstyRjjd95dvoPfzFlJQtsI3rj9dPok2Bxkxvckto9ix37/m7/MkowJWrW1yj8+Wc+T83MYntyB524cRgebg8z4qG7tWrNs+363wzhm1vFvglJpZTWT3viWJ+fncM2wRF6bcLolGOPTEtu3Jq+w3O+m/LeWjAk6OwvLmPhqFtk7i/n9xf2ZcE6KdfAbn9etfWuqa5XdxeV0bdfa7XC8ZknGBJVFm/Yx6Y1vKa+q5cWb0jm/f4LbIRnjlcT2njvMdhSW+VWSsctlJiioKq8u2MJPX1xE28hw3pl0piUY41e6OYkl18/uMLOWjAl4FdU13P9ONrOytnNev078I2OIDbA0ficovkUsAAAX20lEQVSxvSfJ+NsdZpZkTEDbVVTOna8tZfn2Qv7nvN786oI+9gwY45ciw0PpGN2KXEsyxviGpVsLuPO1bymtqOa5G4Yy+pQubodkzAnp1j6KHYX+lWS86pMRkdEisk5EckTkvgbejxCRWc77i0Qkuc57k53160RkVFN1isjdzjoVkY511o8UkSIRWe78HH78szH1vbFoGxlTF9KmVShvTzrLEowJCIntWgdeS0ZEQoGngQuBXGCJiGSq6uo6xW4D9qtqbxHJAB4FxolIGpABDAC6Ap+ISB9nm8bq/Bp4H/i8gXC+UtVLj+M4TZCoqK7hj++t5o1F2xjZN54nxp1KbJT1v5jAkNi+NfNW76a2Vv3msq83l8uGAzmquglARGYCY4C6SWYM8KCzPAf4p3gGHowBZqpqBbBZRHKc+misTlVd5qw7keMyQWhHYRk/e20pK3KLuGtkL379k76E+skfojHe6Na+NZU1tew9WEGntpFuh+MVby6XdQO213md66xrsIyqVgNFQNxRtvWmzoacISIrROQDERnQUAERmSgiWSKSlZ+f70WVJhB8sT6fS5/8ik35JTx3wzB+O7qfJRgTcA7dYbbdjy6ZeZNkGvpLrT+vQWNljnX90XwL9FDVwcBTwDsNFVLVqaqarqrp8fHxTVRp/N2h+cdufmkxCW0jyfyfsxl9Sme3wzKmRXRr9/2ATH/hzeWyXCCpzutEYGcjZXJFJAyIBQqa2LapOo+gqsV1lueKyDMi0lFV93pxDCYAFZRU8stZy/lyfT5XDe3GI1cMpHWrULfDMqbFdGvvfwMyvWnJLAFSRSRFRFrh6cjPrFcmExjvLI8F5quqOusznLvPUoBUYLGXdR5BRDo7/TyIyHAn9n3eHKQJPMu3F3LZU/9l4cZ9/PnKgfz9msGWYEzAi44Io11UuF8NyGyyJaOq1SJyN/AREApMV9VsEXkIyFLVTGAa8KrTsV+AJ2nglJuN5yaBamCSqtaA51bl+nU6638O3At0BlaKyFxVnYAned0lItVAGZDhJDITRFSV1xZt46H3skloG8m/7zqTgYmxbodlzEnTrV1rv7pcJoH8OZ2enq5ZWVluh2GaSUlFNf/7zireXraDH/eN5/FxQ2gXZdPzm+AyYUYWuftL+fCX57bYPkRkqaqmN0ddNuLf+IU1ecVMeuNbNu8t4f9d2IdJP+7tN+MEjGlOXWIjWbzZf3oKLMkYn6aqvLF4G398bzXtWofz+oTTObNXx6Y3NCZAdWkXSXF5NSUV1bSJ8P2PcN+P0ASt4vIqJr/1Hf9Zmce5feKZcu1gOkZHuB2WMa7qEusZhJlXVE7vTtEuR9M0SzLGJ63YXsj/vLmMHYVl/HZ0P+44t6ddHjMG6BLruY05r6jMkowxx0pVmfbfzTz64Vo6xUQy+44RDOvRwe2wjPEZXQ8nmXKXI/GOJRnjM/aXVPKbOSv4ZM0eLkxL4G9jB9ndY8bUkxDruWScV2hJxhivLd5cwC9mLmPvwQoeuCyNm89MtklSjWlARJjn4WW7iv1jrIwlGeOqqppanvhkA898nkNShyj+fdeZDEps53ZYxvi0LrGt2WktGWOObsveEn4xazkrthdyzbBEHrh8ANF+cEumMW7rHBvJ1n0lbofhFfuLNiedqvKvpbk8mJlNWIjw9PVDuWSQPbnSGG91jY1k4Sb/GJBpScacVIWllfzu7e+Y+90uRvTswJRrh9C1XWu3wzLGr3SObc2B8moOVlT7fOvft6MzAeWbjXv5f7NXkH+ggvsu6sft5/S0B4sZcxy6tvMMyNxVVEbvTjEuR3N0lmRMi6usrmXKvPU8/+VGUuLa8PbPzrKZk405AYcGZO4sLLckY4Lb+t0HuGf2clbtKOa64d35w6X9iWplv3bGnIjvp5bx/duY7a/dtIiaWmXafzfx2MfriYkI4/kbhzFqgD0W2ZjmkNA2EhH/GPVvScY0u637Svj1v1awZMt+Rg1I4JErB9rElsY0o1ZhIXSMjvCLUf+WZEyzUVVeX7SNP89dQ2iIMOXawVx5ajcbuW9MC+gSG0lesSUZEyR2FZVz779X8uX6fM5J7cijVw+yW5ONaUFdYiPZlO/7AzJDvCkkIqNFZJ2I5IjIfQ28HyEis5z3F4lIcp33Jjvr14nIqKbqFJG7nXUqIh3rrBcRedJ5b6WIDD3egzbNR1V5e1kuP3n8C5ZsLuBPYwbwyq3DLcEY08K6xLYOjD4ZEQkFngYuBHKBJSKSqaqr6xS7Ddivqr1FJAN4FBgnImlABjAA6Ap8IiJ9nG0aq/Nr4H3g83qhXASkOj+nA886/xqX7DtYwe/fXsWH2bsY1qM9f79mMMkd27gdljFBIaFtJAcrfH9ApjeRDQdyVHUTgIjMBMYAdZPMGOBBZ3kO8E/xXIgfA8xU1Qpgs4jkOPXRWJ2qusxZVz+OMcArqqrAQhFpJyJdVDXvWA7YnDhV5b2VeTyYmc3B8mobWGmMCzo7U/7vLi4nOt53H17mTZLpBmyv8zqXH7YgDpdR1WoRKQLinPUL623bzVluqk5v4ugGHJFkRGQiMBGge/fuTVRpjtWe4nJ+/84q5q3ezeDEWP46djB9O/v2YDBjAlFCW89Ymd3F5fTy8yTT0NdT9bJMY+sb6guqX+fxxIGqTgWmAqSnpzdVp/GSqjJnaS5/en81FdW1/O7iftx6VgphoV516xljmlndJOPLvEkyuUBSndeJwM5GyuSKSBgQCxQ0sW1TdR5PHKYF7CgsY/Jb3/Hl+nyGJ3fg0bGDSLG+F2NcdSjJ7CqqcDmSo/Pma+gSIFVEUkSkFZ6O/Mx6ZTKB8c7yWGC+03eSCWQ4d5+l4Om0X+xlnfVlAjc5d5mNAIqsP6Zl1dYqry3cyk+mfEHWlgIeGjOAmRNHWIIxxgdER4QRHRHm/y0Zp4/lbuAjIBSYrqrZIvIQkKWqmcA04FWnY78AT9LAKTcbz00C1cAkVa0Bz63K9et01v8cuBfoDKwUkbmqOgGYC1wM5AClwC3NdRLMD23dV8Jv/72ShZsKOLt3R/7vqoEkdYhyOyxjTB0JbSN8PsmIp8ERmNLT0zUrK8vtMPxKdU0t07/ezJR56wkPCeF/L+3PtelJNmrfGB90/QsLKa+q4a2fndWs9YrIUlVNb466fPfmanPSrdheyOS3vmN1XjEX9E/g4StOobMz26sxxvd0bhvJos0FbodxVJZkDAcrqnnso3W8smAL8TERPHfDUEYN6GytF2N8XEJsJHsOlFNbq4T46Dg1SzJB7uPsXTyQmc2u4nJuHNGDX4/qS9vIcLfDMsZ4ISEmgqoapaC00mdnOrckE6Tyisp44N1sPl69m36dY3j6p0MZ2r2922EZY47BocvZu4rKLckY31BTq7y6YAuPfbye6tpafju6HxPOSSHcBlUa43c6OWNl9hwoxzM80fdYkgkiq3YU8ft3VrFieyHnpHbkkSsG0j3Obks2xl919oMBmZZkgkBRWRVTPl7Hqwu30qFNK57IGMLlg7tax74xfi4+JgIR355axpJMAFNV3vp2B//3wRoKSiq5cUQP7vlJX2JbW8e+MYEgPDSEuDa+PSDTkkyAWpNXzP3vrmLJlv2c2r0dL98ynFO6+eY1W2PM8esca0nGnEQHyqt4fN4GZizYQmzrcP569SDGDkv02XvojTEnJiEmkp0+/IRMSzIBQlXJXLGTh/+zhr0HK7h+eHd+M6ov7aJauR2aMaYFJcRGsmx7odthNMqSTABYu6uYBzOzWbipgEGJsbx4UzqDk9q5HZYx5iRIiImkoKSSiuoaIsJC3Q7nByzJ+LGCkkoen7ee1xdtJSYynEeuPIWM07rbY5CNCSKHHsO8p7jCJ2dKtyTjh6pqanlt4VYen7eeksoabhzRg19e0If2bezSmDHBpu6ATEsy5oR9sT6fP72/mpw9BzkntSN/uDSNPgkxbodljHGJrw/ItCTjJzbvLeHh91fz6do9JMdF8cJN6VzQv5MNqDQmyB1KMr56G7MlGR9XXF7FP+fn8NLXm4kIC2XyRf24+axkn+zgM8acfO2iwmkVFmJJxhybqppa3ly8jSc/3cC+kkquGZbIr0f1pVOMPUTMGPM9ESGhbQS7fDTJeDX1roiMFpF1IpIjIvc18H6EiMxy3l8kIsl13pvsrF8nIqOaqlNEUpw6Njh1tnLW3ywi+SKy3PmZcCIH7qtUlfdX7uTCKV9w/7vZ9IqPJnPS2fx17GBLMMaYBiXERPpvS0ZEQoGngQuBXGCJiGSq6uo6xW4D9qtqbxHJAB4FxolIGpABDAC6Ap+ISB9nm8bqfBR4XFVnishzTt3POtvMUtW7T/CYfdY3G/fylw/WsjK3iL4JMbx082mM7Btv/S7GmKNKiI1k9c5it8NokDeXy4YDOaq6CUBEZgJjgLpJZgzwoLM8B/ineD4ZxwAzVbUC2CwiOU59NFSniKwBzgOud8rMcOo9lGQC0pq8Yv7ywVq+WJ9P19hIHrtmMFee2s3GuxhjvJIQE8lnxXtQVZ/7UupNkukGbK/zOhc4vbEyqlotIkVAnLN+Yb1tuznLDdUZBxSqanUD5QGuFpFzgfXAr1S1bh1+J3d/KVM+Xs/by3fQNjKc313cj5vOSCYy3Dr1jTHe6xwbQWllDQcqqn3u8eneJJmG0qJ6Waax9Q31BR2tPMB7wJuqWiEid+Jp5Zz3g2BFJgITAbp3795Ade7bX1LJM5/nMOObrSAw8dye/OxHvYmN8q1fDmOMf0g4NCCzuNwvk0wukFTndSKws5EyuSIShuc5oAVNbNvQ+r1AOxEJc1ozh8ur6r465V/A03fzA6o6FZgKkJ6eXj8Zuqq8qoaXvt7CM5/ncLCimrFDE/nVhX3o2q6126EZY/xYQp0Bmb07+dbgbG+SzBIgVURSgB14OvKvr1cmExgPLADGAvNVVUUkE3hDRKbg6fhPBRbjabH8oE5nm8+cOmY6db4LICJdVDXP2d/lwJrjPOaTrqZW+ffSXKbMW8+u4nLO79eJe0f3o29n3/plMMb4p8Oj/n3wDrMmk4zTx3I38BEQCkxX1WwReQjIUtVMYBrwqtOxX4AnaeCUm43nJoFqYJKq1gA0VKezy98CM0XkYWCZUzfAz0XkcqeeAuDmEz76FqaqfLpmD49+uJYNew4yJKkd/8gYwoiecW6HZowJIAk+POpfVH3qilKzSk9P16ysLFf2vXTrfh79YC2LtxSQ0rEN947qy+hTOvvcnR/GmMAw5KGPuXRQFx6+YuAJ1yUiS1U1vRnCshH/zW1j/kH+9uE6PszeRcfoCB6+4hTGnZZEeKhX416NMea4dIltTV6h77VkLMk0k70HK3jikw28sXgbkWEh3HNhH247O4U2EXaKjTEtr0tsJHk++Bhm+wQ8QWWVNUz/ejPPfr6Rsqoarh/enV9ckErH6Ai3QzPGBJEusZEs27bf7TB+wJLMCfgoexd/zMxmZ1E5F6YlcN9F/egVH+12WMaYINQlNpL9pVWUV9X41IBuSzLHYWdhGfe/m80na3bTr3MMU8bZHWPGGHd1ifWMt8srKielYxuXo/meJZlj9FH2Lu6ds5LK6lomX9SPW89OsU59Y4zrusR6bmPOKyqzJOOPVJUp89bz1PwcBnaL5anrTiXZh/4jjTHBrYszc4iv3WFmScYLqsr/vrOK1xdtY1x6Eg9dMcCeTGmM8Sm+OurfkowXnpqfw+uLtnHHj3py3+h+NqDSGONzWrcKpX1UODsLy9wO5QjWmdCEhZv2MWXeeq46tZslGGOMT+sc25pdPjZWxpLMUVTV1DL5re/o3iGKh688xRKMMcandY2NZKclGf/xzrIdbN5bwv2XphHVyq4sGmN8W5d2keQV2eUyvzH96y3079KW8/t3cjsUY4xpUpfY1hSWVlFWWeN2KIdZkmlEzp4DrMkr5tr0RLtMZozxC3XHyvgKSzKN+HTNHgAuHtjF5UiMMcY7h56ym7vfkozPW7p1P8lxUYcfBmSMMb7u0Ej/LftKXI7ke5ZkGrEit5BTu7d3OwxjjPFap5gIolqFsnmvJRmfVl5Vw+7iCp+a/8cYY5oiIiTHtbEk4+sOXc9M6tDa5UiMMebYpHRswxZ/SzIiMlpE1olIjojc18D7ESIyy3l/kYgk13lvsrN+nYiMaqpOEUlx6tjg1NmqqX00t4KSSgDio60/xhjjX5I7RrF9fxlVNbVuhwJ4kWREJBR4GrgISAOuE5G0esVuA/aram/gceBRZ9s0IAMYAIwGnhGR0CbqfBR4XFVTgf1O3Y3uoyWUVlYDnrmAjDHGn/Tv0paaWiV7Z7HboQDetWSGAzmquklVK4GZwJh6ZcYAM5zlOcD54hlcMgaYqaoVqroZyHHqa7BOZ5vznDpw6ryiiX00u/Iqz0CmKEsyxhg/c3qK5wGKH67a5XIkHt7MldIN2F7ndS5wemNlVLVaRIqAOGf9wnrbdnOWG6ozDihU1eoGyje2j711AxGRicBE5+VBEdlXv4y30lqsreSajhznuQhAdi487Dx8L6DOxeRHYfLxbdoR6NFccXiTZBpqLaiXZRpb31AL6mjlvY0DVZ0KTD0cmEiWqqY3sG3QsXPxPTsXHnYevmfnwsM5D8nNVZ83l8tygaQ6rxOBnY2VEZEwIBYoOMq2ja3fC7Rz6qi/r8b2YYwxxkd5k2SWAKnOXV+t8HTkZ9YrkwmMd5bHAvNVVZ31Gc6dYSlAKrC4sTqdbT5z6sCp890m9mGMMcZHNXm5zOn/uBv4CAgFpqtqtog8BGSpaiYwDXhVRHLwtC4ynG2zRWQ2sBqoBiapag1AQ3U6u/wtMFNEHgaWOXXT2D68MLXpIkHDzsX37Fx42Hn4np0Lj2Y9D2KNAWOMMS3FRvwbY4xpMZZkjDHGtJiATjJNTYfj70RkuojsEZFVddZ1EJF5zrQ880SkvbNeRORJ51ysFJGhdbYZ75TfICLjG9qXrxORJBH5TETWiEi2iPzCWR9U50NEIkVksYiscM7DH531xzxdU2NTQvkbZ5aRZSLyvvM6KM+FiGwRke9EZLmIZDnrWv7vQ1UD8gfPDQUbgZ5AK2AFkOZ2XM18jOcCQ4FVddb9FbjPWb4PeNRZvhj4AM94oxHAImd9B2CT8297Z7m928d2HOeiCzDUWY4B1uOZsiiozodzPNHOcjiwyDm+2UCGs/454C5n+WfAc85yBjDLWU5z/mYigBTnbynU7eM7znNyD/AG8L7zOijPBbAF6FhvXYv/fQRyS8ab6XD8mqp+yQ/HCtWdfqf+tDyvqMdCPOORugCjgHmqWqCq+4F5eOaZ8yuqmqeq3zrLB4A1eGaJCKrz4RzPQedluPOjHPt0TY1NCeVXRCQRuAR40Xl9PFNXBcS5aESL/30EcpJpaDqcbo2UDSQJqpoHng9eoJOzvrHzEXDnybnMcSqeb/FBdz6cy0PLgT14PgQ24uV0TUDdKaH8+jw4/gHcCxyaktjrqasIvHOhwMcislQ802/BSfj78GZaGX/l1TQ0QeRYp/7xSyISDfwb+KWqFkvjc6gG7PlQz1i0ISLSDngb6N9QMeffgD0PInIpsEdVl4rIyEOrGyga8OfCcZaq7hSRTsA8EVl7lLLNdi4CuSXjzXQ4gWi306zF+XePs/5Yp/jxOyISjifBvK6qbzmrg/Z8qGoh8Dmea+rHOl1TIJyHs4DLRWQLnsvl5+Fp2QTjuUBVdzr/7sHz5WM4J+HvI5CTjDfT4QSiutPv1J+W5ybnrpERQJHTPP4I+ImItHfuLPmJs86vONfOpwFrVHVKnbeC6nyISLzTgkFEWgMX4OmfOtbpmhqbEspvqOpkVU1Uz2SPGXiO7acE4bkQkTYiEnNoGc/v9SpOxt+H23c8tOQPnjsk1uO5Jv17t+NpgeN7E8gDqvB8w7gNzzXkT4ENzr8dnLKC50FxG4HvgPQ69dyKpzMzB7jF7eM6znNxNp5m+0pgufNzcbCdD2AQnumYVjofIvc763vi+WDMAf4FRDjrI53XOc77PevU9Xvn/KwDLnL72E7wvIzk+7vLgu5cOMe8wvnJPvR5eDL+PmxaGWOMMS0mkC+XGWOMcZklGWOMMS3GkowxxpgWY0nGGGNMi7EkY4wxpsVYkjHGGNNiLMkYY4xpMf8fCbWrKjzpHYAAAAAASUVORK5CYII=\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAD8CAYAAACl69mTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl4VdW5+PHvS+YQSCAJJJAAAQIYBBHCJNaLVQuOtBYrjsggXivXtvZXh3tvbeutvdp6tdWqiKKiBYGitqkTdawTUxgljGEOUyCEhBAynOT9/XE2EGJCDpCwz/B+nieP+6yz9trv2ZLzZu219tqiqhhjjDEtoZXbARhjjAlelmSMMca0GEsyxhhjWowlGWOMMS3GkowxxpgWY0nGGGNMi/EpyYjIaBHZICL5IvJgA+9Hichc5/3FItKtznsPOeUbRGTUabT5jIiU+XIMY4wx/qnJJCMiYcCzwJVAFnCTiGTVqzYJKFbVnsBTwOPOvlnAOKAvMBp4TkTCmmpTRLKBBF+OYYwxxn/50pMZAuSr6hZVrQLmAGPq1RkDzHS25wOXiYg45XNUtVJVtwL5TnuNtukkoD8A9/t4DGOMMX4q3Ic6nYGddV4XAEMbq6OqHhEpARKd8kX19u3sbDfW5lQgR1X31MshjR3jQN1KIjIFmALQunXrQX369PHhIxpjjtl/uJK9pRUkxUWRGh/tdjgn8dQq6/aUAtCvc7zL0QSvZcuWHVDV5OZoy5ck01Bvof5aNI3Vaay8oR6Uikgn4AZg5BnGgapOB6YDZGdna25ubgO7GWMaM+1fm3ns/fVMuaQ7/3nVeW6Hc5L9hysZ/OhHAOQ+drXL0QQvEdneXG35crmsAEiv8zoN2N1YHREJB+KBg6fYt7HyC4GeQL6IbANiRSS/iWMYY5qRP1+D1m//XWn8nC9JZimQKSIZIhKJdyA/p16dHGC8sz0W+ES9K2/mAOOcmWEZQCawpLE2VfVdVU1R1W6q2g0odwb6T3UMY0yosN/4gNPk5TJn/GMqsAAIA15W1TwReQTIVdUcYAbwutPrOIg3aeDUmwesBTzAPapaA9BQm02E0uAxjDGhw3JM4PFlTAZVfQ94r17Zw3W2K/COpTS076PAo7602UCdOF+OcTqqq6spKCigoqLibJtyXXR0NGlpaURERLgdiglC/nihwA9DMk3wKckEk4KCAtq0aUO3bt0I5BnQqkpRUREFBQVkZGS4HY4x54SNyQSekFtWpqKigsTExIBOMAAiQmJiYlD0yIzxlfVkAk/IJRkg4BPMMcHyOYx/8ed/VpZjAk9IJhljTGDyx3Eic2qWZPzUZ599xjXXXANAZWUll19+OQMGDGDu3LkuR2ZChT9+n/tjTObUQm7gPxCtWLGC6upqVq5c6XYoJgSIH9+OaUkm8FhPxgXbtm2jT58+jB8/nv79+zN27FjKy8v54IMP6NOnDxdffDFvvfUWAIWFhdx6662sXLmSAQMGsHnzZpejN6HCH7/PPbW1bodgTlNI92R+84881u4ubdY2szq15VfX9m2y3oYNG5gxYwYjRoxg4sSJPPnkk7zwwgt88skn9OzZkxtvvBGADh068NJLL/HEE0/wzjvvNGusxgSamlp/TH3mVEI6ybgpPT2dESNGAHDrrbfy9NNPk5GRQWZm5vGy6dOnuxmiMX7H0wxJ5nBFNSt2HGLjvsNs3n+EkqNVlFXWEBXeirbREXRPbs15qW0YmpFI6yj7ijxbIX0GfelxtJT6049LSkpsSrLxC8f+Gfrj+MeZ9mQKSyt4e8UuPl5XyPIdxceTVfvWkbSLjSAuKpxKTy2Hyqt5c3kBAJFhrfhOZhK3De/KJZnJtGplv59nIqSTjJt27NjBwoULGT58OG+88QaXX345L7zwAps3b6ZHjx688cYbbodoQtSx5OKPf/OcTk9GVfkqv4gZX27hXxv3U6twfue2TLmkOxf1SOK81DYkxkV9a7/DFdWsLijhk/WF5KzazR2vLCUrtS2/vCaL4T0Sm/PjhARLMi4577zzmDlzJnfddReZmZn86U9/YtCgQVx99dUkJSVx8cUXs2bNGrfDNCHMD3MMNT4M/KsqC/L28dxn+awuKCG5TRR3j+zBDwem0T05rsn920RHMKJnEiN6JvHA6D68s3o3//fPjdz04iJ+lJ3Gr6/rS2ykfXX6ys6US1q1asW0adNOKhs9ejTr16//Vt2RI0cycuTIcxSZMf6r0nPqJLNy5yF++85acrcXk5HUmseu78cPBnYmKjzsjI4XGd6K6wemcVW/VJ7+eBPP/2szuduLeXn8YLoltT6jNkONJRljzEmOLULpj5fLyio8ANQfHik+UsX/vLuWt5bvIikuiseu78cN2emENdM4SnREGPeP7sN3MpO5Z/Zyrn/+a2aMz+bCLu2apf1gZvfJuKBbt252Kcz4rRNjMv6XZcoqvUmmdZ3LVe99s4crnvoXOSt3c/fIHnz2i5GMG9Kl2RJMXcN7JPLm3RcRFxXObTOW8E1BSbMfI9iEZJIJlvWPguVzGOOrA2WVAERFhFFe5eG+uSv58azlpMRHkzP1Yh4Y3Ye4Fp52nJHUmrl3DSM+JoLbX15MfuHhFj1eoAu5JBMdHU1RUVHAf0Efe55MdHS026GYIOPPvxkb95UB3mTz/We/4u2Vu/jJZZn87ccjyOrU9pzFkRofw+w7hxLWqhWTZ+ZScrT6nB070ITcmExaWhoFBQXs37/f7VDO2rEnYxrTnPz1769D5VV8tG7f8df7D1fy2sQhfCcz2ZV4uia25vlbB3LT9EX8bO5KXro92+6laYBPSUZERgN/AsKAl1T1sXrvRwGvAYOAIuBGVd3mvPcQMAmoAe5V1QWnalNEZgDZeGdQbgTuUNUyEbkD+AOwyznsn1X1pdP9wBEREfYkSWNO4fjAv8tx1HXwSBUTXl3K4QoPsyYPZfHWg9w4OJ3OCTGuxjW4W3sevjaLh/+ex6tfb2PixfbdUl+Tl8tEJAx4FrgSyAJuEpGsetUmAcWq2hN4Cnjc2TcLGAf0BUYDz4lIWBNt/kxVL1DV/sAOYGqd48xV1QHOz2knGGNM0/ytJ7PzYDljn/+a9XtKef6WgYzomcR9V/RyPcEcc9uwrlzWpwOPf7CezfvL3A7H7/gyJjMEyFfVLapaBcwBxtSrMwaY6WzPBy4T79SUMcAcVa1U1a1AvtNeo22qaimAs38M/n2J2BjTgtbuLuX657/mQFklsyYP5Xt9U9wO6VtEhP+9vh8xkWH8fN4qam0Rz5P4kmQ6AzvrvC5wyhqso6oeoARIPMW+p2xTRF4B9gJ9gGfq1PuhiKwWkfkiku5D7MaY03R8UozL18sWbi7ixhcWEt5KmH/3RWR3a+9uQKfQoW00D1+Txcqdh5i/rMDtcPyKL0mmoX9q9VN1Y3VOt9y7oToB6ASsA250iv8BdHMuo33EiZ7TyYGITBGRXBHJDYbBfWPOtdrjOca9LPPeN3sY//ISOsZH8+bdF9GrYxvXYvHVDy7szKCu7fj9gvWUVthss2N8STIFQN1eQxqwu7E6IhIOxAMHT7Fvk22qag0wF/ih87pIVSudt1/EO8ngW1R1uqpmq2p2crI7s06MCWRuj8m8vnAb98xeTr+0eOb/+3A6+cnYS1NEhF9f25eiI1X8+ZN8t8PxG74kmaVApohkiEgk3oH8nHp1coDxzvZY4BP19rlzgHEiEiUiGUAmsKSxNsWrJxwfk7kWWO+8Tq1zvOvw9nKMMc3MrWVlVJUn/7mBX/49j8v6dOAvk4aSEBt5boM4S/3S4rn+wjRmfr2NwtIKt8PxC00mGWeMZSqwAO8X+zxVzRORR0TkOqfaDCBRRPKB+4AHnX3zgHnAWuAD4B5VrWmsTbyX0WaKyDfAN0Aq8IhzjHtFJE9EVgH3Anec9ac3xnyLGz2ZmlrlP99ew9Of5POj7DSm3TqImMgzW9TSbT+5LBNPrfLcZ/aodPDxPhlVfQ94r17Zw3W2K4AbGtn3UeBRH9usBUY00s5DwEO+xGuMOXPHcsy5SjYV1TX8ZM4KFuTt455Le/D/vtfbL9dN81WXxFhuGJTG7MU7uOvfupMaHxiX+1pKyC0rY4xpwjnsypQcreb2l5ewIG8fv7o2i1+M6hPQCeaYey7tSY0qM77Y6nYorrMkY4w5yfGeTAvfolZYWsGNLyxkxY5i/jRuABNGBM/d8untY7m6Xypzlu7kcIjPNLMkY4w5Se2xnkwL5pidB8sZO20hOw6WM2P8YMYMqH/rXeCb/J0Myio9zF26s+nKQcySjDHmJC19tWzrgSPc+MJCDpVXMWvyUC7pFZy3GvRPS2BIt/a88tU2PDVNPzY6WFmSMcacROv9tzlt3HeYH72wkApPLW9MGRb0T5aceHE3dh06ymcbQvfGcEsyxpiTHL9a1sxdmjW7Shg3fREAc6cMo2+n+GZt3x9ddl5HkuKimJsbupfMLMkYY07SEgP+K3YUc/OLi4gOb8W8u4aTGQDLxDSHiLBWjB2UxifrC0P25kxLMsaYkx3vyTRPc4u3FHHrS4tJiI1k3r8PJyOpdfM0HCBuHJxOTa0yf3loLpxpScYYc5Jjs8uaI8d8uekA419ZQkp8NPPuGk5au9hmaDWwZCS1ZmhGe/6aWxDwj30/E5ZkjDEn8TjLMJ/t9+GXmw4waeZSuiW2Zu5dw0mJj26G6ALTDwemsfXAEdbsKnU7lHPOkowx5iQ1zfDQra/yvQkmI6k1s+8cRlJcVDNEFrhG9U0hIkz4x+r6C9gHP0syxpiTHO/JnOEFs6/rJJhZk4fSvnVgraTcEuJjI7gkM5l3Vu0OuSdnWpIxxpykpubML5d9nX+AiTOX0rW9N8EkhngPpq5rL+jE7pIKlu8odjuUc8qSjDHmJDVnOBizcHMRE2cupUv7WGbdaQmmvsuzOhIV3oqcVaF1ycySjDHmJGcyJrNoSxETX11KertYG4NpRFxUOCN7J/PPvH0hNcvMkowx5iTHxmRqffwiXLSliAmvLCWtXYwlmCZcfl5H9pZWkLc7dGaZWZIxxpykpta7mKPHhx7NYifBdHYSTHIbSzCncmmfDojAR+v2uR3KOWNJxhhzEo8z8F/tOfXKwct3FDPh1aV0Sohm9p1DLcH4ICkuioFd2lmSMcaErqqapnsyebtLuOPlJSS3ieKNO4fRoU3o3mh5ui4/ryNrdpWyp+So26GcEz4lGREZLSIbRCRfRB5s4P0oEZnrvL9YRLrVee8hp3yDiIxqqk0RmSEiq0RktYjMF5G4po5hjGk+R6tqgBPJpr78wsPcPmMJcVHhzJo8lA5tLcGcjsvP6wDAx+sKXY7k3GgyyYhIGPAscCWQBdwkIln1qk0CilW1J/AU8LizbxYwDugLjAaeE5GwJtr8mapeoKr9gR3A1FMdwxjTvI5We5NMQ5fLdhSVc8tLixERZt05LCTXIjtbPTvEkdYuhs83hsYzZnzpyQwB8lV1i6pWAXOAMfXqjAFmOtvzgctERJzyOapaqapbgXynvUbbVNVSAGf/GE6s09fYMYwxzehYT6a6Xk9m96Gj3PzSIio9tcyaPDTkVlNuLiLCdzKTWLilKCSemOlLkukM1H3iToFT1mAdVfUAJUDiKfY9ZZsi8gqwF+gDPNPEMU4iIlNEJFdEcvfvD42/FIxpTuXHk8yJMZn9hyu59aXFlJRX8/rEofROCY3nwbSUET2TOFzhYfWuErdDaXG+JJmGegv1RwQbq3O65d4N1QlAJ2AdcONpxIGqTlfVbFXNTk4OzmeHG9NSVJWiI5XAiTGZQ+VV3DZjMXtKKnh5wmD6pQX/Ey1b2ogeSYh4V6oOdr4kmQIgvc7rNKD+ugjH64hIOBAPHDzFvk22qao1wFzgh00cwxjTTMoqPVRUe5NLeZWHwxXVjH95CVv2H+HF27MZ3K29yxEGh3atIzm/U7wlGcdSIFNEMkQkEu9Afk69OjnAeGd7LPCJetdNyAHGOTPDMoBMYEljbYpXTzg+JnMtsL6JYxhjmsmOg+WAdwmULfuPMOnVXPJ2l/LcLQO5ODPJ5eiCy4ieSSzfUcyRSo/bobSoJpOMM/4xFViA9/LVPFXNE5FHROQ6p9oMIFFE8oH7gAedffOAecBa4APgHlWtaaxNvJfEZorIN8A3QCrwyKmOYYxpPmucMYKr+6VSXlXDkm0H+eO4AVye1dHlyILPdzKT8NQqi7cWuR1Kiwr3pZKqvge8V6/s4TrbFcANjez7KPCoj23WAiMaaafRYxhjmseCvH2kxkfz0FV9iIkMY2TvZEb27uB2WEFpUNd2RIQJi7ce5Lt9gjeJ+5RkjDHBb82uEj7dUMg9I3uSEBvJr6/r63ZIQS06Ioz+aQks3RrcQ8u2rIwxhorqGh54czUJMRHceUl3t8MJGYO7teebXSVUODfABiNLMsaEOFXlN/9YS97uUp644QLiYyLcDilkDO7WjuoaZcWOQ26H0mIsyRgT4l74fAtvLNnB3SN7cNl5wTs24I+yu7ZHBJZuC95LZpZkjAlh83J38tj767n2gk784nu93Q4n5MTHRtC7YxtLMsaY4PPh2n089NY3fCczif+74QJatbKlAN0wuFt7lm8vDtp1zCzJGBOClmw9yNTZyzm/U1um3TqIyHD7KnDL4Iz2HKmqYd2ew26H0iLsX5YxIWblzkNMfHUpnRNiePmOwbSOsjsZ3HRhegIAK3cWuxxJy7AkY0wIWbOrhNtnLKZd6whm3TmUxDh7ZLLb0trFkBQXycqdwbkisyUZY0LE+r2l3DpjMW2iI5g9eRip8TFuh2TwPl/mgrQEVhUE5zRmSzLGhID8wsPc8uJiosJbMfvOoaS3tyda+pML0hPYvL+M0opqt0NpdpZkjAlyWw8c4eYXvY9Mnn3nMLom2hMt/c0F6QmowjcFwXfJzJKMMUFsR1E5N7+4CE+tMvvOofRIjnM7JNOAC5wHwa3cGXyXzCzJGBOkCorLuenFRRytruEvk4bSq6M9MtlfJcRGkpHUmlWWZIwxgWBvSQU3v7iY0opq/jJpKFmd2rodkmnCBWnxQTn4b0nGmCBTeLiCm19cxMEjVbw2cQjnd453OyTjgwvSE9hXWsm+0gq3Q2lWlmSMCSL7D1dyy4uL2VtawSsTBnNhl3Zuh2R81LeT94+BtbtLXY6keVmSMSZIFB6u4KYXF1FQfJQZ4wczuFt7t0Myp+G8VO+YWd7u4Jph5lOSEZHRIrJBRPJF5MEG3o8SkbnO+4tFpFud9x5yyjeIyKim2hSRWU75GhF5WUQinPKRIlIiIiudn4cxxgBQWFrBTdMXsav4KK9MGMzwHoluh2ROU5voCLomxrJ2T4j1ZEQkDHgWuBLIAm4Skax61SYBxaraE3gKeNzZNwsYB/QFRgPPiUhYE23OAvoA/YAYYHKd43yhqgOcn0fO5AMbE2wKSysY9+Ii9pRU8OqEwQzrbgkmUGWltg3Jy2VDgHxV3aKqVcAcYEy9OmOAmc72fOAyERGnfI6qVqrqViDfaa/RNlX1PXUAS4C0s/uIxgSvfaUVjJu+iL0lFbw6YQhDLcEEtKzUtmwrKqes0uN2KM3GlyTTGdhZ53WBU9ZgHVX1ACVA4in2bbJN5zLZbcAHdYqHi8gqEXlfRPo2FKyITBGRXBHJ3b9/vw8fz5jAtM+5RLavtIKZE4cwJMPGYALdsanm64LokpkvSaahJxmpj3VOt7yu54DPVfUL5/VyoKuqXgA8A/ytoWBVdbqqZqtqdnJyckNVjAl4e0u8PZhjCcYG+YNDMM4w8yXJFADpdV6nAbsbqyMi4UA8cPAU+56yTRH5FZAM3HesTFVLVbXM2X4PiBCRJB/iNyao7Ck5yrjpC9l/uJLXJg0h2xJM0OjYNor2rSNDLsksBTJFJENEIvEO5OfUq5MDjHe2xwKfOGMqOcA4Z/ZZBpCJd5yl0TZFZDIwCrhJVY8/j1REUpxxHkRkiBN70Zl8aGMC1e5DRxk3fREHyqqYOXEIg7paggkmIkJWalvy9gTPNOYmH4mnqh4RmQosAMKAl1U1T0QeAXJVNQeYAbwuIvl4ezDjnH3zRGQesBbwAPeoag1AQ206h5wGbAcWOjnlLWcm2VjgbhHxAEeBcU4iMyYkHEswxUeqeG3SEAbajZZBKatTW179ahuemlrCwwL/VkYJ5u/p7Oxszc3NdTsMY87azoPl3PzSIg4dqea1SUPsTv4gNn9ZAf/vr6v4+Of/5tqq2SKyTFWzm6OtwE+TxgS5rQeOcOMLCykpr+b1yUMtwQS5Xh29iWXTvsMuR9I8LMkY48c27TvMj15YSIWnljlThjMgPcHtkEwL69nBm2Q27itzOZLm0eSYjDHGHXm7S7htxhLCWglzpwwj054HExJiI8NJbx/DRuvJGGNaysqdh7hp+iKiw1sx767hlmBCTK8ObdgUJD0ZSzLG+Jml2w5y60uLSYiNZO5dw8lIau12SOYcy+zYhi0HyqiuqW26sp+zJGOMH/kq/wC3z1hCh7ZRzLtrOOntY90OybigV8c4qmuUbQeOuB3KWbMkY4yf+HR9IRNeXUrXxFjmThlOSny02yEZl/RyLo8Gw+C/JRlj/MAHa/Yy5fVcenWM4407h5HcJsrtkIyLenaIo5UQFIP/NrvMGJf9bcUufv7XVVyQFs8rE4YQHxPhdkjGZdERYXRpH8umwsBPMtaTMcZFry3cxk/nrmRIt/a8NmmoJRhzXGbHNna5zBhzZlSVZz7exMN/z+OKrI68MmEwcVF2YcGc0LNDHNuLjuAJ8BlmlmSMOcdqa5XfvruO//twI9cP7MzztwwkOiLM7bCMn8lIak11jVJQfNTtUM6KJRljziFPTS33v7maGV9u5Y6LuvHE2AuCYqVd0/x6JHvvj9pyILAvmdm/bmPOkYrqGu6ZvZz5ywr42eW9+NW1WbRq1dBDYo2BjCTvGmZb9gf2vTJ2EdiYc6Cs0sNdr+fyVX4Rv7o2iwkjMtwOyfi59q0jSYiNYEuA35BpScaYFlZ8pIo7Xl3Kml0lPPmjC7h+YJrbIZkAkZHUmq0B3pOxy2XGtKB9pRXcOH0h6/aUMu3WQZZgzGnpnhRnYzLGmIblF5Zx/XNfs6v4KK9OGMwVWR3dDskEmO7JrdlXWklZpcftUM6YT0lGREaLyAYRyReRBxt4P0pE5jrvLxaRbnXee8gp3yAio5pqU0RmOeVrRORlEYlwykVEnnbqrxaRgWfzwY1pScu2FzN22tdUemqZe9dwLuqR5HZIJgB1d1bgDuSFMptMMiISBjwLXAlkATeJSFa9apOAYlXtCTwFPO7smwWMA/oCo4HnRCSsiTZnAX2AfkAMMNkpvxLIdH6mAM+fyQc2pqV9vG4ft7y0iISYCN66+yLO7xzvdkgmQGU405g37w/cS2a+9GSGAPmqukVVq4A5wJh6dcYAM53t+cBlIiJO+RxVrVTVrUC+016jbarqe+oAlgBpdY7xmvPWIiBBRFLP8HMb0yLmLt3BlNeX0btjG+bffRFdEm2pfnPmuiW2RgS2BnNPBugM7KzzusApa7COqnqAEiDxFPs22aZzmew24IPTiAMRmSIiuSKSu3//fh8+njFnT1V5+uNNPPDmN1zcM4nZdw4jKc5WUjZnJzoijE7xMQF9r4wvSaahu8XUxzqnW17Xc8DnqvrFacSBqk5X1WxVzU5OTm5gF2OaV02t8t9/W8OTzjIxL43PprWtQ2aaSffk1gE9w8yX34QCIL3O6zRgdyN1CkQkHIgHDjaxb6NtisivgGTgrtOMw5hzqqK6hnvfWME/1+7j7pE9uH9Ub7xXio1pHt0SW7Ny5yFUNSD/bfnSk1kKZIpIhohE4h3Iz6lXJwcY72yPBT5xxlRygHHO7LMMvIP2S07VpohMBkYBN6lqbb1j3O7MMhsGlKjqnjP4zMY0i5Lyam6bsZgP1+3jV9dm8cDoPgH5JWD8W9fEWA5XeDhUXu12KGekyZ6MqnpEZCqwAAgDXlbVPBF5BMhV1RxgBvC6iOTj7cGMc/bNE5F5wFrAA9yjqjUADbXpHHIasB1Y6PzCvqWqjwDvAVfhnTxQDkxojhNgzJnYUVTOHa8uoeDgUZ656UKu6d/J7ZBMkEpv7508suNgOe1aR7oczenz6cKxqr6H90u+btnDdbYrgBsa2fdR4FFf2nTKG4zJ6Rnd40u8xrSkFTuKmTwzF0+t8vqkIQztnuh2SCaIdXVmKG4/WM4F6QkuR3P6bHTSmNPw/jd7+OnclXRsG80rEwbTIznO7ZBMkOvi9GR2Hix3OZIzY0nGGB+oKi99sZXfvb+OAekJvHR7Nok2RdmcA7GR4STFRbG9KDCnMVuSMaYJnppafv2PPP6yaAdX9UvhyR8NsCdZmnOqa2Is24usJ2NM0DlS6WHq7OV8umE/d/1bdx4Y1cceNGbOua7tY1m0pcjtMM6IrcJsTCP2llRww7SFfL7pAL/7QT8euvI8SzDGFentY9lTWkGlp8btUE6b9WSMacDa3aVMmrmU0qPVzBifzcjeHdwOyYSwromxqMLOg0fp2SGwJptYT8aYehbk7WXstK8B+Ou/X2QJxrju2DTmQJxhZj0ZYxyqynOfbeYPCzYwID2B6bcNokPbaLfDMub4DZmBOMPMkowxeNcge/DN1fxt5W7GDOjE4z/sbzPIjN9IjosiNjKMHQePuh3KabMkY0Je4eEK7np9GSt2HOIXo3rz45E9bA0y41dEhC7tY9lx0HoyxgSUNbtKmPJaLsXl1Uy7dRCjz09xOyRjGpTePjYgH8NsA/8mZH2wZg83TFsIwPy7h1uCMX4tvV0suw4dxbuMY+CwnowJOarKs5/m88Q/N3oH+G8fRIc2NsBv/FvndjGUV9VQXF5N+wBajdmSjAkpRyo93D9/Ne9+s4fvD+jEYzbAbwJE54QYAHYVH7UkY4w/2l50hCmvLWNT4WEeurIPUy7pbgP8JmCktXOSzKFy+qXFuxyN7yzJmJDw2YZC7n1jBSLCzIlD+E5mstshGXNajvVkCooDaxqzJRkT1FSV5//lvcGyd8c2TL8tmy7O3dPGBJKE2AhiI8PYdciSjDF+oe74yzVaZifVAAAXXElEQVT9U/n92P7ERto/eROYRIS0djEB15PxaQqziIwWkQ0iki8iDzbwfpSIzHXeXywi3eq895BTvkFERjXVpohMdcpURJLqlI8UkRIRWen8HH/8szH1bS86wvXPfc37a/bw0JV9eOamCy3BmIDXOSGGXQGWZJr8rRORMOBZ4AqgAFgqIjmqurZOtUlAsar2FJFxwOPAjSKSBYwD+gKdgI9EpJezT2NtfgW8A3zWQDhfqOo1Z/A5TQix8RcTrDq3i2H5jkNuh3FafOnJDAHyVXWLqlYBc4Ax9eqMAWY62/OBy8Q7bWcMMEdVK1V1K5DvtNdom6q6QlW3neXnMiGotlZ5+uNNTHh1KZ0SYvjH1IstwZig0jkhlpKj1ZRVetwOxWe+JJnOwM46rwucsgbrqKoHKAEST7GvL202ZLiIrBKR90Wkb0MVRGSKiOSKSO7+/ft9aNIEg+IjVUx4dSlPfriRMRd04q0fX2QD/CbodG534l6ZQOHLReqGbiSov65BY3UaK28ouTW1VsJyoKuqlonIVcDfgMxvNaI6HZgOkJ2dHVjrL5gzsnLnIe6ZtZz9hyv57ffP55ahXez+FxOUTkxjLqd3ShuXo/GNLz2ZAiC9zus0YHdjdUQkHIgHDp5iX1/aPImqlqpqmbP9HhBRd2KACT2qyusLt3GD84Cx+XcP59ZhXS3BmKCVfvyGzMDpyfiSZJYCmSKSISKReAfyc+rVyQHGO9tjgU/Uu4pbDjDOmX2WgbfnscTHNk8iIinOOA8iMsSJvciXD2mCz5FKDz+du5Jf/j2Pi3sm8e69F9M/LcHtsIxpUUlxUUSGtQquy2Wq6hGRqcACIAx4WVXzROQRIFdVc4AZwOsiko+3BzPO2TdPROYBawEPcI+q1oB3qnL9Np3ye4H7gRRgtYi8p6qT8Savu0XEAxwFxmmgLUdqmkV+4WHu/styNu8v4/99rxc/HtmTVq2s92KCX6tWQqeEaAoCqCcjwfw9nZ2drbm5uW6HYZrRP1bt5oE3VxMTEcbTN13IiJ52xdSEllteWsSRyhr+ds+IFjuGiCxT1ezmaMvuTjMBoaK6hv95Zy2zFu9gUNd2/PnmC0mNj3E7LGPOuc4JMXyyPnBmzlqSMX4vv7CMqbOXs37vYaZc0p1fjOpNRJg9b8+Eps4JsRwoq6TSU0NUuP8/psKSjPFbqsr8ZQU8/Pc8YiLDeGXCYC7t3cHtsIxxVWq89wF7haWVpLf3/3vBLMkYv1RW6eGXf1vD2yt2Max7e/5444WkxNvTK41JTfD+Huw+dNSSjDFnIm93CVNnr2B70RF+enkm//HdTMJs9pgxwImezN7SCpcj8Y0lGeM3VJXXFm7n0XfX0a51BLPvHMaw7oluh2WMX0lxJrzsKbEkY4zPDpVX8cCbq1mQt49LeyfzxA0XkBgX5XZYxviduKhw2kSHsydA7pWxJGNc9/XmA9w3dxUHyir576vPY+KIDLu50phTSI2Ptp6MMU2p8tTyf//cwPQvtpCR2Jq3fnyRLQ1jjA9S42NsTMaYU8kvPMxP5qwkb3cpNw/twn9ffZ49udIYH6XGR5O3u9TtMHxiv9XmnFJV/rJ4B4++u5bYyHCm3zaI7/VNcTssYwJKSnw0B8oqqfLUEhnu3zcmW5Ix58yBskoemL+aj9cXckmvZJ4Y258Obe3eF2NOVydnhtm+0gq/v1fGkow5Jz7dUMgv/rqK0goPv7o2i/HDu9ngvjFn6NiNyXtKLMmYEFde5eHx99czc+F2+qS0YdbkYQHzRD9j/FWnhGNJxv+nMVuSMS1m2faD/HzeKrYVlTNxRAb3j+5NdIT/L+hnjL8LpBsyLcmYZlfpqeGpDzcx/fPNpMbH8Madwxjew+7cN6a5HLshc68lGRNq1uwq4efzVrFh32FuGpLOf12dRVyU/TMzprl5b8i0y2UmRFTX1PLcp5t55pNNtG8dySt3DObSPrYsvzEtJSU+JiAul/k0wVpERovIBhHJF5EHG3g/SkTmOu8vFpFudd57yCnfICKjmmpTRKY6ZSoiSXXKRUSedt5bLSIDz/RDm+a1ad9hrn/ua576aCNX90/lnz+7xBKMMS2sU4AsLdNkT0ZEwoBngSuAAmCpiOSo6to61SYBxaraU0TGAY8DN4pIFjAO6At0Aj4SkV7OPo21+RXwDvBZvVCuBDKdn6HA885/jUtqapUZX27hiX9uJC4qnOduGchV/VLdDsuYkBAoN2T6crlsCJCvqlsARGQOMAaom2TGAL92tucDfxYRccrnqGolsFVE8p32aKxNVV3hlNWPYwzwmqoqsEhEEkQkVVX3nM4HNs0jv/AwD7z5Dcu2F3NFVkd+94N+JLexVZONOVc6xceg6v83ZPqSZDoDO+u8LuDbPYjjdVTVIyIlQKJTvqjevp2d7aba9CWOzsBJSUZEpgBTALp06dJEk+Z0VdfUMv3zLfzpo03ERoXx5I8u4AcXdm7ojwJjTAtKqfPwskBPMg19e6iPdRorb6hvV7/NM4kDVZ0OTAfIzs5uqk1zGvJ2l3D//NXk7S7lqn4p/Oa68633YoxLjj0hc7efP1fGlyRTAKTXeZ0G7G6kToGIhAPxwMEm9m2qzTOJw7SASk8Nz3ycz7R/bSYhNpJptw5k9Pk29mKMm471ZPb5+ZL/vowWLQUyRSRDRCLxDuTn1KuTA4x3tscCnzhjJznAOGf2WQbeQfslPrZZXw5wuzPLbBhQYuMxLW/Z9mKufvpL/vxpPtcN6MRH911iCcYYP9AmOoK4qHC/n2HWZE/GGWOZCiwAwoCXVTVPRB4BclU1B5gBvO4M7B/EmzRw6s3DO0nAA9yjqjXgnapcv02n/F7gfiAFWC0i76nqZOA94CogHygHJjTXSTDfVl7l4YkFG3nl662kto3mlQmDubS3TUs2xp+kxEf7/V3/4u1wBKfs7GzNzc11O4yA8+WmA/zn29+w42A5tw7rwgOj+9AmOsLtsIwx9dw2YzFllR7e/vGIZm1XRJapanZztGV3/JvjDpRV8ui763h7xS66JcYyZ8owhnW3NceM8VcpbaP5Mv+A22GckiUZg6ry19wCfvf+Oo5Uerj3uz358aU9bcVkY/xcSnw0hYcr8dTUEh7mnzdkWpIJcfmFZfzn29+wZOtBBndrx+9+0I/Mjva8F2MCQUp8NDW1yoGyquOzzfyNJZkQVVFdw3Ofbeb5z/KJiQjjsev78aPsdHtapTEBJDX+xMPLLMkYv/H15gP899tr2HLgCN8f0In/ujrLbqo0JgCltPU+vMyfZ5hZkgkhRWWV/O699by5vICuibG8PmkI38lMdjssY8wZOtGTsSRjXFRTq8xesoMnFmzgSKWHey7twX98N9MG9o0JcAmxEUSFt/Lru/4tyQS5FTuK+eXf17BmVynDuyfyyJi+NrBvTJAQEVL8/LkylmSCVFFZJb//YANzc3fSsW0Uz9x0Idf0T7XVko0JMilt/fuuf0syQab+pbG7LunOf1yWSVyU/a82JhilxkezbEex22E0yr55gohdGjMm9KTEx7CvZC+1teqXtyBYkgkCRWWV/GHBBuYs9V4ae/qmC7nWLo0ZExJS46OpqqnlYHkVSXH+dyuCJZkAVuWp5bWF2/jTx5s4WlXDlEu6c69dGjMmpBx/QmZJhSUZ0zxUlU/WF/Lou+vYcuAI/9YrmV9ecx49O9ilMWNCTUrbE0nm/M7xLkfzbZZkAsymfYf5n3fX8fnG/XRPbs0rdwzm0j72nBdjQtXxGzL99F4ZSzIB4lB5FX/8aBOvL9pObGQYv7wmi9uHdyXCT1deNcacG4lxUYS3EvaWHHU7lAZZkvFznppaZi3ewVMfbaT0aDU3D+3CfVf0pn3rSLdDM8b4gbBWQse2/ntDpiUZP6WqfLqhkP99bz2bCsu4qEciD1+bRZ+Utm6HZozxM/78GGafrrWIyGgR2SAi+SLyYAPvR4nIXOf9xSLSrc57DznlG0RkVFNtikiG08Ymp81Ip/wOEdkvIiudn8ln88H92fIdxdw4fRETX82luqaWF24bxKzJQy3BGGMalBIfzd5AHZMRkTDgWeAKoABYKiI5qrq2TrVJQLGq9hSRccDjwI0ikgWMA/oCnYCPRKSXs09jbT4OPKWqc0RkmtP2884+c1V16ll+Zr+VX1jGHxasZ0HePpLiovif75/PuMHpNu5ijDmllLbRfLq+EFX1u/vjfLlcNgTIV9UtACIyBxgD1E0yY4BfO9vzgT+L95OOAeaoaiWwVUTynfZoqE0RWQd8F7jZqTPTafdYkglK+0or+ONHG5mXW0B0eCvuu6IXky7OoLXd72KM8UFqfDTlVTWUVniIj4lwO5yT+PIt1hnYWed1ATC0sTqq6hGREiDRKV9Ub9/OznZDbSYCh1TV00B9gB+KyCXARuBnqlq3jYBTcrSaF/61mZe/2kpNrXL78K5MvbQniX54Q5Uxxn/VvSEzEJNMQ30v9bFOY+UNXf85VX2AfwBvqGqliPw73l7Od78VrMgUYApAly5dGmjOfRXVNfxl0Xb+/Gk+h8qr+f6ATtx3RW+6JMa6HZoxJgDVfQxz7xT/uinblyRTAKTXeZ0G7G6kToGIhAPxwMEm9m2o/ACQICLhTm/meH1VLapT/0W8YzffoqrTgekA2dnZ9ZOhq2pqlb+t2MWTH25k16GjXNIrmftH9fbLu3SNMYEjJd5/H8PsS5JZCmSKSAawC+9A/s316uQA44GFwFjgE1VVEckBZovIk3gH/jOBJXh7LN9q09nnU6eNOU6bfwcQkVRV3eMc7zpg3Rl+5nNOVflsw34e/2A96/cepl/neH4/tj8jeia5HZoxJgh0aBOFCH45w6zJJOOMsUwFFgBhwMuqmicijwC5qpoDzABedwb2D+JNGjj15uGdJOAB7lHVGoCG2nQO+QAwR0R+C6xw2ga4V0Suc9o5CNxx1p/+HFixo5jH3l/P4q0H6ZoYy59vvpCrzk/1yyW5jTGBKSKsFUlxUX7ZkxFVv7qi1Kyys7M1NzfXlWNv3l/GEws28P6avSTFRfKTyzK5cXAXIsNtOrIxpvld9+cvaRcbycyJQ5qu3AQRWaaq2c0Qlt3x39wKSyv408ebmLN0J9HhrfjZ5b2Y/B2bjmyMaVkpbaPZXlTudhjfYt98zeRIpYcXPt/Ci59vwVNby23DujL1uz398vkOxpjgkxofzaItRU1XPMcsyZwlT00tc3N38tSHmzhQVsnV/VO5f1Rvuia2djs0Y0wISYmPobTCw5FKj19dOfGfSALQws1FPPz3NWwqLGNwt3a8ePsgLuzSzu2wjDEh6Ni9MntLK+iRHOdyNCdYkjkDB49U8ei763hzeQHp7WOYdusgRvXt6HdrBhljQkdH5wmZ+0osyQS0r/MP8NO5Kykur+KeS3sw9dJMYiLD3A7LGBPiTtz171/TmC3J+EhVefGLLfzv++vpntSaVycMIauTLb1vjPEPKXUul/kTSzI+UFUee389L3y+hav7p/KHsf2JjbRTZ4zxH9ERYbSLjWCPnz2G2b4pffDSF1t54fMt3DasK7+5rq/drW+M8Usp8TF+d9e/3X7ehOU7ivnf99dxVb8USzDGGL+WGh/td2MylmROoaZWeWD+alLaRvP4D/tbgjHG+LWObaPZ52djMpZkTuGd1bvZVFjGL6/Jok20fz0IyBhj6kuNj+ZAWRWVnhq3QznOkswpzPhyK706xjGqb4rboRhjTJOOzTArLK10OZITLMk0Ysv+MlYXlPCj7HS7TGaMCQj+eK+MJZlGfLyuEICr+6e6HIkxxvim7mOY/YUlmUYs215Ml/axpDqPNTXGGH/nj49htiTTiFUFhxjYJcHtMIwxxmdxUeHERYX71V3/lmQaUFFdw56SCrr70SJzxhjji5T4aHYV2+Uyv1bg/A9Kb2+XyowxgaVnchz5hWVuh3GcT0lGREaLyAYRyReRBxt4P0pE5jrvLxaRbnXee8gp3yAio5pqU0QynDY2OW1GNnWM5nbwSBUAyXHRLXUIY4xpEb1T2rCt6AhHq/zjXpkmk4yIhAHPAlcCWcBNIpJVr9okoFhVewJPAY87+2YB44C+wGjgOREJa6LNx4GnVDUTKHbabvQYLaG8ygNgS/gbYwLOealtqFXYVHjY7VAA33oyQ4B8Vd2iqlXAHGBMvTpjgJnO9nzgMvE+wWsMMEdVK1V1K5DvtNdgm84+33XawGnz+00co9lVVHv/Aoi1JGOMCTD907wTlr7YdMDlSLx8WYW5M7CzzusCYGhjdVTVIyIlQKJTvqjevp2d7YbaTAQOqaqngfqNHeOkMykiU4ApzssyESmqX8dXWS3WV3JNEmd4LoKQnQsvOw8nBNW5mPo4TD2zXZOArs0Vhy9JpqHegvpYp7HyhnpQp6rvaxyo6nRg+vHARHJVNbuBfUOOnYsT7Fx42Xk4wc6Fl3MeujVXe75cLisA0uu8TgN2N1ZHRMKBeODgKfZtrPwAkOC0Uf9YjR3DGGOMn/IlySwFMp1ZX5F4B/Jz6tXJAcY722OBT1RVnfJxzsywDCATWNJYm84+nzpt4LT59yaOYYwxxk81ebnMGf+YCiwAwoCXVTVPRB4BclU1B5gBvC4i+Xh7F+OcffNEZB6wFvAA96hqDUBDbTqHfACYIyK/BVY4bdPYMXwwvekqIcPOxQl2LrzsPJxg58KrWc+DWGfAGGNMS7E7/o0xxrQYSzLGGGNaTFAnmaaWwwl0IvKyiBSKyJo6Ze1F5ENnWZ4PRaSdUy4i8rRzLlaLyMA6+4x36m8SkfENHcvfiUi6iHwqIutEJE9EfuKUh9T5EJFoEVkiIquc8/Abp/y0l2tqbEmoQOOsMrJCRN5xXofkuRCRbSLyjYisFJFcp6zlfz9UNSh/8E4o2Ax0ByKBVUCW23E182e8BBgIrKlT9nvgQWf7QeBxZ/sq4H289xsNAxY75e2BLc5/2znb7dz+bGdwLlKBgc52G2Aj3iWLQup8OJ8nztmOABY7n28eMM4pnwbc7Wz/GJjmbI8D5jrbWc7vTBSQ4fwuhbn9+c7wnNwHzAbecV6H5LkAtgFJ9cpa/PcjmHsyviyHE9BU9XO+fa9Q3eV36i/L85p6LcJ7P1IqMAr4UFUPqmox8CHedeYCiqruUdXlzvZhYB3eVSJC6nw4n+fYErwRzo9y+ss1NbYkVEARkTTgauAl5/WZLF0VFOeiES3++xHMSaah5XA6N1I3mHRU1T3g/eIFOjjljZ2PoDtPzmWOC/H+FR9y58O5PLQSKMT7JbAZH5drAuouCRXQ58HxR+B+oNZ57fPSVQTfuVDgnyKyTLzLb8E5+P3wZVmZQOXTMjQh5HSX/glIIhIHvAn8VFVLpfE1VIP2fKj3XrQBIpIAvA2c11A1579Bex5E5BqgUFWXicjIY8UNVA36c+EYoaq7RaQD8KGIrD9F3WY7F8Hck/FlOZxgtM/p1uL8t9ApP90lfgKOiETgTTCzVPUtpzhkz4eqHgI+w3tN/XSXawqG8zACuE5EtuG9XP5dvD2bUDwXqOpu57+FeP/4GMI5+P0I5iTjy3I4waju8jv1l+W53Zk1MgwocbrHC4DviUg7Z2bJ95yygOJcO58BrFPVJ+u8FVLnQ0SSnR4MIhIDXI53fOp0l2tqbEmogKGqD6lqmnoXexyH97PdQgieCxFpLSJtjm3j/Xe9hnPx++H2jIeW/ME7Q2Ij3mvS/+V2PC3w+d4A9gDVeP/CmIT3GvLHwCbnv+2duoL3QXGbgW+A7DrtTMQ7mJkPTHD7c53hubgYb7d9NbDS+bkq1M4H0B/vckyrnS+Rh53y7ni/GPOBvwJRTnm08zrfeb97nbb+yzk/G4Ar3f5sZ3leRnJidlnInQvnM69yfvKOfR+ei98PW1bGGGNMiwnmy2XGGGNcZknGGGNMi7EkY4wxpsVYkjHGGNNiLMkYY4xpMZZkjDHGtBhLMsYYY1rM/we5gZknDF8RXwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] @@ -554,31 +550,37 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 10, "metadata": {}, "outputs": [], "source": [ - "# total_f.update_integration_options(draws_per_dim=10000000, mc_sampler=None)\n", - "# inte = total_f.integrate(limits = (3000, 3200), norm_range=False)\n", - "# print(zfit.run(inte))\n", - "# print(pdg[\"jpsi_BR\"]/pdg[\"NR_BR\"], zfit.run(inte)/pdg[\"NR_auc\"])" + "# total_f.update_integration_options(draws_per_dim=2000000, mc_sampler=None)\n", + "# inte = total_f.integrate(limits = (3090, 3102), norm_range=False)\n", + "# inte_fl = zfit.run(inte)\n", + "# print(inte_fl)\n", + "# print(pdg[\"jpsi_BR\"]/pdg[\"NR_BR\"], inte_fl/pdg[\"NR_auc\"])" ] }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 11, "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": 13, + "execution_count": 12, "metadata": {}, "outputs": [], "source": [ @@ -610,7 +612,7 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 13, "metadata": {}, "outputs": [ { @@ -627,7 +629,7 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 14, "metadata": {}, "outputs": [], "source": [ @@ -648,58 +650,53 @@ }, { "cell_type": "code", - "execution_count": 16, - "metadata": {}, - "outputs": [], - "source": [ - "nevents = int(pdg[\"number_of_decays\"])\n", - "event_stack = 100000\n", - "\n", - "calls = int(nevents/event_stack + 1)\n", - "\n", - "# total_samp = []\n", - "\n", - "# start = time.time()\n", - "\n", - "# for call in range(calls):\n", - "# samp = total_f.sample(n=event_stack)\n", - "# sam = samp.unstack_x()\n", - "# sam = zfit.run(sam)\n", - "# clear_output(wait=True)\n", - " \n", - "# print(\"{0}/{1}\".format(call + 1, calls))\n", - "# print(\"Time taken: {}\".format(display_time(int(time.time() - start))))\n", - "# c = call + 1\n", - "# print(\"Projected time left: {}\".format(display_time(int((time.time() - start)/c*(calls-c)))))\n", - " \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": 17, - "metadata": {}, - "outputs": [], - "source": [ - "# print(\"Time to generate full toy: {} s\".format(int(time.time()-start)))" - ] - }, - { - "cell_type": "code", - "execution_count": 18, + "execution_count": null, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "(5404696,)\n" + "1/1081\n", + "Time taken: 5 min, 13 s\n", + "Projected time left: 3 d, 21 h\n" ] } ], "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", + "\n", + "for call in range(calls):\n", + " sam = samp.unstack_x()\n", + " sam = zfit.run(sam)\n", + " clear_output(wait=True)\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", @@ -716,22 +713,9 @@ }, { "cell_type": "code", - "execution_count": 19, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAD8CAYAAABDwhLXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAFslJREFUeJzt3X+MXeWd3/H3p3YgJFliAyalNls7jdUuQduGjAhtqlUUtmAgivmDSEZRsbJIVlPSZptWiWmkoiZdibTVskVKWKHgxlRpCGWzwgoQr0WIokqBMOQHP0KIJ4TCLGzsxMCyjZos2W//uM+Ey+TODJ4Znrn2vF/S1T3ne55znnMe+frjc+7DJVWFJEmvtr+10icgSVodDBxJUhcGjiSpCwNHktSFgSNJ6sLAkSR1YeBIkrowcCRJXRg4kqQu1q70CfRy2mmn1ebNm1f6NCTpmPLAAw/8pKo2LMexVk3gbN68mcnJyZU+DUk6piT5P8t1LB+pSZK6MHAkSV0YOJKkLgwcSVIXBo4kqQsDR5LUhYEjSepiwcBJsifJoSQPj9j275JUktPaepJcn2QqyYNJzhlquzPJwfbaOVR/e5KH2j7XJ0mrn5LkQGt/IMn6hfqQJI2vV3KH8zlg2+xikjOBfwY8OVS+CNjaXruAG1rbU4BrgHcA5wLXzARIa7NraL+ZvnYDd1fVVuDutj5nH5Kk8bZg4FTV14EjIzZdB3wUqKHaduDmGrgXWJfkDOBC4EBVHamqZ4EDwLa27eSq+kZVFXAzcOnQsfa25b2z6qP6kPQq2Lz7jpU+BR0nFvUdTpL3An9eVd+dtWkj8NTQ+nSrzVefHlEHeFNVPQPQ3k9foA9J0hg76t9SS/I64OPABaM2j6jVIurznsIr3SfJLgaP3fjN3/zNBQ4rSXo1LeYO5+8BW4DvJnkC2AR8K8nfZnC3ceZQ203A0wvUN42oA/x45lFZez/U6nMd69dU1Y1VNVFVExs2LMuPnUqSFumoA6eqHqqq06tqc1VtZhAA51TVXwD7gCvaTLLzgOfb47D9wAVJ1rfJAhcA+9u2F5Kc12anXQHc3rraB8zMZts5qz6qD0nSGFvwkVqSLwDvAk5LMg1cU1U3zdH8TuBiYAr4GfABgKo6kuSTwP2t3SeqamYiwgcZzIQ7CbirvQCuBW5NciWDmXDvm68PSdJ4WzBwquryBbZvHlou4Ko52u0B9oyoTwJnj6j/FDh/RH3OPiRJ48tfGpAkdWHgSJK6MHAkSV0YOJKkLgwcSVIXBo4kqQsDR5LUhYEjSerCwJEkdWHgSJK6MHAkSV0YOJKkLgwcSVIXBo4kqQsDR5LUhYEjSerCwJEkdWHgSJK6MHAkSV0YOJKkLhYMnCR7khxK8vBQ7b8k+X6SB5P8aZJ1Q9uuTjKV5LEkFw7Vt7XaVJLdQ/UtSe5LcjDJF5Oc0OontvWptn3zQn1IksbXK7nD+RywbVbtAHB2Vf028APgaoAkZwE7gLe2fT6TZE2SNcCngYuAs4DLW1uATwHXVdVW4Fngyla/Eni2qt4CXNfazdnHUV63JKmzBQOnqr4OHJlV+7OqerGt3gtsasvbgVuq6udV9SNgCji3vaaq6vGq+gVwC7A9SYB3A7e1/fcClw4da29bvg04v7Wfqw9J0hhbju9wfg+4qy1vBJ4a2jbdanPVTwWeGwqvmfrLjtW2P9/az3UsSdIYW1LgJPk48CLw+ZnSiGa1iPpijjXq/HYlmUwyefjw4VFNJEmdLDpwkuwE3gO8v6pm/sKfBs4carYJeHqe+k+AdUnWzqq/7Fht+xsZPNqb61i/pqpurKqJqprYsGHDYi5TkrRMFhU4SbYBHwPeW1U/G9q0D9jRZphtAbYC3wTuB7a2GWknMPjSf18LqnuAy9r+O4Hbh461sy1fBny1tZ+rD0nSGFu7UIMkXwDeBZyWZBq4hsGstBOBA4Pv8bm3qv5FVT2S5FbgewwetV1VVb9sx/kQsB9YA+ypqkdaFx8Dbknyn4BvAze1+k3A/0gyxeDOZgfAfH1IksZXXnoadnybmJioycnJlT4N6ZizefcdPHHtJSt9GlohSR6oqonlOJa/NCBJ6sLAkSR1YeBIkrowcCRJXRg4kqQuDBxJUhcGjiSpCwNHktSFgSNJ6sLAkSR1YeBIkrowcCRJXRg4kqQuDBxJUhcGjiSpCwNHktSFgSNJ6sLAkSR1YeBIkrowcCRJXRg4kqQuFgycJHuSHEry8FDtlCQHkhxs7+tbPUmuTzKV5MEk5wzts7O1P5hk51D97UkeavtcnySL7UOSNL5eyR3O54Bts2q7gburaitwd1sHuAjY2l67gBtgEB7ANcA7gHOBa2YCpLXZNbTftsX0IUkabwsGTlV9HTgyq7wd2NuW9wKXDtVvroF7gXVJzgAuBA5U1ZGqehY4AGxr206uqm9UVQE3zzrW0fQhSRpji/0O501V9QxAez+91TcCTw21m261+erTI+qL6ePXJNmVZDLJ5OHDh4/qAiVJy2u5Jw1kRK0WUV9MH79erLqxqiaqamLDhg0LHFaS9GpabOD8eOYxVns/1OrTwJlD7TYBTy9Q3zSivpg+JEljbLGBsw+YmWm2E7h9qH5Fm0l2HvB8exy2H7ggyfo2WeACYH/b9kKS89rstCtmHeto+pAkjbG1CzVI8gXgXcBpSaYZzDa7Frg1yZXAk8D7WvM7gYuBKeBnwAcAqupIkk8C97d2n6iqmYkIH2QwE+4k4K724mj7kCSNtwUDp6oun2PT+SPaFnDVHMfZA+wZUZ8Ezh5R/+nR9iFJGl/+0oAkqQsDR5LUhYEjSerCwJEkdWHgSJK6MHAkSV0YOJKkLgwcSVIXBo4kqQsDR5LUhYEjSerCwJEkdWHgSJK6MHAkSV0YOJKkLgwcSVIXBo4kqQsDR5LUhYEjSerCwJEkdbGkwEnyb5I8kuThJF9I8tokW5Lcl+Rgki8mOaG1PbGtT7Xtm4eOc3WrP5bkwqH6tlabSrJ7qD6yD0nS+Fp04CTZCPxrYKKqzgbWADuATwHXVdVW4FngyrbLlcCzVfUW4LrWjiRntf3eCmwDPpNkTZI1wKeBi4CzgMtbW+bpQ5I0ppb6SG0tcFKStcDrgGeAdwO3te17gUvb8va2Ttt+fpK0+i1V9fOq+hEwBZzbXlNV9XhV/QK4Bdje9pmrD0nSmFp04FTVnwP/FXiSQdA8DzwAPFdVL7Zm08DGtrwReKrt+2Jrf+pwfdY+c9VPnacPSdKYWsojtfUM7k62AH8HeD2Dx1+z1cwuc2xbrvqoc9yVZDLJ5OHDh0c1kSR1spRHar8L/KiqDlfVXwNfAv4JsK49YgPYBDzdlqeBMwHa9jcCR4brs/aZq/6Tefp4maq6saomqmpiw4YNS7hUSdJSLSVwngTOS/K69r3K+cD3gHuAy1qbncDtbXlfW6dt/2pVVavvaLPYtgBbgW8C9wNb24y0ExhMLNjX9pmrD0nSmFrKdzj3Mfji/lvAQ+1YNwIfAz6SZIrB9y03tV1uAk5t9Y8Au9txHgFuZRBWXwGuqqpftu9oPgTsBx4Fbm1tmacPSdKYyuCG4fg3MTFRk5OTK30a0jFn8+47eOLaS1b6NLRCkjxQVRPLcSx/aUCS1IWBI0nqwsCRJHVh4EiSujBwJEldGDiSpC4MHElSFwaOJKkLA0eS1IWBI0nqwsCRJHVh4Eg6Zm3efcdKn4KOgoEjSerCwJEkdWHgSJK6MHAkSV0YOJKkLgwcSVIXBo4kqQsDR5LUxZICJ8m6JLcl+X6SR5P84ySnJDmQ5GB7X9/aJsn1SaaSPJjknKHj7GztDybZOVR/e5KH2j7XJ0mrj+xDkjS+lnqH89+Ar1TVPwD+IfAosBu4u6q2Ane3dYCLgK3ttQu4AQbhAVwDvAM4F7hmKEBuaG1n9tvW6nP1IUkaU4sOnCQnA78D3ARQVb+oqueA7cDe1mwvcGlb3g7cXAP3AuuSnAFcCByoqiNV9SxwANjWtp1cVd+oqgJunnWsUX1IksbUUu5w3gwcBv57km8n+WyS1wNvqqpnANr76a39RuCpof2nW22++vSIOvP0IUkaU0sJnLXAOcANVfU24P8y/6OtjKjVIuqvWJJdSSaTTB4+fPhodpUkLbOlBM40MF1V97X12xgE0I/b4zDa+6Gh9mcO7b8JeHqB+qYRdebp42Wq6saqmqiqiQ0bNizqIiVJy2PRgVNVfwE8leTvt9L5wPeAfcDMTLOdwO1teR9wRZutdh7wfHscth+4IMn6NlngAmB/2/ZCkvPa7LQrZh1rVB+SpDG1don7/yvg80lOAB4HPsAgxG5NciXwJPC+1vZO4GJgCvhZa0tVHUnySeD+1u4TVXWkLX8Q+BxwEnBXewFcO0cfkqQxtaTAqarvABMjNp0/om0BV81xnD3AnhH1SeDsEfWfjupDkjS+/KUBSVIXBo4kqQsDR5LUhYEjSerCwJEkdWHgSJK6MHAkSV0YOJKkLgwcSVIXBo4kqQsDR5LUhYEjSerCwJEkdWHgSJK6MHAkSV0YOJKkLgwcSVIXBo4kqQsDR5LUhYEjSerCwJEkdbHkwEmyJsm3k3y5rW9Jcl+Sg0m+mOSEVj+xrU+17ZuHjnF1qz+W5MKh+rZWm0qye6g+sg9J0vhajjucDwOPDq1/CriuqrYCzwJXtvqVwLNV9RbgutaOJGcBO4C3AtuAz7QQWwN8GrgIOAu4vLWdrw9J0phaUuAk2QRcAny2rQd4N3Bba7IXuLQtb2/rtO3nt/bbgVuq6udV9SNgCji3vaaq6vGq+gVwC7B9gT4kSWNqqXc4fwR8FPibtn4q8FxVvdjWp4GNbXkj8BRA2/58a/+r+qx95qrP18fLJNmVZDLJ5OHDhxd7jZKkZbDowEnyHuBQVT0wXB7RtBbYtlz1Xy9W3VhVE1U1sWHDhlFNJEmdrF3Cvu8E3pvkYuC1wMkM7njWJVnb7kA2AU+39tPAmcB0krXAG4EjQ/UZw/uMqv9knj4kSWNq0Xc4VXV1VW2qqs0MvvT/alW9H7gHuKw12wnc3pb3tXXa9q9WVbX6jjaLbQuwFfgmcD+wtc1IO6H1sa/tM1cfkqQx9Wr8dzgfAz6SZIrB9y03tfpNwKmt/hFgN0BVPQLcCnwP+ApwVVX9st29fAjYz2AW3K2t7Xx9SJLG1FIeqf1KVX0N+FpbfpzBDLPZbf4f8L459v8D4A9G1O8E7hxRH9mHJGl8+UsDkqQuDBxJUhcGjiSpCwNHktSFgSNJ6sLAkSR1YeBIkrowcCRJXRg4kqQuDBxJUhcGjiSpCwNHktSFgSNJ6sLAkSR1YeBIkrowcCRJXRg4kqQuDBxJUhcGjiSpCwNHktTFogMnyZlJ7knyaJJHkny41U9JciDJwfa+vtWT5PokU0keTHLO0LF2tvYHk+wcqr89yUNtn+uTZL4+JEnjayl3OC8C/7aqfgs4D7gqyVnAbuDuqtoK3N3WAS4CtrbXLuAGGIQHcA3wDuBc4JqhALmhtZ3Zb1urz9WHJGlMLTpwquqZqvpWW34BeBTYCGwH9rZme4FL2/J24OYauBdYl+QM4ELgQFUdqapngQPAtrbt5Kr6RlUVcPOsY43qQ5I0ppblO5wkm4G3AfcBb6qqZ2AQSsDprdlG4Kmh3aZbbb769Ig68/Qx+7x2JZlMMnn48OHFXp4kaRksOXCSvAH4E+D3q+ov52s6olaLqL9iVXVjVU1U1cSGDRuOZldJ0jJbUuAkeQ2DsPl8VX2plX/cHofR3g+1+jRw5tDum4CnF6hvGlGfrw9J0phayiy1ADcBj1bVHw5t2gfMzDTbCdw+VL+izVY7D3i+PQ7bD1yQZH2bLHABsL9teyHJea2vK2Yda1QfkqQxtXYJ+74T+OfAQ0m+02r/HrgWuDXJlcCTwPvatjuBi4Ep4GfABwCq6kiSTwL3t3afqKojbfmDwOeAk4C72ot5+pAkjalFB05V/W9Gf88CcP6I9gVcNcex9gB7RtQngbNH1H86qg9J0vjylwYkSV0YOJKkLgwcSVIXBo4kqQsDR5LUhYEjSerCwJEkdWHgSJK6MHAkSV0YOJKkLgwcSVIXBo4kqQsDR5LUhYEjSerCwJEkdWHgSJK6MHAkSV0YOJKkLgwcSVIXBo4kqYtjOnCSbEvyWJKpJLtX+nwkSXM7ZgMnyRrg08BFwFnA5UnOWtmzkiTNZe1Kn8ASnAtMVdXjAEluAbYD31vRs5K0JJt338ET117C5t13vOL24+6Jay9Z6VMYC8dy4GwEnhpanwbesULncsw6Fj6sWnm9/5wcb38ue17POIfbsRw4GVGrlzVIdgG72upfJfkp8JNX+8SOEafhWMxwLAYch5ccs2ORTy3r4U4D/u5yHexYDpxp4Myh9U3A08MNqupG4MaZ9SSTVTXR5/TGm2PxEsdiwHF4iWMx0MZh83Id75idNADcD2xNsiXJCcAOYN8Kn5MkaQ7H7B1OVb2Y5EPAfmANsKeqHlnh05IkzeGYDRyAqroTuPModrlx4SarhmPxEsdiwHF4iWMxsKzjkKpauJUkSUt0LH+HI0k6hqyawDnefwYnyZ4kh5I8PFQ7JcmBJAfb+/pWT5Lr21g8mOScoX12tvYHk+xciWtZqiRnJrknyaNJHkny4VZfVeOR5LVJvpnku20c/mOrb0lyX7umL7ZJNyQ5sa1Pte2bh451das/luTClbmipUuyJsm3k3y5ra/KsUjyRJKHknwnyWSrvfqfj6o67l8MJhX8EHgzcALwXeCslT6vZb7G3wHOAR4eqv1nYHdb3g18qi1fDNzF4L9lOg+4r9VPAR5v7+vb8vqVvrZFjMUZwDlt+TeAHzD4+aNVNR7tet7Qll8D3Neu71ZgR6v/MfDBtvwvgT9uyzuAL7bls9pn5kRgS/ssrVnp61vkmHwE+J/Al9v6qhwL4AngtFm1V/3zsVrucH71MzhV9Qtg5mdwjhtV9XXgyKzydmBvW94LXDpUv7kG7gXWJTkDuBA4UFVHqupZ4ACw7dU/++VVVc9U1bfa8gvAowx+mWJVjUe7nr9qq69prwLeDdzW6rPHYWZ8bgPOT5JWv6Wqfl5VPwKmGHymjilJNgGXAJ9t62GVjsUcXvXPx2oJnFE/g7Nxhc6lpzdV1TMw+EsYOL3V5xqP426c2qOQtzH41/2qG4/2COk7wCEGfyH8EHiuql5sTYav6VfX27Y/D5zKcTAOzR8BHwX+pq2fyuodiwL+LMkDGfwiC3T4fBzT06KPwoI/g7PKzDUex9U4JXkD8CfA71fVXw7+gTq66YjacTEeVfVL4B8lWQf8KfBbo5q19+N2HJK8BzhUVQ8keddMeUTT434smndW1dNJTgcOJPn+PG2XbSxWyx3Ogj+Dc5z6cbv1pb0favW5xuO4Gackr2EQNp+vqi+18qodj6p6Dvgag2fw65LM/GNz+Jp+db1t+xsZPKY9HsbhncB7kzzB4JH6uxnc8azGsaCqnm7vhxj8Q+RcOnw+VkvgrNafwdkHzMwc2QncPlS/os0+OQ94vt1C7wcuSLK+zVC5oNWOKe1Z+03Ao1X1h0ObVtV4JNnQ7mxIchLwuwy+z7oHuKw1mz0OM+NzGfDVGnw7vA/Y0WZubQG2At/scxXLo6qurqpNNfhdsB0Mru39rMKxSPL6JL8xs8zgz/XD9Ph8rPRsiV4vBjMtfsDgGfbHV/p8XoXr+wLwDPDXDP7lcSWDZ853Awfb+ymtbRj8z+t+CDwETAwd5/cYfBE6BXxgpa9rkWPxTxnc2j8IfKe9Ll5t4wH8NvDtNg4PA/+h1d/M4C/JKeB/ASe2+mvb+lTb/uahY328jc9jwEUrfW1LHJd38dIstVU3Fu2av9tej8z8fdjj8+EvDUiSulgtj9QkSSvMwJEkdWHgSJK6MHAkSV0YOJKkLgwcSVIXBo4kqQsDR5LUxf8Htmuht2g0nfcAAAAASUVORK5CYII=\n", - "text/plain": [ - "
" - ] - }, - "metadata": { - "needs_background": "light" - }, - "output_type": "display_data" - } - ], + "outputs": [], "source": [ "bins = int((x_max-x_min)/7)\n", "\n", @@ -741,10 +725,10 @@ "\n", "# plt.plot(sam, calcs, '.')\n", "# plt.plot(test_q, calcs_test)\n", - "# plt.ylim(6000, 10000)\n", - "# plt.xlim(3000, 3750)\n", + "plt.ylim(4000, 12000)\n", + "plt.xlim(3000, 3750)\n", "\n", - "plt.savefig('test.png')" + "plt.savefig('test2.png')" ] }, { @@ -756,7 +740,7 @@ }, { "cell_type": "code", - "execution_count": 20, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -786,316 +770,9 @@ }, { "cell_type": "code", - "execution_count": 21, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": { - "text/html": [ - "
" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/html": [ - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
FCN = 19224270.58100143TOTAL NCALL = 58NCALLS = 58
EDM = 1.5229666301917304e-07GOAL EDM = 5e-06\n", - " UP = 0.5
\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
ValidValid ParamAccurate CovarPosDefMade PosDef
TrueTrueTrueTrueFalse
Hesse FailHasCovAbove EDMReach calllim
FalseTrueFalseFalse
" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/html": [ - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
+NameValueHesse ErrorMinos Error-Minos Error+Limit-Limit+Fixed?
0psi2s_p-9.424780.00651368No
1jpsi_s444.4950.220516No
2jpsi_p-9.424780.00631191No
3psi2s_s74.73490.048852No
\n", - "
\n",
-       "\n",
-       "
" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/html": [ - "
" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/html": [ - "Minos status for psi2s_p: VALID\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
Error-0.0065151409160649580.006517797812598556
ValidTrueTrue
At LimitFalseFalse
Max FCNFalseFalse
New MinFalseFalse
" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/html": [ - "Minos status for jpsi_s: VALID\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
Error-0.220659286900404150.22053488031823795
ValidTrueTrue
At LimitFalseFalse
Max FCNFalseFalse
New MinFalseFalse
" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/html": [ - "Minos status for jpsi_p: VALID\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
Error-0.0063110017718012870.006310920098439091
ValidTrueTrue
At LimitFalseFalse
Max FCNFalseFalse
New MinFalseFalse
" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/html": [ - "Minos status for psi2s_s: VALID\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
Error-0.049227653394393410.04926101726478015
ValidTrueTrue
At LimitFalseFalse
Max FCNFalseFalse
New MinFalseFalse
" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "psi2s_p: ^{+0.006517797812598556}_{-0.006515140916064958}\n", - "jpsi_s: ^{+0.22053488031823795}_{-0.22065928690040415}\n", - "jpsi_p: ^{+0.006310920098439091}_{-0.006311001771801287}\n", - "psi2s_s: ^{+0.04926101726478015}_{-0.04922765339439341}\n", - "Function minimum: 19224270.58100143\n" - ] - } - ], + "outputs": [], "source": [ "nll = zfit.loss.UnbinnedNLL(model=total_f, data=data2, fit_range = (x_min, x_max))\n", "\n", @@ -1113,22 +790,11 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "-3.1420346928204133" - ] - }, - "execution_count": 22, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ - "(-9.42522+2*np.pi)#/np.pi" + "(-3.14+2*np.pi)/np.pi" ] }, { @@ -1140,37 +806,18 @@ }, { "cell_type": "code", - "execution_count": 23, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "'5 h, 55 min'" - ] - }, - "execution_count": 23, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "display_time(int(395*pdg[\"number_of_decays\"]/100000))" ] }, { "cell_type": "code", - "execution_count": 24, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "6 h, 12 min\n" - ] - } - ], + "outputs": [], "source": [ "print(display_time(22376))" ] @@ -1180,6 +827,36 @@ "execution_count": null, "metadata": {}, "outputs": [], + "source": [ + "probs = total_f.pdf(test_q)\n", + "\n", + "calcs_test = zfit.run(probs)\n", + "res_y = zfit.run(jpsi_res(test_q))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.clf()\n", + "# plt.plot(x_part, calcs, '.')\n", + "plt.plot(test_q, calcs_test, label = 'pdf')\n", + "# plt.plot(test_q, res_y, label = 'res')\n", + "plt.legend()\n", + "plt.ylim(0.0, 4e-4)\n", + "# plt.yscale('log')\n", + "# plt.xlim(3080, 3110)\n", + "plt.savefig('test3.png')\n", + "print(jpsi_width)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], "source": [] } ], diff --git a/test.png b/test.png index 5829e42..4d21890 100644 --- a/test.png +++ b/test.png Binary files differ diff --git a/test2.png b/test2.png new file mode 100644 index 0000000..e9c5be6 --- /dev/null +++ b/test2.png Binary files differ diff --git a/test3.png b/test3.png new file mode 100644 index 0000000..a86e763 --- /dev/null +++ b/test3.png Binary files differ