From 96d3003311ff72e0f07b45eb36eca73a317f66ab Mon Sep 17 00:00:00 2001 From: mmarchegiani <matmarcheg@gmail.com> Date: Mon, 16 Jan 2023 17:00:45 +0100 Subject: [PATCH 01/11] Update notebook --- notebooks/jec.ipynb | 182 ++++++++++++++++++++++++++++++-------------- 1 file changed, 124 insertions(+), 58 deletions(-) diff --git a/notebooks/jec.ipynb b/notebooks/jec.ipynb index 28fb1e17..2993d7a9 100644 --- a/notebooks/jec.ipynb +++ b/notebooks/jec.ipynb @@ -10,7 +10,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "/work/dvalsecc/miniconda3/envs/pocket-coffea/lib/python3.9/site-packages/coffea/util.py:154: FutureWarning: In coffea version v0.8.0 (target date: 31 Dec 2022), this will be an error.\n", + "/work/mmarcheg/miniconda3/envs/pocket-coffea/lib/python3.8/site-packages/coffea/util.py:154: FutureWarning: In coffea version v0.8.0 (target date: 31 Dec 2022), this will be an error.\n", "(Set coffea.deprecations_as_errors = True to get a stack trace now.)\n", "ImportError: coffea.hist is deprecated\n", " warnings.warn(message, FutureWarning)\n" @@ -27,7 +27,7 @@ }, { "cell_type": "code", - "execution_count": 26, + "execution_count": 2, "id": "b0fbce2b", "metadata": {}, "outputs": [], @@ -39,17 +39,17 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 3, "id": "fc676d4c", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "<Array [[0.0635, 0.04, ... 0.336, 0.158]] type='100 * var * float32[parameters={...'>" + "<Array [[0.0635, 0.04, ... -0.00543, -0.664]] type='400000 * var * float32[param...'>" ] }, - "execution_count": 4, + "execution_count": 3, "metadata": {}, "output_type": "execute_result" } @@ -60,17 +60,17 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 4, "id": "7f52a8fb", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "<Array [0.937, 0.96, 0.958, ... 0.664, 0.842] type='1098 * float32'>" + "<Array [0.937, 0.96, 0.958, ... 1.01, 1.66] type='4343693 * float32'>" ] }, - "execution_count": 5, + "execution_count": 4, "metadata": {}, "output_type": "execute_result" } @@ -82,20 +82,18 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 5, "id": "3628b6ad", "metadata": {}, "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAARK0lEQVR4nO3dbYxcV33H8e+POAEEFOdh61q2i6mwQGlVglmlRlSIYlHlocKRClFQRUzkylWbtiAqFZcXRVR9Ed5ASVsFWYTWqSiQBmjcEGgtJwj1RQIbCAESaJaIyLaSeAnEPKSAQv99sccwWdae2d3ZXe/h+5FGc+655+78T67z890zd8apKiRJfXnGahcgSRo/w12SOmS4S1KHDHdJ6pDhLkkdWrfaBQBccMEFtXXr1tUuQ5LWlHvuuedbVTUx374zIty3bt3K1NTUapchSWtKkodPtc9lGUnqkOEuSR0y3CWpQ4a7JHXIcJekDhnuktQhw12SOjQ03JO8OMm9A4/vJnlrkvOSHEryYHs+t41PkuuTTCe5L8n25Z+GJGnQ0HCvqq9X1UVVdRHwcuBJ4BPAPuBwVW0DDrdtgEuBbe2xF7hhGeqWJJ3GQj+huhP4RlU9nGQX8OrWfwD4DPB2YBdwU83+KyB3JVmfZGNVPTKmmrWKtu775KKP/eZ1l4+xEkmns9A196uAD7f2hoHAfhTY0NqbgCMDxxxtfU+TZG+SqSRTMzMzCyxDknQ6I4d7knOA1wH/Nndfu0pf0L/XV1X7q2qyqiYnJub93htJ0iIt5Mr9UuALVfVY234syUaA9ny89R8Dtgwct7n1SZJWyELC/Y38bEkG4CCwu7V3A7cO9F/d7prZAZxwvV2SVtZIb6gmeQ7wWuCPBrqvA25Osgd4GLiy9d8OXAZMM3tnzTVjq1aSNJKRwr2qfgCcP6fvcWbvnpk7toBrx1KdJGlR/ISqJHXIcJekDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6NFK4J1mf5JYkX0vyQJJXJDkvyaEkD7bnc9vYJLk+yXSS+5JsX94pSJLmGvXK/X3Ap6vqJcBLgQeAfcDhqtoGHG7bAJcC29pjL3DDWCuWJA01NNyTPB94FXAjQFX9uKqeAHYBB9qwA8AVrb0LuKlm3QWsT7JxzHVLkk5jlCv3FwIzwD8l+WKSDyR5DrChqh5pYx4FNrT2JuDIwPFHW9/TJNmbZCrJ1MzMzOJnIEn6OaOE+zpgO3BDVb0M+AE/W4IBoKoKqIW8cFXtr6rJqpqcmJhYyKGSpCFGCfejwNGqurtt38Js2D92crmlPR9v+48BWwaO39z6JEkrZGi4V9WjwJEkL25dO4H7gYPA7ta3G7i1tQ8CV7e7ZnYAJwaWbyRJK2DdiOP+DPhQknOAh4BrmP2L4eYke4CHgSvb2NuBy4Bp4Mk2VpK0gkYK96q6F5icZ9fOecYWcO3SypIkLYWfUJWkDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nq0EjhnuSbSb6c5N4kU63vvCSHkjzYns9t/UlyfZLpJPcl2b6cE5Ak/byFXLn/TlVdVFWTbXsfcLiqtgGH2zbApcC29tgL3DCuYiVJo1nKsswu4EBrHwCuGOi/qWbdBaxPsnEJryNJWqBRw72A/0pyT5K9rW9DVT3S2o8CG1p7E3Bk4Nijre9pkuxNMpVkamZmZhGlS5JOZd2I4367qo4l+WXgUJKvDe6sqkpSC3nhqtoP7AeYnJxc0LGSpNMb6cq9qo615+PAJ4CLgcdOLre05+Nt+DFgy8Dhm1ufJGmFDA33JM9J8ryTbeB3ga8AB4Hdbdhu4NbWPghc3e6a2QGcGFi+kSStgFGWZTYAn0hycvy/VtWnk3weuDnJHuBh4Mo2/nbgMmAaeBK4ZuxVS5JOa2i4V9VDwEvn6X8c2DlPfwHXjqU6SdKi+AlVSeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR0y3CWpQ4a7JHXIcJekDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1aORwT3JWki8mua1tvzDJ3Ummk3w0yTmt/5lte7rt37pMtUuSTmEhV+5vAR4Y2H438N6qehHwHWBP698DfKf1v7eNkyStoJHCPclm4HLgA207wGuAW9qQA8AVrb2rbdP272zjJUkrZNQr978D/hL4v7Z9PvBEVT3Vto8Cm1p7E3AEoO0/0cY/TZK9SaaSTM3MzCyueknSvIaGe5LfA45X1T3jfOGq2l9Vk1U1OTExMc4fLUm/8NaNMOaVwOuSXAY8C/gl4H3A+iTr2tX5ZuBYG38M2AIcTbIOeD7w+NgrlySd0tAr96r6q6raXFVbgauAO6rqD4A7gde3YbuBW1v7YNum7b+jqmqsVUuSTmsp97m/HXhbkmlm19RvbP03Aue3/rcB+5ZWoiRpoUZZlvmpqvoM8JnWfgi4eJ4xPwTeMIbaJEmL5CdUJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR0y3CWpQ4a7JHXIcJekDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUoaHhnuRZST6X5EtJvprkXa3/hUnuTjKd5KNJzmn9z2zb023/1mWegyRpjlGu3H8EvKaqXgpcBFySZAfwbuC9VfUi4DvAnjZ+D/Cd1v/eNk6StIKGhnvN+n7bPLs9CngNcEvrPwBc0dq72jZt/84kGVfBkqThRlpzT3JWknuB48Ah4BvAE1X1VBtyFNjU2puAIwBt/wng/DHWLEkaYqRwr6qfVNVFwGbgYuAlS33hJHuTTCWZmpmZWeqPkyQNWNDdMlX1BHAn8ApgfZJ1bddm4FhrHwO2ALT9zwcen+dn7a+qyaqanJiYWFz1kqR5jXK3zESS9a39bOC1wAPMhvzr27DdwK2tfbBt0/bfUVU1xpolSUOsGz6EjcCBJGcx+5fBzVV1W5L7gY8k+Vvgi8CNbfyNwL8kmQa+DVy1DHVLkk5jaLhX1X3Ay+bpf4jZ9fe5/T8E3jCW6iRJi+InVCWpQ4a7JHXIcJekDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6NDTck2xJcmeS+5N8NclbWv95SQ4lebA9n9v6k+T6JNNJ7kuyfbknIUl6ulGu3J8C/qKqLgR2ANcmuRDYBxyuqm3A4bYNcCmwrT32AjeMvWpJ0mkNDfeqeqSqvtDa3wMeADYBu4ADbdgB4IrW3gXcVLPuAtYn2TjuwiVJp7agNfckW4GXAXcDG6rqkbbrUWBDa28CjgwcdrT1zf1Ze5NMJZmamZlZaN2SpNMYOdyTPBf4GPDWqvru4L6qKqAW8sJVtb+qJqtqcmJiYiGHSpKGGCnck5zNbLB/qKo+3rofO7nc0p6Pt/5jwJaBwze3PknSChnlbpkANwIPVNV7BnYdBHa39m7g1oH+q9tdMzuAEwPLN5KkFbBuhDGvBN4EfDnJva3vHcB1wM1J9gAPA1e2fbcDlwHTwJPANeMsWJI03NBwr6r/BnKK3TvnGV/AtUusS5K0BH5CVZI6ZLhLUocMd0nq0ChvqKojW/d9crVLkLQCvHKXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR3yK3+1YpbydcPfvO7yMVYi9c8rd0nq0NBwT/LBJMeTfGWg77wkh5I82J7Pbf1Jcn2S6ST3Jdm+nMVLkuY3ypX7PwOXzOnbBxyuqm3A4bYNcCmwrT32AjeMp0xJ0kIMDfeq+izw7Tndu4ADrX0AuGKg/6aadRewPsnGMdUqSRrRYtfcN1TVI639KLChtTcBRwbGHW19PyfJ3iRTSaZmZmYWWYYkaT5LfkO1qgqoRRy3v6omq2pyYmJiqWVIkgYsNtwfO7nc0p6Pt/5jwJaBcZtbnyRpBS023A8Cu1t7N3DrQP/V7a6ZHcCJgeUbSdIKGfohpiQfBl4NXJDkKPBO4Drg5iR7gIeBK9vw24HLgGngSeCaZahZkjTE0HCvqjeeYtfOecYWcO1Si5IkLY2fUJWkDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHRr6lb/SmWDrvk8u+thvXnf5GCuR1gav3CWpQ165q3te9esXkVfuktQhr9yl01jKVT945a/V45W7JHXIcJekDi3LskySS4D3AWcBH6iq65bjdX5RLXWpQCvHN3O1WsZ+5Z7kLOAfgUuBC4E3Jrlw3K8jSTq15bhyvxiYrqqHAJJ8BNgF3L8MryV1a63+huZvHGeG5Qj3TcCRge2jwG/NHZRkL7C3bX4/ydcX+XoXAN9a5LFnGudy5ullHrBCc8m7l/sVAM/LSS841Y5VuxWyqvYD+5f6c5JMVdXkGEpadc7lzNPLPMC5nKmWay7LcbfMMWDLwPbm1idJWiHLEe6fB7YleWGSc4CrgIPL8DqSpFMY+7JMVT2V5E+B/2T2VsgPVtVXx/06A5a8tHMGcS5nnl7mAc7lTLUsc0lVLcfPlSStIj+hKkkdMtwlqUNrJtyTXJLk60mmk+ybZ/8zk3y07b87ydZVKHMkI8zlzUlmktzbHn+4GnUOk+SDSY4n+cop9ifJ9W2e9yXZvtI1jmqEubw6yYmBc/LXK13jKJJsSXJnkvuTfDXJW+YZsybOy4hzWSvn5VlJPpfkS20u75pnzHgzrKrO+Aezb8x+A/g14BzgS8CFc8b8CfD+1r4K+Ohq172EubwZ+IfVrnWEubwK2A585RT7LwM+BQTYAdy92jUvYS6vBm5b7TpHmMdGYHtrPw/4n3n+fK2J8zLiXNbKeQnw3NY+G7gb2DFnzFgzbK1cuf/0Kw2q6sfAya80GLQLONDatwA7k2QFaxzVKHNZE6rqs8C3TzNkF3BTzboLWJ9k48pUtzAjzGVNqKpHquoLrf094AFmPzU+aE2clxHnsia0/9bfb5tnt8fcu1nGmmFrJdzn+0qDuSf5p2Oq6ingBHD+ilS3MKPMBeD326/MtyTZMs/+tWDUua4Vr2i/Vn8qya+vdjHDtF/rX8bsVeKgNXdeTjMXWCPnJclZSe4FjgOHquqU52UcGbZWwv0XzX8AW6vqN4FD/Oxvc62eLwAvqKqXAn8P/PvqlnN6SZ4LfAx4a1V9d7XrWYohc1kz56WqflJVFzH7qf2Lk/zGcr7eWgn3Ub7S4KdjkqwDng88viLVLczQuVTV41X1o7b5AeDlK1TbuHXzVRRV9d2Tv1ZX1e3A2UkuWOWy5pXkbGbD8ENV9fF5hqyZ8zJsLmvpvJxUVU8AdwKXzNk11gxbK+E+ylcaHAR2t/brgTuqvTNxhhk6lznrn69jdq1xLToIXN3uztgBnKiqR1a7qMVI8isn1z+TXMzs/ztn3MVDq/FG4IGqes8phq2J8zLKXNbQeZlIsr61nw28FvjanGFjzbA18Q9k1ym+0iDJ3wBTVXWQ2T8E/5Jkmtk3xq5avYpPbcS5/HmS1wFPMTuXN69awaeR5MPM3q1wQZKjwDuZfaOIqno/cDuzd2ZMA08C16xOpcONMJfXA3+c5Cngf4GrztCLh1cCbwK+3NZ3Ad4B/CqsufMyylzWynnZCBzI7D9m9Azg5qq6bTkzzK8fkKQOrZVlGUnSAhjuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUP/D2YsydU15EV0AAAAAElFTkSuQmCC\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGsCAYAAAAPJKchAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAfIklEQVR4nO3de2yV533A8d8BgkkX7ISm2AYcoKGhuRpKuZisBTYaj7Io1rSMRVphbcnWCqpQtlWwS6K0m5woIWHqWGiUBSuNGCnNAA3aBAoFlOCo4qZClrKSC5AEO4mS2MHtHGa/+yOKOzeY+DiYx5fPR3r/8OvnPec5j46Ov3rPe45zWZZlAQCQyIDUEwAA+jcxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJNWrYmT37t1x4403xogRIyKXy8XGjRvzvo0sy+Lee++NK664IgoKCmLkyJHxT//0T+d+sgBApwxKPYF8NDU1RXl5eXzlK1+JP/qjP+rSbdx2222xdevWuPfee+Paa6+NN998M958881zPFMAoLNyvfUf5eVyudiwYUNUVVW17Wtubo6/+7u/i3//93+Pt99+O6655pq4++67Y+bMmRER8dxzz8V1110Xhw8fjvHjx6eZOADQTq96m+bDLF68OGpra2PdunXx85//PG6++eb4gz/4g/jlL38ZERH/+Z//GZ/85Cdj8+bNMXbs2BgzZkwsXLjQmREASKjPxMjx48djzZo1sX79+vjc5z4Xl19+efz1X/91/O7v/m6sWbMmIiJeeOGFOHbsWKxfvz4eeeSRqKmpiX379sUf//EfJ549APRfveqakbM5dOhQtLS0xBVXXNFuf3Nzc3z84x+PiIjW1tZobm6ORx55pG3cv/3bv8WkSZPiyJEj3roBgAT6TIycOnUqBg4cGPv27YuBAwe2+91FF10UERGlpaUxaNCgdsFy5ZVXRsR7Z1bECACcf30mRiZOnBgtLS3x2muvxec+97kzjrn++uvjf//3f+P555+Pyy+/PCIi/vu//zsiIkaPHn3e5goA/Eav+jTNqVOn4ujRoxHxXnzcd999MWvWrBg2bFhcdtll8Wd/9mfx9NNPx4oVK2LixInx+uuvx/bt2+O6666LuXPnRmtra0yePDkuuuiiWLlyZbS2tsaiRYuisLAwtm7dmvjRAUD/1KtiZOfOnTFr1qwP7F+wYEHU1NTE6dOn4x//8R/jkUceiVdeeSUuvfTSmDZtWtx5551x7bXXRkTEq6++Gt/4xjdi69at8Tu/8zsxZ86cWLFiRQwbNux8PxwAIHpZjAAAfU+f+WgvANA7iREAIKle8Wma1tbWePXVV2Po0KGRy+VSTwcA6IQsy+Kdd96JESNGxIABHZ//6BUx8uqrr0ZZWVnqaQAAXXDixIkYNWpUh7/vFTEydOjQiHjvwRQWFiaeDQDQGY2NjVFWVtb2d7wjvSJG3n9rprCwUIwAQC/zYZdYuIAVAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJDUoNQTgJ5szLIt3XK7L901t1tuF6A3cmYEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApPKKkerq6pg8eXIMHTo0hg8fHlVVVXHkyJGzHlNTUxO5XK7dNmTIkI80aQCg78grRnbt2hWLFi2KZ555JrZt2xanT5+OG264IZqams56XGFhYZw8ebJtO3bs2EeaNADQdwzKZ/ATTzzR7ueampoYPnx47Nu3Lz7/+c93eFwul4uSkpKuzRAA6NM+0jUjDQ0NERExbNiws447depUjB49OsrKyuKmm26KZ5999qzjm5ubo7Gxsd0GAPRNXY6R1tbWWLJkSVx//fVxzTXXdDhu/Pjx8fDDD8emTZvi0UcfjdbW1pg+fXq8/PLLHR5TXV0dRUVFbVtZWVlXpwkA9HC5LMuyrhz49a9/PX784x/HU089FaNGjer0cadPn44rr7wybrnllvjOd75zxjHNzc3R3Nzc9nNjY2OUlZVFQ0NDFBYWdmW60CVjlm3pltt96a653XK7AD1JY2NjFBUVfejf77yuGXnf4sWLY/PmzbF79+68QiQi4oILLoiJEyfG0aNHOxxTUFAQBQUFXZkaANDL5PU2TZZlsXjx4tiwYUPs2LEjxo4dm/cdtrS0xKFDh6K0tDTvYwGAvievMyOLFi2KtWvXxqZNm2Lo0KFRV1cXERFFRUVx4YUXRkTE/PnzY+TIkVFdXR0REd/+9rdj2rRpMW7cuHj77bfjnnvuiWPHjsXChQvP8UMBAHqjvGLkgQceiIiImTNnttu/Zs2a+PM///OIiDh+/HgMGPCbEy5vvfVW3HrrrVFXVxeXXHJJTJo0Kfbs2RNXXXXVR5s5ANAndPkC1vOpsxfAwLnmAlaAruvs32//mwYASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUnnFSHV1dUyePDmGDh0aw4cPj6qqqjhy5MiHHrd+/fr49Kc/HUOGDIlrr702fvSjH3V5wgBA35JXjOzatSsWLVoUzzzzTGzbti1Onz4dN9xwQzQ1NXV4zJ49e+KWW26Jr371q3HgwIGoqqqKqqqqOHz48EeePADQ++WyLMu6evDrr78ew4cPj127dsXnP//5M46ZN29eNDU1xebNm9v2TZs2LSZMmBCrV6/u1P00NjZGUVFRNDQ0RGFhYVenC3kbs2xLt9zuS3fN7ZbbBehJOvv3+yNdM9LQ0BAREcOGDetwTG1tbcyePbvdvsrKyqitre3wmObm5mhsbGy3AQB9U5djpLW1NZYsWRLXX399XHPNNR2Oq6uri+Li4nb7iouLo66ursNjqquro6ioqG0rKyvr6jQBgB6uyzGyaNGiOHz4cKxbt+5cziciIpYvXx4NDQ1t24kTJ875fQAAPcOgrhy0ePHi2Lx5c+zevTtGjRp11rElJSVRX1/fbl99fX2UlJR0eExBQUEUFBR0ZWoAQC+T15mRLMti8eLFsWHDhtixY0eMHTv2Q4+pqKiI7du3t9u3bdu2qKioyG+mAECflNeZkUWLFsXatWtj06ZNMXTo0LbrPoqKiuLCCy+MiIj58+fHyJEjo7q6OiIibrvttpgxY0asWLEi5s6dG+vWrYu9e/fGgw8+eI4fCgDQG+V1ZuSBBx6IhoaGmDlzZpSWlrZtjz32WNuY48ePx8mTJ9t+nj59eqxduzYefPDBKC8vjx/+8IexcePGs170CgD0H3mdGenMV5Ls3LnzA/tuvvnmuPnmm/O5KwCgn/C/aQCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAklXeM7N69O2688cYYMWJE5HK52Lhx41nH79y5M3K53Ae2urq6rs4ZAOhD8o6RpqamKC8vj1WrVuV13JEjR+LkyZNt2/Dhw/O9awCgDxqU7wFz5syJOXPm5H1Hw4cPj4svvjjv4wCAvu28XTMyYcKEKC0tjS984Qvx9NNPn3Vsc3NzNDY2ttsAgL6p22OktLQ0Vq9eHY8//ng8/vjjUVZWFjNnzoz9+/d3eEx1dXUUFRW1bWVlZd09TQAgkVyWZVmXD87lYsOGDVFVVZXXcTNmzIjLLrssvv/975/x983NzdHc3Nz2c2NjY5SVlUVDQ0MUFhZ2dbqQtzHLtnTL7b5019xuuV2AnqSxsTGKioo+9O933teMnAtTpkyJp556qsPfFxQUREFBwXmcEQCQSpLvGTl48GCUlpamuGsAoIfJ+8zIqVOn4ujRo20/v/jii3Hw4MEYNmxYXHbZZbF8+fJ45ZVX4pFHHomIiJUrV8bYsWPj6quvjv/5n/+Jhx56KHbs2BFbt249d48CAOi18o6RvXv3xqxZs9p+Xrp0aURELFiwIGpqauLkyZNx/Pjxtt+/++678Vd/9VfxyiuvxMc+9rG47rrr4ic/+Um72wAA+q+PdAHr+dLZC2DgXHMBK0DXdfbvt/9NAwAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASeUdI7t3744bb7wxRowYEblcLjZu3Pihx+zcuTM+85nPREFBQYwbNy5qamq6MFUAoC/KO0aampqivLw8Vq1a1anxL774YsydOzdmzZoVBw8ejCVLlsTChQvjySefzHuyAEDfMyjfA+bMmRNz5szp9PjVq1fH2LFjY8WKFRERceWVV8ZTTz0V999/f1RWVuZ79wBAH9Pt14zU1tbG7Nmz2+2rrKyM2traDo9pbm6OxsbGdhsA0Dd1e4zU1dVFcXFxu33FxcXR2NgYv/71r894THV1dRQVFbVtZWVl3T1NACCRHvlpmuXLl0dDQ0PbduLEidRTAgC6Sd7XjOSrpKQk6uvr2+2rr6+PwsLCuPDCC894TEFBQRQUFHT31ACAHqDbz4xUVFTE9u3b2+3btm1bVFRUdPddAwC9QN4xcurUqTh48GAcPHgwIt776O7Bgwfj+PHjEfHeWyzz589vG/+1r30tXnjhhfjWt74Vv/jFL+Jf//Vf4wc/+EF885vfPDePAADo1fKOkb1798bEiRNj4sSJERGxdOnSmDhxYtx+++0REXHy5Mm2MImIGDt2bGzZsiW2bdsW5eXlsWLFinjooYd8rBcAiIiIXJZlWepJfJjGxsYoKiqKhoaGKCwsTD0d+pExy7Z0y+2+dNfcbrldgJ6ks3+/e+SnaQCA/kOMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUoNSTwA+qjHLtqSeAgAfgTMjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJBUl2Jk1apVMWbMmBgyZEhMnTo1fvazn3U4tqamJnK5XLttyJAhXZ4wANC35B0jjz32WCxdujTuuOOO2L9/f5SXl0dlZWW89tprHR5TWFgYJ0+ebNuOHTv2kSYNAPQdg/I94L777otbb701vvzlL0dExOrVq2PLli3x8MMPx7Jly854TC6Xi5KSko82U+hDxizb0m23/dJdc7vttgG6Q15nRt59993Yt29fzJ49+zc3MGBAzJ49O2prazs87tSpUzF69OgoKyuLm266KZ599tmz3k9zc3M0Nja22wCAvimvGHnjjTeipaUliouL2+0vLi6Ourq6Mx4zfvz4ePjhh2PTpk3x6KOPRmtra0yfPj1efvnlDu+nuro6ioqK2raysrJ8pgkA9CLd/mmaioqKmD9/fkyYMCFmzJgR//Ef/xGf+MQn4nvf+16HxyxfvjwaGhrathMnTnT3NAGARPK6ZuTSSy+NgQMHRn19fbv99fX1nb4m5IILLoiJEyfG0aNHOxxTUFAQBQUF+UwNAOil8jozMnjw4Jg0aVJs3769bV9ra2ts3749KioqOnUbLS0tcejQoSgtLc1vpgBAn5T3p2mWLl0aCxYsiM9+9rMxZcqUWLlyZTQ1NbV9umb+/PkxcuTIqK6ujoiIb3/72zFt2rQYN25cvP3223HPPffEsWPHYuHChef2kQAAvVLeMTJv3rx4/fXX4/bbb4+6urqYMGFCPPHEE20XtR4/fjwGDPjNCZe33norbr311qirq4tLLrkkJk2aFHv27Imrrrrq3D0KAKDXymVZlqWexIdpbGyMoqKiaGhoiMLCwtTToYfpzu/s6I18zwjQU3T277f/TQMAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSg1JPADi3xizb0i23+9Jdc7vldgGcGQEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUr4OHuiU7vqa+QhfNQ/9nTMjAEBSYgQASEqMAABJuWYESK67rkdxLQr0Ds6MAABJOTMC9Fk+AQS9gzMjAEBSXTozsmrVqrjnnnuirq4uysvL47vf/W5MmTKlw/Hr16+Pf/iHf4iXXnopPvWpT8Xdd98dX/ziF7s8aYDUXOcC507eMfLYY4/F0qVLY/Xq1TF16tRYuXJlVFZWxpEjR2L48OEfGL9nz5645ZZborq6Ov7wD/8w1q5dG1VVVbF///645pprzsmDoOfrztPl0Jd4a4n+KJdlWZbPAVOnTo3JkyfHv/zLv0RERGtra5SVlcU3vvGNWLZs2QfGz5s3L5qammLz5s1t+6ZNmxYTJkyI1atXd+o+Gxsbo6ioKBoaGqKwsDCf6dJDiBFIT4xwvnX273deZ0befffd2LdvXyxfvrxt34ABA2L27NlRW1t7xmNqa2tj6dKl7fZVVlbGxo0bO7yf5ubmaG5ubvu5oaEhIt57UPROrc2/Sj0F6Pcu++b61FPI2+E7K1NPgY/g/b/bH3beI68YeeONN6KlpSWKi4vb7S8uLo5f/OIXZzymrq7ujOPr6uo6vJ/q6uq48847P7C/rKwsn+kC0MsVrUw9A86Fd955J4qKijr8fY/8aO/y5cvbnU1pbW2NN998Mz7+8Y9HLpc7Z/fT2NgYZWVlceLECW//fAhrlR/r1XnWqvOsVedZq87rzrXKsizeeeedGDFixFnH5RUjl156aQwcODDq6+vb7a+vr4+SkpIzHlNSUpLX+IiIgoKCKCgoaLfv4osvzmeqeSksLPRk7SRrlR/r1XnWqvOsVedZq87rrrU62xmR9+X1PSODBw+OSZMmxfbt29v2tba2xvbt26OiouKMx1RUVLQbHxGxbdu2DscDAP1L3m/TLF26NBYsWBCf/exnY8qUKbFy5cpoamqKL3/5yxERMX/+/Bg5cmRUV1dHRMRtt90WM2bMiBUrVsTcuXNj3bp1sXfv3njwwQfP7SMBAHqlvGNk3rx58frrr8ftt98edXV1MWHChHjiiSfaLlI9fvx4DBjwmxMu06dPj7Vr18bf//3fx9/+7d/Gpz71qdi4cWOP+I6RgoKCuOOOOz7wlhAfZK3yY706z1p1nrXqPGvVeT1hrfL+nhEAgHPJ/6YBAJISIwBAUmIEAEhKjAAASfX5GFm1alWMGTMmhgwZElOnTo2f/exnZx2/fv36+PSnPx1DhgyJa6+9Nn70ox+dp5mml89a1dTURC6Xa7cNGTLkPM42nd27d8eNN94YI0aMiFwud9b/s/S+nTt3xmc+85koKCiIcePGRU1NTbfPsyfId6127tz5gedVLpc767+P6Cuqq6tj8uTJMXTo0Bg+fHhUVVXFkSNHPvS4/via1ZW16q+vWQ888EBcd911bV9oVlFRET/+8Y/PekyK51SfjpHHHnssli5dGnfccUfs378/ysvLo7KyMl577bUzjt+zZ0/ccsst8dWvfjUOHDgQVVVVUVVVFYcPHz7PMz//8l2riPe+re/kyZNt27Fjx87jjNNpamqK8vLyWLVqVafGv/jiizF37tyYNWtWHDx4MJYsWRILFy6MJ598sptnml6+a/W+I0eOtHtuDR8+vJtm2HPs2rUrFi1aFM8880xs27YtTp8+HTfccEM0NTV1eEx/fc3qylpF9M/XrFGjRsVdd90V+/bti71798bv/d7vxU033RTPPvvsGccne05lfdiUKVOyRYsWtf3c0tKSjRgxIquurj7j+D/5kz/J5s6d227f1KlTs7/8y7/s1nn2BPmu1Zo1a7KioqLzNLueKyKyDRs2nHXMt771rezqq69ut2/evHlZZWVlN86s5+nMWv30pz/NIiJ76623zsucerLXXnsti4hs165dHY7pz69Z/19n1spr1m9ccskl2UMPPXTG36V6TvXZMyPvvvtu7Nu3L2bPnt22b8CAATF79uyora094zG1tbXtxkdEVFZWdji+r+jKWkVEnDp1KkaPHh1lZWVnLe3+rr8+rz6KCRMmRGlpaXzhC1+Ip59+OvV0kmhoaIiIiGHDhnU4xnPrPZ1ZqwivWS0tLbFu3bpoamrq8F+ypHpO9dkYeeONN6KlpaXtm2HfV1xc3OH7z3V1dXmN7yu6slbjx4+Phx9+ODZt2hSPPvpotLa2xvTp0+Pll18+H1PuVTp6XjU2Nsavf/3rRLPqmUpLS2P16tXx+OOPx+OPPx5lZWUxc+bM2L9/f+qpnVetra2xZMmSuP7668/6bdX99TXr/+vsWvXn16xDhw7FRRddFAUFBfG1r30tNmzYEFddddUZx6Z6TuX9dfAQ8d4/QPz/ZT19+vS48sor43vf+1585zvfSTgzerPx48fH+PHj236ePn16PP/883H//ffH97///YQzO78WLVoUhw8fjqeeeir1VHq8zq5Vf37NGj9+fBw8eDAaGhrihz/8YSxYsCB27drVYZCk0GfPjFx66aUxcODAqK+vb7e/vr4+SkpKznhMSUlJXuP7iq6s1W+74IILYuLEiXH06NHumGKv1tHzqrCwMC688MJEs+o9pkyZ0q+eV4sXL47NmzfHT3/60xg1atRZx/bX16z35bNWv60/vWYNHjw4xo0bF5MmTYrq6uooLy+Pf/7nfz7j2FTPqT4bI4MHD45JkybF9u3b2/a1trbG9u3bO3yvrKKiot34iIht27Z1OL6v6Mpa/baWlpY4dOhQlJaWdtc0e63++rw6Vw4ePNgvnldZlsXixYtjw4YNsWPHjhg7duyHHtNfn1tdWavf1p9fs1pbW6O5ufmMv0v2nOrWy2MTW7duXVZQUJDV1NRk//Vf/5X9xV/8RXbxxRdndXV1WZZl2Ze+9KVs2bJlbeOffvrpbNCgQdm9996bPffcc9kdd9yRXXDBBdmhQ4dSPYTzJt+1uvPOO7Mnn3wye/7557N9+/Zlf/qnf5oNGTIke/bZZ1M9hPPmnXfeyQ4cOJAdOHAgi4jsvvvuyw4cOJAdO3Ysy7IsW7ZsWfalL32pbfwLL7yQfexjH8v+5m/+JnvuueeyVatWZQMHDsyeeOKJVA/hvMl3re6///5s48aN2S9/+cvs0KFD2W233ZYNGDAg+8lPfpLqIZw3X//617OioqJs586d2cmTJ9u2X/3qV21jvGa9pytr1V9fs5YtW5bt2rUre/HFF7Of//zn2bJly7JcLpdt3bo1y7Ke85zq0zGSZVn23e9+N7vsssuywYMHZ1OmTMmeeeaZtt/NmDEjW7BgQbvxP/jBD7IrrrgiGzx4cHb11VdnW7ZsOc8zTieftVqyZEnb2OLi4uyLX/xitn///gSzPv/e//jpb2/vr8+CBQuyGTNmfOCYCRMmZIMHD84++clPZmvWrDnv804h37W6++67s8svvzwbMmRINmzYsGzmzJnZjh070kz+PDvTOkVEu+eK16z3dGWt+utr1le+8pVs9OjR2eDBg7NPfOIT2e///u+3hUiW9ZznVC7Lsqx7z70AAHSsz14zAgD0DmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgqf8DgyUK1RABovAAAAAASUVORK5CYII=\n", "text/plain": [ - "<Figure size 432x288 with 1 Axes>" + "<Figure size 640x480 with 1 Axes>" ] }, - "metadata": { - "needs_background": "light" - }, + "metadata": {}, "output_type": "display_data" } ], @@ -105,7 +103,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 6, "id": "6c2a4a1a", "metadata": {}, "outputs": [], @@ -115,17 +113,17 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 7, "id": "a7231e19", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "<Array [1.07, 1.04, 1.04, ... 1.21, 1.51, 1.19] type='1098 * float32'>" + "<Array [1.07, 1.04, 1.04, ... 0.995, 0.601] type='4343693 * float32'>" ] }, - "execution_count": 8, + "execution_count": 7, "metadata": {}, "output_type": "execute_result" } @@ -137,20 +135,18 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 8, "id": "abfe684d", "metadata": {}, "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAOzUlEQVR4nO3db4xcV33G8e9DnAAqFEO8dSPbYVMRqUqrAqkVjKiqlIgqfyocqQEFVcREriy1qQqiVTG8KKLqC/OGlLQVyCKoDqKQNEDjhtDWCkGoUhPYQAj5A2WJEsVWwCYkBpRCZfrriz1Oh2XXM7s7u+s5+X6k0Zx7zpmZ38l1Ht+9e+c6VYUkqS/PW+8CJEnjZ7hLUocMd0nqkOEuSR0y3CWpQxvWuwCATZs21fT09HqXIUkT5d577/1eVU0tNHZahPv09DQzMzPrXYYkTZQkjy025mkZSeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUodGCvckjyb5epL7ksy0vpclOZTkW+35pa0/SW5IMpvk/iQXruYCJEk/bynfUP2dqvrewPZe4M6q2pdkb9t+F3AZcH57vAb4UHtWB6b3fnbZr3103xVjrETSqazktMxO4EBrHwCuHOi/qebcDWxMcs4KPkeStESjhnsB/57k3iR7Wt/mqnqitb8DbG7tLcDjA6893PokSWtk1NMyv1VVR5L8EnAoyTcGB6uqkizpH2Ntf0nsATj33HOX8lJJ0hAjHblX1ZH2fBT4DHAR8N2Tp1va89E2/QiwbeDlW1vf/PfcX1Xbq2r71NSCd6yUJC3T0HBP8gtJXnyyDfwu8ABwENjVpu0Cbmvtg8A17aqZHcDxgdM3kqQ1MMppmc3AZ5KcnP+PVfWvSb4M3JJkN/AY8OY2/w7gcmAWeAa4duxVS5JOaWi4V9UjwCsX6H8SuGSB/gKuG0t1kqRl8RuqktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR0y3CWpQ4a7JHXIcJekDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA6NHO5Jzkjy1SS3t+3zktyTZDbJzUnOav3Pb9uzbXx6lWqXJC1iKUfubwceHth+P3B9Vb0CeArY3fp3A0+1/uvbPEnSGhop3JNsBa4APtK2A7weuLVNOQBc2do72zZt/JI2X5K0RkY9cv8b4C+A/23bZwNPV9WJtn0Y2NLaW4DHAdr48Tb/ZyTZk2QmycyxY8eWV70kaUFDwz3J7wFHq+recX5wVe2vqu1VtX1qamqcby1Jz3kbRpjzOuCNSS4HXgD8IvBBYGOSDe3ofCtwpM0/AmwDDifZALwEeHLslUuSFjU03Kvq3cC7AZJcDPx5Vf1Bkn8CrgI+CewCbmsvOdi2/7ONf76qauyVa+JM7/3ssl/76L4rxliJ1L+VXOf+LuCdSWaZO6d+Y+u/ETi79b8T2LuyEiVJSzXKaZlnVdUXgC+09iPARQvM+THwpjHUJklaJr+hKkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR0y3CWpQ4a7JHXIcJekDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SerQ0HBP8oIkX0rytSQPJnlf6z8vyT1JZpPcnOSs1v/8tj3bxqdXeQ2SpHlGOXL/CfD6qnol8Crg0iQ7gPcD11fVK4CngN1t/m7gqdZ/fZsnSVpDQ8O95vyobZ7ZHgW8Hri19R8ArmztnW2bNn5JkoyrYEnScCOdc09yRpL7gKPAIeDbwNNVdaJNOQxsae0twOMAbfw4cPYC77knyUySmWPHjq1oEZKknzVSuFfVT6vqVcBW4CLgV1f6wVW1v6q2V9X2qamplb6dJGnAkq6WqaqngbuA1wIbk2xoQ1uBI619BNgG0MZfAjw5jmIlSaMZ5WqZqSQbW/uFwBuAh5kL+avatF3Aba19sG3Txj9fVTXGmiVJQ2wYPoVzgANJzmDuL4Nbqur2JA8Bn0zy18BXgRvb/BuBjyWZBb4PXL0KdUuSTmFouFfV/cCrF+h/hLnz7/P7fwy8aSzVSZKWxW+oSlKHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR0y3CWpQ4a7JHXIcJekDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjo0NNyTbEtyV5KHkjyY5O2t/2VJDiX5Vnt+aetPkhuSzCa5P8mFq70ISdLPGuXI/QTwZ1V1AbADuC7JBcBe4M6qOh+4s20DXAac3x57gA+NvWpJ0ikNDfeqeqKqvtLaPwQeBrYAO4EDbdoB4MrW3gncVHPuBjYmOWfchUuSFrekc+5JpoFXA/cAm6vqiTb0HWBza28BHh942eHWN/+99iSZSTJz7NixpdYtSTqFkcM9yYuATwHvqKofDI5VVQG1lA+uqv1Vtb2qtk9NTS3lpZKkIUYK9yRnMhfsH6+qT7fu75483dKej7b+I8C2gZdvbX2SpDUyytUyAW4EHq6qDwwMHQR2tfYu4LaB/mvaVTM7gOMDp28kSWtgwwhzXge8Ffh6kvta33uAfcAtSXYDjwFvbmN3AJcDs8AzwLXjLFiSNNzQcK+q/wCyyPAlC8wv4LoV1iVJWgG/oSpJHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR0y3CWpQ4a7JHXIcJekDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nq0NBwT/LRJEeTPDDQ97Ikh5J8qz2/tPUnyQ1JZpPcn+TC1SxekrSwUY7c/wG4dF7fXuDOqjofuLNtA1wGnN8ee4APjadMSdJSDA33qvoi8P153TuBA619ALhyoP+mmnM3sDHJOWOqVZI0ouWec99cVU+09neAza29BXh8YN7h1vdzkuxJMpNk5tixY8ssQ5K0kBX/QrWqCqhlvG5/VW2vqu1TU1MrLUOSNGC54f7dk6db2vPR1n8E2DYwb2vrkyStoeWG+0FgV2vvAm4b6L+mXTWzAzg+cPpGkrRGNgybkOQTwMXApiSHgfcC+4BbkuwGHgPe3KbfAVwOzALPANeuQs2SpCGGhntVvWWRoUsWmFvAdSstSpK0Mn5DVZI6ZLhLUoeGnpZRX6b3fna9S5C0Bjxyl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR0y3CWpQ4a7JHXIcJekDhnuktQhw12SOmS4S1KHDHdJ6pB3hZxAz8U7O65kzY/uu2KMlUiTwSN3SeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR3y9gPr5Ll4CwFJa8cjd0nqkOEuSR0y3CWpQ4a7JHVoVX6hmuRS4IPAGcBHqmrfanyONArvBa/norEfuSc5A/h74DLgAuAtSS4Y9+dIkha3GkfuFwGzVfUIQJJPAjuBh1bhs9aVlzP2b6X72CN/rZfVCPctwOMD24eB18yflGQPsKdt/ijJN5f5eZuA7y3ztacb13L6WdE68v4xVrJyvewTcC0nvXyxgXX7ElNV7Qf2r/R9ksxU1fYxlLTuXMvpp5d1gGs5Xa3WWlbjapkjwLaB7a2tT5K0RlYj3L8MnJ/kvCRnAVcDB1fhcyRJixj7aZmqOpHkT4B/Y+5SyI9W1YPj/pwBKz61cxpxLaefXtYBruV0tSprSVWtxvtKktaR31CVpA4Z7pLUoYkJ9ySXJvlmktkkexcYf36Sm9v4PUmm16HMkYywlrclOZbkvvb4w/Woc5gkH01yNMkDi4wnyQ1tnfcnuXCtaxzVCGu5OMnxgX3yl2td4yiSbEtyV5KHkjyY5O0LzJmI/TLiWiZlv7wgyZeSfK2t5X0LzBlvhlXVaf9g7hez3wZ+BTgL+Bpwwbw5fwx8uLWvBm5e77pXsJa3AX+33rWOsJbfBi4EHlhk/HLgc0CAHcA9613zCtZyMXD7etc5wjrOAS5s7RcD/7XAn6+J2C8jrmVS9kuAF7X2mcA9wI55c8aaYZNy5P7sLQ2q6n+Ak7c0GLQTONDatwKXJMka1jiqUdYyEarqi8D3TzFlJ3BTzbkb2JjknLWpbmlGWMtEqKonquorrf1D4GHmvjU+aCL2y4hrmQjtv/WP2uaZ7TH/apaxZtikhPtCtzSYv5OfnVNVJ4DjwNlrUt3SjLIWgN9vPzLfmmTbAuOTYNS1TorXth+rP5fk19a7mGHaj/WvZu4ocdDE7ZdTrAUmZL8kOSPJfcBR4FBVLbpfxpFhkxLuzzX/AkxX1W8Ah/j/v821fr4CvLyqXgn8LfDP61vOqSV5EfAp4B1V9YP1rmclhqxlYvZLVf20ql7F3Lf2L0ry66v5eZMS7qPc0uDZOUk2AC8BnlyT6pZm6Fqq6smq+knb/Ajwm2tU27h1cyuKqvrByR+rq+oO4Mwkm9a5rAUlOZO5MPx4VX16gSkTs1+GrWWS9stJVfU0cBdw6byhsWbYpIT7KLc0OAjsau2rgM9X+83EaWboWuad/3wjc+caJ9FB4Jp2dcYO4HhVPbHeRS1Hkl8+ef4zyUXM/b9z2h08tBpvBB6uqg8sMm0i9ssoa5mg/TKVZGNrvxB4A/CNedPGmmHrdlfIpahFbmmQ5K+Amao6yNwfgo8lmWXuF2NXr1/FixtxLX+a5I3ACebW8rZ1K/gUknyCuasVNiU5DLyXuV8UUVUfBu5g7sqMWeAZ4Nr1qXS4EdZyFfBHSU4A/w1cfZoePLwOeCvw9XZ+F+A9wLkwcftllLVMyn45BziQuX/M6HnALVV1+2pmmLcfkKQOTcppGUnSEhjuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUP/B+tC2+0i88UDAAAAAElFTkSuQmCC\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGsCAYAAAD+L/ysAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAArDUlEQVR4nO3de3BUZZ7/8U8nmA6oaUDITSMJF8OIkDAoMag/YG0N2SxldmoUqBmJGdEdFi2YeJnEVdDV3aAiojNZMioQmFkFWQVrhIkw0UAhQYpASnCQBQw3SYeL0k16xuCkz+8Py3ZbEsgJJv3Qeb+qTmk/53uefM+ppvtTp0+fdliWZQkAAMBgUeFuAAAA4HwILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeBEXWDZu3KiJEycqOTlZDodDq1evtj2HZVmaN2+errnmGjmdTl155ZX6j//4jx++WQAA0C49wt3AD83v9ysjI0O/+MUv9JOf/KRDc8ycOVPr1q3TvHnzNHz4cH3xxRf64osvfuBOAQBAezki+ccPHQ6HVq1apfz8/OBYc3Oz/u3f/k1vvPGGTp06peuuu07PPvusxo0bJ0navXu3RowYoV27dik9PT08jQMAgBAR95HQ+TzwwAOqqanR8uXL9fHHH+vOO+/UhAkTtHfvXknSH//4Rw0cOFDvvvuu0tLSlJqaqmnTpnGGBQCAMOpWgeXQoUNasmSJVq5cqVtuuUWDBg3Sww8/rJtvvllLliyRJH322Wc6ePCgVq5cqWXLlqmiokK1tbX66U9/GubuAQDoviLuGpZz2blzp1paWnTNNdeEjDc3N+uKK66QJAUCATU3N2vZsmXBukWLFmnUqFHas2cPHxMBABAG3SqwNDU1KTo6WrW1tYqOjg5Zd9lll0mSkpKS1KNHj5BQ86Mf/UjSN2doCCwAAHS9bhVYRo4cqZaWFh07dky33HJLqzU33XST/v73v2v//v0aNGiQJOl///d/JUkDBgzosl4BAMB3Iu5bQk1NTdq3b5+kbwLK/PnzNX78ePXt21dXX321fv7zn+vDDz/UCy+8oJEjR+r48eOqqqrSiBEjlJeXp0AgoBtuuEGXXXaZFixYoEAgoBkzZiguLk7r1q0L894BANA9RVxgqa6u1vjx488aLygoUEVFhb7++ms988wzWrZsmT7//HP169dPN954o5566ikNHz5cknT06FE9+OCDWrdunS699FLl5ubqhRdeUN++fbt6dwAAgCIwsAAAgMjTrb7WDAAALk4EFgAAYLyI+JZQIBDQ0aNHdfnll8vhcIS7HQAA0A6WZen06dNKTk5WVNS5z6FERGA5evSoUlJSwt0GAADogMOHD+uqq646Z01EBJbLL79c0jc7HBcXF+ZuAABAe/h8PqWkpATfx88lIgLLtx8DxcXFEVgAALjItOdyDi66BQAAxiOwAAAA4xFYAACA8WwFltLSUt1www26/PLLFR8fr/z8fO3Zs+e8261cuVJDhw5VbGyshg8frrVr14astyxLs2fPVlJSknr27Cm32629e/fa2xMAABCxbAWWDRs2aMaMGdqyZYvWr1+vr7/+Wrfffrv8fn+b22zevFlTpkzRvffeqx07dig/P1/5+fnatWtXsOa5557Tyy+/rPLycn300Ue69NJLlZOTo6+++qrjewYAACLGBf2W0PHjxxUfH68NGzbo//2//9dqzaRJk+T3+/Xuu+8Gx2688UZlZmaqvLxclmUpOTlZDz30kB5++GFJktfrVUJCgioqKjR58uTz9uHz+eRyueT1evmWEAAAFwk7798XdA2L1+uVpHP+inFNTY3cbnfIWE5OjmpqaiRJ9fX18ng8ITUul0tZWVnBmu9rbm6Wz+cLWQAAQOTqcGAJBAKaNWuWbrrpJl133XVt1nk8HiUkJISMJSQkyOPxBNd/O9ZWzfeVlpbK5XIFF+5yCwBAZOtwYJkxY4Z27dql5cuX/5D9tEtJSYm8Xm9wOXz4cJf3AAAAuk6H7nT7wAMP6N1339XGjRvPe+//xMRENTY2how1NjYqMTExuP7bsaSkpJCazMzMVud0Op1yOp0daR0AAFyEbJ1hsSxLDzzwgFatWqX3339faWlp590mOztbVVVVIWPr169Xdna2JCktLU2JiYkhNT6fTx999FGwBgAAdG+2zrDMmDFDr7/+ut555x1dfvnlwWtMXC6XevbsKUmaOnWqrrzySpWWlkqSZs6cqbFjx+qFF15QXl6eli9frm3btumVV16R9M3vB8yaNUvPPPOMhgwZorS0ND3xxBNKTk5Wfn7+D7irAADgYmUrsCxcuFCSNG7cuJDxJUuW6J577pEkHTp0SFFR3524GTNmjF5//XU9/vjjeuyxxzRkyBCtXr065ELdRx99VH6/X/fff79OnTqlm2++WZWVlYqNje3gbgEAgEhyQfdhMQX3YQEA4OJj5/27QxfdAvhGavGaTpv7wNy8TpsbAC42/PghAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADCe7cCyceNGTZw4UcnJyXI4HFq9evU56++55x45HI6zlmHDhgVrnnzyybPWDx061PbOAACAyGQ7sPj9fmVkZKisrKxd9S+99JIaGhqCy+HDh9W3b1/deeedIXXDhg0Lqdu0aZPd1gAAQITqYXeD3Nxc5ebmtrve5XLJ5XIFH69evVpffvmlCgsLQxvp0UOJiYl22wEAAN1Al1/DsmjRIrndbg0YMCBkfO/evUpOTtbAgQP1s5/9TIcOHWpzjubmZvl8vpAFAABEri4NLEePHtWf/vQnTZs2LWQ8KytLFRUVqqys1MKFC1VfX69bbrlFp0+fbnWe0tLS4Jkbl8ullJSUrmgfAACESZcGlqVLl6p3797Kz88PGc/NzdWdd96pESNGKCcnR2vXrtWpU6f05ptvtjpPSUmJvF5vcDl8+HAXdA8AAMLF9jUsHWVZlhYvXqy7775bMTEx56zt3bu3rrnmGu3bt6/V9U6nU06nszPaBAAABuqyMywbNmzQvn37dO+99563tqmpSfv371dSUlIXdAYAAExnO7A0NTWprq5OdXV1kqT6+nrV1dUFL5ItKSnR1KlTz9pu0aJFysrK0nXXXXfWuocfflgbNmzQgQMHtHnzZv3zP/+zoqOjNWXKFLvtAQCACGT7I6Ft27Zp/PjxwcdFRUWSpIKCAlVUVKihoeGsb/h4vV699dZbeumll1qd88iRI5oyZYpOnjyp/v376+abb9aWLVvUv39/u+0BAIAI5LAsywp3ExfK5/PJ5XLJ6/UqLi4u3O2gG0ktXtNpcx+Ym9dpcwOACey8f/NbQgAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgvB52N9i4caOef/551dbWqqGhQatWrVJ+fn6b9dXV1Ro/fvxZ4w0NDUpMTAw+Lisr0/PPPy+Px6OMjAz95je/0ejRo+22B0SM1OI1nTLvgbl5nTIvAHQm22dY/H6/MjIyVFZWZmu7PXv2qKGhIbjEx8cH161YsUJFRUWaM2eOtm/froyMDOXk5OjYsWN22wMAABHI9hmW3Nxc5ebm2v5D8fHx6t27d6vr5s+fr/vuu0+FhYWSpPLycq1Zs0aLFy9WcXGx7b8FAAAiS5ddw5KZmamkpCTddttt+vDDD4PjZ86cUW1trdxu93dNRUXJ7Xarpqam1bmam5vl8/lCFgAAELk6PbAkJSWpvLxcb731lt566y2lpKRo3Lhx2r59uyTpxIkTamlpUUJCQsh2CQkJ8ng8rc5ZWloql8sVXFJSUjp7NwAAQBjZ/kjIrvT0dKWnpwcfjxkzRvv379eLL76o3//+9x2as6SkREVFRcHHPp+P0AIAQATr9MDSmtGjR2vTpk2SpH79+ik6OlqNjY0hNY2NjSHfIvq/nE6nnE5np/cJAADMEJb7sNTV1SkpKUmSFBMTo1GjRqmqqiq4PhAIqKqqStnZ2eFoDwAAGMb2GZampibt27cv+Li+vl51dXXq27evrr76apWUlOjzzz/XsmXLJEkLFixQWlqahg0bpq+++kqvvfaa3n//fa1bty44R1FRkQoKCnT99ddr9OjRWrBggfx+f/BbQwAAoHuzHVi2bdsWciO4b68lKSgoUEVFhRoaGnTo0KHg+jNnzuihhx7S559/rl69emnEiBH685//HDLHpEmTdPz4cc2ePVsej0eZmZmqrKw860JcAADQPTksy7LC3cSF8vl8crlc8nq9iouLC3c76EY66260nYk73QIwhZ33b35LCAAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMZzuwbNy4URMnTlRycrIcDodWr159zvq3335bt912m/r376+4uDhlZ2frvffeC6l58skn5XA4QpahQ4fabQ0AAEQo24HF7/crIyNDZWVl7arfuHGjbrvtNq1du1a1tbUaP368Jk6cqB07doTUDRs2TA0NDcFl06ZNdlsDAAARqofdDXJzc5Wbm9vu+gULFoQ8/s///E+98847+uMf/6iRI0d+10iPHkpMTLTbDgAA6Aa6/BqWQCCg06dPq2/fviHje/fuVXJysgYOHKif/exnOnToUJtzNDc3y+fzhSwAACBydXlgmTdvnpqamnTXXXcFx7KyslRRUaHKykotXLhQ9fX1uuWWW3T69OlW5ygtLZXL5QouKSkpXdU+AAAIgy4NLK+//rqeeuopvfnmm4qPjw+O5+bm6s4779SIESOUk5OjtWvX6tSpU3rzzTdbnaekpERerze4HD58uKt2AQAAhIHta1g6avny5Zo2bZpWrlwpt9t9ztrevXvrmmuu0b59+1pd73Q65XQ6O6NNAABgoC45w/LGG2+osLBQb7zxhvLy8s5b39TUpP379yspKakLugMAAKazfYalqakp5MxHfX296urq1LdvX1199dUqKSnR559/rmXLlkn65mOggoICvfTSS8rKypLH45Ek9ezZUy6XS5L08MMPa+LEiRowYICOHj2qOXPmKDo6WlOmTPkh9hEAAFzkbJ9h2bZtm0aOHBn8SnJRUZFGjhyp2bNnS5IaGhpCvuHzyiuv6O9//7tmzJihpKSk4DJz5sxgzZEjRzRlyhSlp6frrrvu0hVXXKEtW7aof//+F7p/AAAgAjgsy7LC3cSF8vl8crlc8nq9iouLC3c76EZSi9eEuwXbDsw9/8eyANAV7Lx/81tCAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA824Fl48aNmjhxopKTk+VwOLR69erzblNdXa0f//jHcjqdGjx4sCoqKs6qKSsrU2pqqmJjY5WVlaWtW7fabQ0AAEQo24HF7/crIyNDZWVl7aqvr69XXl6exo8fr7q6Os2aNUvTpk3Te++9F6xZsWKFioqKNGfOHG3fvl0ZGRnKycnRsWPH7LYHAAAikMOyLKvDGzscWrVqlfLz89us+fWvf601a9Zo165dwbHJkyfr1KlTqqyslCRlZWXphhtu0G9/+1tJUiAQUEpKih588EEVFxeftw+fzyeXyyWv16u4uLiO7g5gW2rxmnC3YNuBuXnhbgEAJNl7/+70a1hqamrkdrtDxnJyclRTUyNJOnPmjGpra0NqoqKi5Ha7gzXf19zcLJ/PF7IAAIDI1emBxePxKCEhIWQsISFBPp9Pf/vb33TixAm1tLS0WuPxeFqds7S0VC6XK7ikpKR0Wv8AACD8LspvCZWUlMjr9QaXw4cPh7slAADQiXp09h9ITExUY2NjyFhjY6Pi4uLUs2dPRUdHKzo6utWaxMTEVud0Op1yOp2d1jMAADBLp59hyc7OVlVVVcjY+vXrlZ2dLUmKiYnRqFGjQmoCgYCqqqqCNQAAoHuzHViamppUV1enuro6Sd98bbmurk6HDh2S9M3HNVOnTg3W//KXv9Rnn32mRx99VJ9++qn+67/+S2+++aZ+9atfBWuKior06quvaunSpdq9e7emT58uv9+vwsLCC9w9AAAQCWx/JLRt2zaNHz8++LioqEiSVFBQoIqKCjU0NATDiySlpaVpzZo1+tWvfqWXXnpJV111lV577TXl5OQEayZNmqTjx49r9uzZ8ng8yszMVGVl5VkX4gIAgO7pgu7DYgruw4Jw4T4sANBxRt2HBQAA4EIRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA43UosJSVlSk1NVWxsbHKysrS1q1b26wdN26cHA7HWUteXl6w5p577jlr/YQJEzrSGgAAiEA97G6wYsUKFRUVqby8XFlZWVqwYIFycnK0Z88excfHn1X/9ttv68yZM8HHJ0+eVEZGhu68886QugkTJmjJkiXBx06n025rAAAgQtk+wzJ//nzdd999Kiws1LXXXqvy8nL16tVLixcvbrW+b9++SkxMDC7r169Xr169zgosTqczpK5Pnz4d2yMAABBxbAWWM2fOqLa2Vm63+7sJoqLkdrtVU1PTrjkWLVqkyZMn69JLLw0Zr66uVnx8vNLT0zV9+nSdPHnSTmsAACCC2fpI6MSJE2ppaVFCQkLIeEJCgj799NPzbr9161bt2rVLixYtChmfMGGCfvKTnygtLU379+/XY489ptzcXNXU1Cg6OvqseZqbm9Xc3Bx87PP57OwGAAC4yNi+huVCLFq0SMOHD9fo0aNDxidPnhz8/+HDh2vEiBEaNGiQqqurdeutt541T2lpqZ566qlO7xcAAJjB1kdC/fr1U3R0tBobG0PGGxsblZiYeM5t/X6/li9frnvvvfe8f2fgwIHq16+f9u3b1+r6kpISeb3e4HL48OH27wQAALjo2AosMTExGjVqlKqqqoJjgUBAVVVVys7OPue2K1euVHNzs37+85+f9+8cOXJEJ0+eVFJSUqvrnU6n4uLiQhYAABC5bH9LqKioSK+++qqWLl2q3bt3a/r06fL7/SosLJQkTZ06VSUlJWdtt2jRIuXn5+uKK64IGW9qatIjjzyiLVu26MCBA6qqqtIdd9yhwYMHKycnp4O7BQAAIonta1gmTZqk48ePa/bs2fJ4PMrMzFRlZWXwQtxDhw4pKio0B+3Zs0ebNm3SunXrzpovOjpaH3/8sZYuXapTp04pOTlZt99+u55++mnuxQIAACRJDsuyrHA3caF8Pp9cLpe8Xi8fD6FLpRavCXcLth2Ym3f+IgDoAnbev/ktIQAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwXocCS1lZmVJTUxUbG6usrCxt3bq1zdqKigo5HI6QJTY2NqTGsizNnj1bSUlJ6tmzp9xut/bu3duR1gAAQASyHVhWrFihoqIizZkzR9u3b1dGRoZycnJ07NixNreJi4tTQ0NDcDl48GDI+ueee04vv/yyysvL9dFHH+nSSy9VTk6OvvrqK/t7BAAAIo7twDJ//nzdd999Kiws1LXXXqvy8nL16tVLixcvbnMbh8OhxMTE4JKQkBBcZ1mWFixYoMcff1x33HGHRowYoWXLluno0aNavXp1h3YKAABEFluB5cyZM6qtrZXb7f5ugqgoud1u1dTUtLldU1OTBgwYoJSUFN1xxx365JNPguvq6+vl8XhC5nS5XMrKympzzubmZvl8vpAFAABELluB5cSJE2ppaQk5QyJJCQkJ8ng8rW6Tnp6uxYsX65133tEf/vAHBQIBjRkzRkeOHJGk4HZ25iwtLZXL5QouKSkpdnYDAABcZDr9W0LZ2dmaOnWqMjMzNXbsWL399tvq37+/fve733V4zpKSEnm93uBy+PDhH7BjAABgGluBpV+/foqOjlZjY2PIeGNjoxITE9s1xyWXXKKRI0dq3759khTczs6cTqdTcXFxIQsAAIhctgJLTEyMRo0apaqqquBYIBBQVVWVsrOz2zVHS0uLdu7cqaSkJElSWlqaEhMTQ+b0+Xz66KOP2j0nAACIbD3sblBUVKSCggJdf/31Gj16tBYsWCC/36/CwkJJ0tSpU3XllVeqtLRUkvTv//7vuvHGGzV48GCdOnVKzz//vA4ePKhp06ZJ+uYbRLNmzdIzzzyjIUOGKC0tTU888YSSk5OVn5//w+0pAAC4aNkOLJMmTdLx48c1e/ZseTweZWZmqrKyMnjR7KFDhxQV9d2Jmy+//FL33XefPB6P+vTpo1GjRmnz5s269tprgzWPPvqo/H6/7r//fp06dUo333yzKisrz7rBHAAA6J4clmVZ4W7iQvl8PrlcLnm9Xq5nQZdKLV4T7hZsOzA3L9wtAIAke+/f/JYQAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHg9wt0A0BVSi9eEuwUAwAXgDAsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYr0OBpaysTKmpqYqNjVVWVpa2bt3aZu2rr76qW265RX369FGfPn3kdrvPqr/nnnvkcDhClgkTJnSkNQAAEIFsB5YVK1aoqKhIc+bM0fbt25WRkaGcnBwdO3as1frq6mpNmTJFH3zwgWpqapSSkqLbb79dn3/+eUjdhAkT1NDQEFzeeOONju0RAACIOLYDy/z583XfffepsLBQ1157rcrLy9WrVy8tXry41fr//u//1r/+678qMzNTQ4cO1WuvvaZAIKCqqqqQOqfTqcTExODSp0+fju0RAACIOLYCy5kzZ1RbWyu32/3dBFFRcrvdqqmpadccf/3rX/X111+rb9++IePV1dWKj49Xenq6pk+frpMnT7Y5R3Nzs3w+X8gCAAAil63AcuLECbW0tCghISFkPCEhQR6Pp11z/PrXv1ZycnJI6JkwYYKWLVumqqoqPfvss9qwYYNyc3PV0tLS6hylpaVyuVzBJSUlxc5uAACAi0yPrvxjc+fO1fLly1VdXa3Y2Njg+OTJk4P/P3z4cI0YMUKDBg1SdXW1br311rPmKSkpUVFRUfCxz+cjtAAAEMFsnWHp16+foqOj1djYGDLe2NioxMTEc247b948zZ07V+vWrdOIESPOWTtw4ED169dP+/bta3W90+lUXFxcyAIAACKXrcASExOjUaNGhVww++0FtNnZ2W1u99xzz+npp59WZWWlrr/++vP+nSNHjujkyZNKSkqy0x4AAIhQtr8lVFRUpFdffVVLly7V7t27NX36dPn9fhUWFkqSpk6dqpKSkmD9s88+qyeeeEKLFy9WamqqPB6PPB6PmpqaJElNTU165JFHtGXLFh04cEBVVVW64447NHjwYOXk5PxAuwkAAC5mtq9hmTRpko4fP67Zs2fL4/EoMzNTlZWVwQtxDx06pKio73LQwoULdebMGf30pz8NmWfOnDl68sknFR0drY8//lhLly7VqVOnlJycrNtvv11PP/20nE7nBe4eAACIBA7LsqxwN3GhfD6fXC6XvF4v17OgVanFa8LdgjEOzM0LdwsAIMne+ze/JQQAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMF6X/vghcC7cKwUA0BYCC9DNdGYw5KZ0ADoLHwkBAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPG4NT9s4zd/AABdjTMsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMbrUGApKytTamqqYmNjlZWVpa1bt56zfuXKlRo6dKhiY2M1fPhwrV27NmS9ZVmaPXu2kpKS1LNnT7ndbu3du7cjrQEAgAjUw+4GK1asUFFRkcrLy5WVlaUFCxYoJydHe/bsUXx8/Fn1mzdv1pQpU1RaWqp/+qd/0uuvv678/Hxt375d1113nSTpueee08svv6ylS5cqLS1NTzzxhHJycvSXv/xFsbGxF76X3VBq8Zpwt4BuqLOedwfm5nXKvAAuHg7Lsiw7G2RlZemGG27Qb3/7W0lSIBBQSkqKHnzwQRUXF59VP2nSJPn9fr377rvBsRtvvFGZmZkqLy+XZVlKTk7WQw89pIcffliS5PV6lZCQoIqKCk2ePPm8Pfl8PrlcLnm9XsXFxdnZnYhFYEEkIbAAkcnO+7etMyxnzpxRbW2tSkpKgmNRUVFyu92qqalpdZuamhoVFRWFjOXk5Gj16tWSpPr6enk8Hrnd7uB6l8ulrKws1dTUtBpYmpub1dzcHHzs9XolfbPj+Eag+a/hbgH4wfBvG4hM3/7bbs+5E1uB5cSJE2ppaVFCQkLIeEJCgj799NNWt/F4PK3Wezye4Ppvx9qq+b7S0lI99dRTZ42npKS0b0cAXFRcC8LdAYDOdPr0ablcrnPW2L6GxQQlJSUhZ20CgYC++OILXXHFFXI4HD/o3/L5fEpJSdHhw4f5uOk8OFbtx7FqP46VPRyv9uNYtV9nHSvLsnT69GklJyeft9ZWYOnXr5+io6PV2NgYMt7Y2KjExMRWt0lMTDxn/bf/bWxsVFJSUkhNZmZmq3M6nU45nc6Qsd69e9vZFdvi4uJ4QrcTx6r9OFbtx7Gyh+PVfhyr9uuMY3W+MyvfsvW15piYGI0aNUpVVVXBsUAgoKqqKmVnZ7e6TXZ2dki9JK1fvz5Yn5aWpsTExJAan8+njz76qM05AQBA92L7I6GioiIVFBTo+uuv1+jRo7VgwQL5/X4VFhZKkqZOnaorr7xSpaWlkqSZM2dq7NixeuGFF5SXl6fly5dr27ZteuWVVyRJDodDs2bN0jPPPKMhQ4YEv9acnJys/Pz8H25PAQDARct2YJk0aZKOHz+u2bNny+PxKDMzU5WVlcGLZg8dOqSoqO9O3IwZM0avv/66Hn/8cT322GMaMmSIVq9eHbwHiyQ9+uij8vv9uv/++3Xq1CndfPPNqqysNOIeLE6nU3PmzDnrIyicjWPVfhyr9uNY2cPxaj+OVfuZcKxs34cFAACgq/FbQgAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AIqmsrEypqamKjY1VVlaWtm7des76lStXaujQoYqNjdXw4cO1du3aLuo0/Owcq4qKCjkcjpDFhG9+dYWNGzdq4sSJSk5OlsPhCP521rlUV1frxz/+sZxOpwYPHqyKiopO79MEdo9VdXX1Wc8rh8PR5k95RJLS0lLdcMMNuvzyyxUfH6/8/Hzt2bPnvNt1x9esjhyr7vqatXDhQo0YMSJ4U7js7Gz96U9/Ouc24XhOdfvAsmLFChUVFWnOnDnavn27MjIylJOTo2PHjrVav3nzZk2ZMkX33nuvduzYofz8fOXn52vXrl1d3HnXs3uspG/uitjQ0BBcDh482IUdh4/f71dGRobKysraVV9fX6+8vDyNHz9edXV1mjVrlqZNm6b33nuvkzsNP7vH6lt79uwJeW7Fx8d3Uofm2LBhg2bMmKEtW7Zo/fr1+vrrr3X77bfL7/e3uU13fc3qyLGSuudr1lVXXaW5c+eqtrZW27Zt0z/8wz/ojjvu0CeffNJqfdieU1Y3N3r0aGvGjBnBxy0tLVZycrJVWlraav1dd91l5eXlhYxlZWVZ//Iv/9KpfZrA7rFasmSJ5XK5uqg7c0myVq1adc6aRx991Bo2bFjI2KRJk6ycnJxO7Mw87TlWH3zwgSXJ+vLLL7ukJ5MdO3bMkmRt2LChzZru/Jr1f7XnWPGa9Z0+ffpYr732WqvrwvWc6tZnWM6cOaPa2lq53e7gWFRUlNxut2pqalrdpqamJqReknJyctqsjxQdOVaS1NTUpAEDBiglJeWcib27667PqwuRmZmppKQk3Xbbbfrwww/D3U5YeL1eSVLfvn3brOG59Y32HCuJ16yWlhYtX75cfr+/zZ/HCddzqlsHlhMnTqilpSV4l95vJSQktPl5uMfjsVUfKTpyrNLT07V48WK98847+sMf/qBAIKAxY8boyJEjXdHyRaWt55XP59Pf/va3MHVlpqSkJJWXl+utt97SW2+9pZSUFI0bN07bt28Pd2tdKhAIaNasWbrppptC7hz+fd31Nev/au+x6s6vWTt37tRll10mp9OpX/7yl1q1apWuvfbaVmvD9ZyyfWt+oL2ys7NDEvqYMWP0ox/9SL/73e/09NNPh7EzXMzS09OVnp4efDxmzBjt379fL774on7/+9+HsbOuNWPGDO3atUubNm0KdyvGa++x6s6vWenp6aqrq5PX69X//M//qKCgQBs2bGgztIRDtz7D0q9fP0VHR6uxsTFkvLGxUYmJia1uk5iYaKs+UnTkWH3fJZdcopEjR2rfvn2d0eJFra3nVVxcnHr27Bmmri4eo0eP7lbPqwceeEDvvvuuPvjgA1111VXnrO2ur1nfsnOsvq87vWbFxMRo8ODBGjVqlEpLS5WRkaGXXnqp1dpwPae6dWCJiYnRqFGjVFVVFRwLBAKqqqpq87O77OzskHpJWr9+fZv1kaIjx+r7WlpatHPnTiUlJXVWmxet7vq8+qHU1dV1i+eVZVl64IEHtGrVKr3//vtKS0s77zbd9bnVkWP1fd35NSsQCKi5ubnVdWF7TnXqJb0XgeXLl1tOp9OqqKiw/vKXv1j333+/1bt3b8vj8ViWZVl33323VVxcHKz/8MMPrR49eljz5s2zdu/ebc2ZM8e65JJLrJ07d4ZrF7qM3WP11FNPWe+99561f/9+q7a21po8ebIVGxtrffLJJ+HahS5z+vRpa8eOHdaOHTssSdb8+fOtHTt2WAcPHrQsy7KKi4utu+++O1j/2WefWb169bIeeeQRa/fu3VZZWZkVHR1tVVZWhmsXuozdY/Xiiy9aq1evtvbu3Wvt3LnTmjlzphUVFWX9+c9/DtcudJnp06dbLpfLqq6uthoaGoLLX//612ANr1nf6Mix6q6vWcXFxdaGDRus+vp66+OPP7aKi4sth8NhrVu3zrIsc55T3T6wWJZl/eY3v7GuvvpqKyYmxho9erS1ZcuW4LqxY8daBQUFIfVvvvmmdc0111gxMTHWsGHDrDVr1nRxx+Fj51jNmjUrWJuQkGD94z/+o7V9+/YwdN31vv3q7feXb49PQUGBNXbs2LO2yczMtGJiYqyBAwdaS5Ys6fK+w8HusXr22WetQYMGWbGxsVbfvn2tcePGWe+//354mu9irR0nSSHPFV6zvtGRY9VdX7N+8YtfWAMGDLBiYmKs/v37W7feemswrFiWOc8ph2VZVueewwEAALgw3foaFgAAcHEgsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeP8f72oB3J5OJH8AAAAASUVORK5CYII=\n", "text/plain": [ - "<Figure size 432x288 with 1 Axes>" + "<Figure size 640x480 with 1 Axes>" ] }, - "metadata": { - "needs_background": "light" - }, + "metadata": {}, "output_type": "display_data" } ], @@ -160,35 +156,38 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 9, "id": "4a38abae", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "(array([ 0., 0., 0., 0., 9., 169., 704., 118., 34., 25., 19.,\n", - " 9., 4., 4., 2., 0., 1., 0., 0., 0.]),\n", - " array([0. , 0.15, 0.3 , 0.45, 0.6 , 0.75, 0.9 , 1.05, 1.2 , 1.35, 1.5 ,\n", - " 1.65, 1.8 , 1.95, 2.1 , 2.25, 2.4 , 2.55, 2.7 , 2.85, 3. ],\n", - " dtype=float32),\n", + "(array([0.000000e+00, 0.000000e+00, 0.000000e+00, 2.050000e+02,\n", + " 2.136400e+04, 7.107010e+05, 2.630348e+06, 5.204000e+05,\n", + " 1.720750e+05, 9.079700e+04, 6.480300e+04, 4.666500e+04,\n", + " 3.338600e+04, 2.284300e+04, 1.320600e+04, 8.632000e+03,\n", + " 4.469000e+03, 1.948000e+03, 1.068000e+03, 5.150000e+02]),\n", + " array([0. , 0.15000001, 0.30000001, 0.44999999, 0.60000002,\n", + " 0.75 , 0.89999998, 1.04999995, 1.20000005, 1.35000002,\n", + " 1.5 , 1.64999998, 1.79999995, 1.95000005, 2.0999999 ,\n", + " 2.25 , 2.4000001 , 2.54999995, 2.70000005, 2.8499999 ,\n", + " 3. ]),\n", " <BarContainer object of 20 artists>)" ] }, - "execution_count": 10, + "execution_count": 9, "metadata": {}, "output_type": "execute_result" }, { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAARK0lEQVR4nO3dbYxcV33H8e+POAEEFOdh61q2i6mwQGlVglmlRlSIYlHlocKRClFQRUzkylWbtiAqFZcXRVR9Ed5ASVsFWYTWqSiQBmjcEGgtJwj1RQIbCAESaJaIyLaSeAnEPKSAQv99sccwWdae2d3ZXe/h+5FGc+655+78T67z890zd8apKiRJfXnGahcgSRo/w12SOmS4S1KHDHdJ6pDhLkkdWrfaBQBccMEFtXXr1tUuQ5LWlHvuuedbVTUx374zIty3bt3K1NTUapchSWtKkodPtc9lGUnqkOEuSR0y3CWpQ4a7JHXIcJekDhnuktQhw12SOjQ03JO8OMm9A4/vJnlrkvOSHEryYHs+t41PkuuTTCe5L8n25Z+GJGnQ0HCvqq9X1UVVdRHwcuBJ4BPAPuBwVW0DDrdtgEuBbe2xF7hhGeqWJJ3GQj+huhP4RlU9nGQX8OrWfwD4DPB2YBdwU83+KyB3JVmfZGNVPTKmmrWKtu775KKP/eZ1l4+xEkmns9A196uAD7f2hoHAfhTY0NqbgCMDxxxtfU+TZG+SqSRTMzMzCyxDknQ6I4d7knOA1wH/Nndfu0pf0L/XV1X7q2qyqiYnJub93htJ0iIt5Mr9UuALVfVY234syUaA9ny89R8Dtgwct7n1SZJWyELC/Y38bEkG4CCwu7V3A7cO9F/d7prZAZxwvV2SVtZIb6gmeQ7wWuCPBrqvA25Osgd4GLiy9d8OXAZMM3tnzTVjq1aSNJKRwr2qfgCcP6fvcWbvnpk7toBrx1KdJGlR/ISqJHXIcJekDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6NFK4J1mf5JYkX0vyQJJXJDkvyaEkD7bnc9vYJLk+yXSS+5JsX94pSJLmGvXK/X3Ap6vqJcBLgQeAfcDhqtoGHG7bAJcC29pjL3DDWCuWJA01NNyTPB94FXAjQFX9uKqeAHYBB9qwA8AVrb0LuKlm3QWsT7JxzHVLkk5jlCv3FwIzwD8l+WKSDyR5DrChqh5pYx4FNrT2JuDIwPFHW9/TJNmbZCrJ1MzMzOJnIEn6OaOE+zpgO3BDVb0M+AE/W4IBoKoKqIW8cFXtr6rJqpqcmJhYyKGSpCFGCfejwNGqurtt38Js2D92crmlPR9v+48BWwaO39z6JEkrZGi4V9WjwJEkL25dO4H7gYPA7ta3G7i1tQ8CV7e7ZnYAJwaWbyRJK2DdiOP+DPhQknOAh4BrmP2L4eYke4CHgSvb2NuBy4Bp4Mk2VpK0gkYK96q6F5icZ9fOecYWcO3SypIkLYWfUJWkDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nq0EjhnuSbSb6c5N4kU63vvCSHkjzYns9t/UlyfZLpJPcl2b6cE5Ak/byFXLn/TlVdVFWTbXsfcLiqtgGH2zbApcC29tgL3DCuYiVJo1nKsswu4EBrHwCuGOi/qWbdBaxPsnEJryNJWqBRw72A/0pyT5K9rW9DVT3S2o8CG1p7E3Bk4Nijre9pkuxNMpVkamZmZhGlS5JOZd2I4367qo4l+WXgUJKvDe6sqkpSC3nhqtoP7AeYnJxc0LGSpNMb6cq9qo615+PAJ4CLgcdOLre05+Nt+DFgy8Dhm1ufJGmFDA33JM9J8ryTbeB3ga8AB4Hdbdhu4NbWPghc3e6a2QGcGFi+kSStgFGWZTYAn0hycvy/VtWnk3weuDnJHuBh4Mo2/nbgMmAaeBK4ZuxVS5JOa2i4V9VDwEvn6X8c2DlPfwHXjqU6SdKi+AlVSeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR0y3CWpQ4a7JHXIcJekDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1aORwT3JWki8mua1tvzDJ3Ummk3w0yTmt/5lte7rt37pMtUuSTmEhV+5vAR4Y2H438N6qehHwHWBP698DfKf1v7eNkyStoJHCPclm4HLgA207wGuAW9qQA8AVrb2rbdP272zjJUkrZNQr978D/hL4v7Z9PvBEVT3Vto8Cm1p7E3AEoO0/0cY/TZK9SaaSTM3MzCyueknSvIaGe5LfA45X1T3jfOGq2l9Vk1U1OTExMc4fLUm/8NaNMOaVwOuSXAY8C/gl4H3A+iTr2tX5ZuBYG38M2AIcTbIOeD7w+NgrlySd0tAr96r6q6raXFVbgauAO6rqD4A7gde3YbuBW1v7YNum7b+jqmqsVUuSTmsp97m/HXhbkmlm19RvbP03Aue3/rcB+5ZWoiRpoUZZlvmpqvoM8JnWfgi4eJ4xPwTeMIbaJEmL5CdUJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR0y3CWpQ4a7JHXIcJekDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUoaHhnuRZST6X5EtJvprkXa3/hUnuTjKd5KNJzmn9z2zb023/1mWegyRpjlGu3H8EvKaqXgpcBFySZAfwbuC9VfUi4DvAnjZ+D/Cd1v/eNk6StIKGhnvN+n7bPLs9CngNcEvrPwBc0dq72jZt/84kGVfBkqThRlpzT3JWknuB48Ah4BvAE1X1VBtyFNjU2puAIwBt/wng/DHWLEkaYqRwr6qfVNVFwGbgYuAlS33hJHuTTCWZmpmZWeqPkyQNWNDdMlX1BHAn8ApgfZJ1bddm4FhrHwO2ALT9zwcen+dn7a+qyaqanJiYWFz1kqR5jXK3zESS9a39bOC1wAPMhvzr27DdwK2tfbBt0/bfUVU1xpolSUOsGz6EjcCBJGcx+5fBzVV1W5L7gY8k+Vvgi8CNbfyNwL8kmQa+DVy1DHVLkk5jaLhX1X3Ay+bpf4jZ9fe5/T8E3jCW6iRJi+InVCWpQ4a7JHXIcJekDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6NDTck2xJcmeS+5N8NclbWv95SQ4lebA9n9v6k+T6JNNJ7kuyfbknIUl6ulGu3J8C/qKqLgR2ANcmuRDYBxyuqm3A4bYNcCmwrT32AjeMvWpJ0mkNDfeqeqSqvtDa3wMeADYBu4ADbdgB4IrW3gXcVLPuAtYn2TjuwiVJp7agNfckW4GXAXcDG6rqkbbrUWBDa28CjgwcdrT1zf1Ze5NMJZmamZlZaN2SpNMYOdyTPBf4GPDWqvru4L6qKqAW8sJVtb+qJqtqcmJiYiGHSpKGGCnck5zNbLB/qKo+3rofO7nc0p6Pt/5jwJaBwze3PknSChnlbpkANwIPVNV7BnYdBHa39m7g1oH+q9tdMzuAEwPLN5KkFbBuhDGvBN4EfDnJva3vHcB1wM1J9gAPA1e2fbcDlwHTwJPANeMsWJI03NBwr6r/BnKK3TvnGV/AtUusS5K0BH5CVZI6ZLhLUocMd0nq0ChvqKojW/d9crVLkLQCvHKXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR3yK3+1YpbydcPfvO7yMVYi9c8rd0nq0NBwT/LBJMeTfGWg77wkh5I82J7Pbf1Jcn2S6ST3Jdm+nMVLkuY3ypX7PwOXzOnbBxyuqm3A4bYNcCmwrT32AjeMp0xJ0kIMDfeq+izw7Tndu4ADrX0AuGKg/6aadRewPsnGMdUqSRrRYtfcN1TVI639KLChtTcBRwbGHW19PyfJ3iRTSaZmZmYWWYYkaT5LfkO1qgqoRRy3v6omq2pyYmJiqWVIkgYsNtwfO7nc0p6Pt/5jwJaBcZtbnyRpBS023A8Cu1t7N3DrQP/V7a6ZHcCJgeUbSdIKGfohpiQfBl4NXJDkKPBO4Drg5iR7gIeBK9vw24HLgGngSeCaZahZkjTE0HCvqjeeYtfOecYWcO1Si5IkLY2fUJWkDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHRr6lb/SmWDrvk8u+thvXnf5GCuR1gav3CWpQ165q3te9esXkVfuktQhr9yl01jKVT945a/V45W7JHXIcJekDi3LskySS4D3AWcBH6iq65bjdX5RLXWpQCvHN3O1WsZ+5Z7kLOAfgUuBC4E3Jrlw3K8jSTq15bhyvxiYrqqHAJJ8BNgF3L8MryV1a63+huZvHGeG5Qj3TcCRge2jwG/NHZRkL7C3bX4/ydcX+XoXAN9a5LFnGudy5ullHrBCc8m7l/sVAM/LSS841Y5VuxWyqvYD+5f6c5JMVdXkGEpadc7lzNPLPMC5nKmWay7LcbfMMWDLwPbm1idJWiHLEe6fB7YleWGSc4CrgIPL8DqSpFMY+7JMVT2V5E+B/2T2VsgPVtVXx/06A5a8tHMGcS5nnl7mAc7lTLUsc0lVLcfPlSStIj+hKkkdMtwlqUNrJtyTXJLk60mmk+ybZ/8zk3y07b87ydZVKHMkI8zlzUlmktzbHn+4GnUOk+SDSY4n+cop9ifJ9W2e9yXZvtI1jmqEubw6yYmBc/LXK13jKJJsSXJnkvuTfDXJW+YZsybOy4hzWSvn5VlJPpfkS20u75pnzHgzrKrO+Aezb8x+A/g14BzgS8CFc8b8CfD+1r4K+Ohq172EubwZ+IfVrnWEubwK2A585RT7LwM+BQTYAdy92jUvYS6vBm5b7TpHmMdGYHtrPw/4n3n+fK2J8zLiXNbKeQnw3NY+G7gb2DFnzFgzbK1cuf/0Kw2q6sfAya80GLQLONDatwA7k2QFaxzVKHNZE6rqs8C3TzNkF3BTzboLWJ9k48pUtzAjzGVNqKpHquoLrf094AFmPzU+aE2clxHnsia0/9bfb5tnt8fcu1nGmmFrJdzn+0qDuSf5p2Oq6ingBHD+ilS3MKPMBeD326/MtyTZMs/+tWDUua4Vr2i/Vn8qya+vdjHDtF/rX8bsVeKgNXdeTjMXWCPnJclZSe4FjgOHquqU52UcGbZWwv0XzX8AW6vqN4FD/Oxvc62eLwAvqKqXAn8P/PvqlnN6SZ4LfAx4a1V9d7XrWYohc1kz56WqflJVFzH7qf2Lk/zGcr7eWgn3Ub7S4KdjkqwDng88viLVLczQuVTV41X1o7b5AeDlK1TbuHXzVRRV9d2Tv1ZX1e3A2UkuWOWy5pXkbGbD8ENV9fF5hqyZ8zJsLmvpvJxUVU8AdwKXzNk11gxbK+E+ylcaHAR2t/brgTuqvTNxhhk6lznrn69jdq1xLToIXN3uztgBnKiqR1a7qMVI8isn1z+TXMzs/ztn3MVDq/FG4IGqes8phq2J8zLKXNbQeZlIsr61nw28FvjanGFjzbA18Q9k1ym+0iDJ3wBTVXWQ2T8E/5Jkmtk3xq5avYpPbcS5/HmS1wFPMTuXN69awaeR5MPM3q1wQZKjwDuZfaOIqno/cDuzd2ZMA08C16xOpcONMJfXA3+c5Cngf4GrztCLh1cCbwK+3NZ3Ad4B/CqsufMyylzWynnZCBzI7D9m9Azg5qq6bTkzzK8fkKQOrZVlGUnSAhjuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUP/D2YsydU15EV0AAAAAElFTkSuQmCC\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGsCAYAAAAPJKchAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAfIklEQVR4nO3de2yV533A8d8BgkkX7ISm2AYcoKGhuRpKuZisBTYaj7Io1rSMRVphbcnWCqpQtlWwS6K0m5woIWHqWGiUBSuNGCnNAA3aBAoFlOCo4qZClrKSC5AEO4mS2MHtHGa/+yOKOzeY+DiYx5fPR3r/8OvnPec5j46Ov3rPe45zWZZlAQCQyIDUEwAA+jcxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJNWrYmT37t1x4403xogRIyKXy8XGjRvzvo0sy+Lee++NK664IgoKCmLkyJHxT//0T+d+sgBApwxKPYF8NDU1RXl5eXzlK1+JP/qjP+rSbdx2222xdevWuPfee+Paa6+NN998M958881zPFMAoLNyvfUf5eVyudiwYUNUVVW17Wtubo6/+7u/i3//93+Pt99+O6655pq4++67Y+bMmRER8dxzz8V1110Xhw8fjvHjx6eZOADQTq96m+bDLF68OGpra2PdunXx85//PG6++eb4gz/4g/jlL38ZERH/+Z//GZ/85Cdj8+bNMXbs2BgzZkwsXLjQmREASKjPxMjx48djzZo1sX79+vjc5z4Xl19+efz1X/91/O7v/m6sWbMmIiJeeOGFOHbsWKxfvz4eeeSRqKmpiX379sUf//EfJ549APRfveqakbM5dOhQtLS0xBVXXNFuf3Nzc3z84x+PiIjW1tZobm6ORx55pG3cv/3bv8WkSZPiyJEj3roBgAT6TIycOnUqBg4cGPv27YuBAwe2+91FF10UERGlpaUxaNCgdsFy5ZVXRsR7Z1bECACcf30mRiZOnBgtLS3x2muvxec+97kzjrn++uvjf//3f+P555+Pyy+/PCIi/vu//zsiIkaPHn3e5goA/Eav+jTNqVOn4ujRoxHxXnzcd999MWvWrBg2bFhcdtll8Wd/9mfx9NNPx4oVK2LixInx+uuvx/bt2+O6666LuXPnRmtra0yePDkuuuiiWLlyZbS2tsaiRYuisLAwtm7dmvjRAUD/1KtiZOfOnTFr1qwP7F+wYEHU1NTE6dOn4x//8R/jkUceiVdeeSUuvfTSmDZtWtx5551x7bXXRkTEq6++Gt/4xjdi69at8Tu/8zsxZ86cWLFiRQwbNux8PxwAIHpZjAAAfU+f+WgvANA7iREAIKle8Wma1tbWePXVV2Po0KGRy+VSTwcA6IQsy+Kdd96JESNGxIABHZ//6BUx8uqrr0ZZWVnqaQAAXXDixIkYNWpUh7/vFTEydOjQiHjvwRQWFiaeDQDQGY2NjVFWVtb2d7wjvSJG3n9rprCwUIwAQC/zYZdYuIAVAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJDUoNQTgJ5szLIt3XK7L901t1tuF6A3cmYEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApPKKkerq6pg8eXIMHTo0hg8fHlVVVXHkyJGzHlNTUxO5XK7dNmTIkI80aQCg78grRnbt2hWLFi2KZ555JrZt2xanT5+OG264IZqams56XGFhYZw8ebJtO3bs2EeaNADQdwzKZ/ATTzzR7ueampoYPnx47Nu3Lz7/+c93eFwul4uSkpKuzRAA6NM+0jUjDQ0NERExbNiws447depUjB49OsrKyuKmm26KZ5999qzjm5ubo7Gxsd0GAPRNXY6R1tbWWLJkSVx//fVxzTXXdDhu/Pjx8fDDD8emTZvi0UcfjdbW1pg+fXq8/PLLHR5TXV0dRUVFbVtZWVlXpwkA9HC5LMuyrhz49a9/PX784x/HU089FaNGjer0cadPn44rr7wybrnllvjOd75zxjHNzc3R3Nzc9nNjY2OUlZVFQ0NDFBYWdmW60CVjlm3pltt96a653XK7AD1JY2NjFBUVfejf77yuGXnf4sWLY/PmzbF79+68QiQi4oILLoiJEyfG0aNHOxxTUFAQBQUFXZkaANDL5PU2TZZlsXjx4tiwYUPs2LEjxo4dm/cdtrS0xKFDh6K0tDTvYwGAvievMyOLFi2KtWvXxqZNm2Lo0KFRV1cXERFFRUVx4YUXRkTE/PnzY+TIkVFdXR0REd/+9rdj2rRpMW7cuHj77bfjnnvuiWPHjsXChQvP8UMBAHqjvGLkgQceiIiImTNnttu/Zs2a+PM///OIiDh+/HgMGPCbEy5vvfVW3HrrrVFXVxeXXHJJTJo0Kfbs2RNXXXXVR5s5ANAndPkC1vOpsxfAwLnmAlaAruvs32//mwYASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUnnFSHV1dUyePDmGDh0aw4cPj6qqqjhy5MiHHrd+/fr49Kc/HUOGDIlrr702fvSjH3V5wgBA35JXjOzatSsWLVoUzzzzTGzbti1Onz4dN9xwQzQ1NXV4zJ49e+KWW26Jr371q3HgwIGoqqqKqqqqOHz48EeePADQ++WyLMu6evDrr78ew4cPj127dsXnP//5M46ZN29eNDU1xebNm9v2TZs2LSZMmBCrV6/u1P00NjZGUVFRNDQ0RGFhYVenC3kbs2xLt9zuS3fN7ZbbBehJOvv3+yNdM9LQ0BAREcOGDetwTG1tbcyePbvdvsrKyqitre3wmObm5mhsbGy3AQB9U5djpLW1NZYsWRLXX399XHPNNR2Oq6uri+Li4nb7iouLo66ursNjqquro6ioqG0rKyvr6jQBgB6uyzGyaNGiOHz4cKxbt+5cziciIpYvXx4NDQ1t24kTJ875fQAAPcOgrhy0ePHi2Lx5c+zevTtGjRp11rElJSVRX1/fbl99fX2UlJR0eExBQUEUFBR0ZWoAQC+T15mRLMti8eLFsWHDhtixY0eMHTv2Q4+pqKiI7du3t9u3bdu2qKioyG+mAECflNeZkUWLFsXatWtj06ZNMXTo0LbrPoqKiuLCCy+MiIj58+fHyJEjo7q6OiIibrvttpgxY0asWLEi5s6dG+vWrYu9e/fGgw8+eI4fCgDQG+V1ZuSBBx6IhoaGmDlzZpSWlrZtjz32WNuY48ePx8mTJ9t+nj59eqxduzYefPDBKC8vjx/+8IexcePGs170CgD0H3mdGenMV5Ls3LnzA/tuvvnmuPnmm/O5KwCgn/C/aQCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAklXeM7N69O2688cYYMWJE5HK52Lhx41nH79y5M3K53Ae2urq6rs4ZAOhD8o6RpqamKC8vj1WrVuV13JEjR+LkyZNt2/Dhw/O9awCgDxqU7wFz5syJOXPm5H1Hw4cPj4svvjjv4wCAvu28XTMyYcKEKC0tjS984Qvx9NNPn3Vsc3NzNDY2ttsAgL6p22OktLQ0Vq9eHY8//ng8/vjjUVZWFjNnzoz9+/d3eEx1dXUUFRW1bWVlZd09TQAgkVyWZVmXD87lYsOGDVFVVZXXcTNmzIjLLrssvv/975/x983NzdHc3Nz2c2NjY5SVlUVDQ0MUFhZ2dbqQtzHLtnTL7b5019xuuV2AnqSxsTGKioo+9O933teMnAtTpkyJp556qsPfFxQUREFBwXmcEQCQSpLvGTl48GCUlpamuGsAoIfJ+8zIqVOn4ujRo20/v/jii3Hw4MEYNmxYXHbZZbF8+fJ45ZVX4pFHHomIiJUrV8bYsWPj6quvjv/5n/+Jhx56KHbs2BFbt249d48CAOi18o6RvXv3xqxZs9p+Xrp0aURELFiwIGpqauLkyZNx/Pjxtt+/++678Vd/9VfxyiuvxMc+9rG47rrr4ic/+Um72wAA+q+PdAHr+dLZC2DgXHMBK0DXdfbvt/9NAwAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASeUdI7t3744bb7wxRowYEblcLjZu3Pihx+zcuTM+85nPREFBQYwbNy5qamq6MFUAoC/KO0aampqivLw8Vq1a1anxL774YsydOzdmzZoVBw8ejCVLlsTChQvjySefzHuyAEDfMyjfA+bMmRNz5szp9PjVq1fH2LFjY8WKFRERceWVV8ZTTz0V999/f1RWVuZ79wBAH9Pt14zU1tbG7Nmz2+2rrKyM2traDo9pbm6OxsbGdhsA0Dd1e4zU1dVFcXFxu33FxcXR2NgYv/71r894THV1dRQVFbVtZWVl3T1NACCRHvlpmuXLl0dDQ0PbduLEidRTAgC6Sd7XjOSrpKQk6uvr2+2rr6+PwsLCuPDCC894TEFBQRQUFHT31ACAHqDbz4xUVFTE9u3b2+3btm1bVFRUdPddAwC9QN4xcurUqTh48GAcPHgwIt776O7Bgwfj+PHjEfHeWyzz589vG/+1r30tXnjhhfjWt74Vv/jFL+Jf//Vf4wc/+EF885vfPDePAADo1fKOkb1798bEiRNj4sSJERGxdOnSmDhxYtx+++0REXHy5Mm2MImIGDt2bGzZsiW2bdsW5eXlsWLFinjooYd8rBcAiIiIXJZlWepJfJjGxsYoKiqKhoaGKCwsTD0d+pExy7Z0y+2+dNfcbrldgJ6ks3+/e+SnaQCA/kOMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUoNSTwA+qjHLtqSeAgAfgTMjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJBUl2Jk1apVMWbMmBgyZEhMnTo1fvazn3U4tqamJnK5XLttyJAhXZ4wANC35B0jjz32WCxdujTuuOOO2L9/f5SXl0dlZWW89tprHR5TWFgYJ0+ebNuOHTv2kSYNAPQdg/I94L777otbb701vvzlL0dExOrVq2PLli3x8MMPx7Jly854TC6Xi5KSko82U+hDxizb0m23/dJdc7vttgG6Q15nRt59993Yt29fzJ49+zc3MGBAzJ49O2prazs87tSpUzF69OgoKyuLm266KZ599tmz3k9zc3M0Nja22wCAvimvGHnjjTeipaUliouL2+0vLi6Ourq6Mx4zfvz4ePjhh2PTpk3x6KOPRmtra0yfPj1efvnlDu+nuro6ioqK2raysrJ8pgkA9CLd/mmaioqKmD9/fkyYMCFmzJgR//Ef/xGf+MQn4nvf+16HxyxfvjwaGhrathMnTnT3NAGARPK6ZuTSSy+NgQMHRn19fbv99fX1nb4m5IILLoiJEyfG0aNHOxxTUFAQBQUF+UwNAOil8jozMnjw4Jg0aVJs3769bV9ra2ts3749KioqOnUbLS0tcejQoSgtLc1vpgBAn5T3p2mWLl0aCxYsiM9+9rMxZcqUWLlyZTQ1NbV9umb+/PkxcuTIqK6ujoiIb3/72zFt2rQYN25cvP3223HPPffEsWPHYuHChef2kQAAvVLeMTJv3rx4/fXX4/bbb4+6urqYMGFCPPHEE20XtR4/fjwGDPjNCZe33norbr311qirq4tLLrkkJk2aFHv27Imrrrrq3D0KAKDXymVZlqWexIdpbGyMoqKiaGhoiMLCwtTToYfpzu/s6I18zwjQU3T277f/TQMAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSg1JPADi3xizb0i23+9Jdc7vldgGcGQEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUr4OHuiU7vqa+QhfNQ/9nTMjAEBSYgQASEqMAABJuWYESK67rkdxLQr0Ds6MAABJOTMC9Fk+AQS9gzMjAEBSXTozsmrVqrjnnnuirq4uysvL47vf/W5MmTKlw/Hr16+Pf/iHf4iXXnopPvWpT8Xdd98dX/ziF7s8aYDUXOcC507eMfLYY4/F0qVLY/Xq1TF16tRYuXJlVFZWxpEjR2L48OEfGL9nz5645ZZborq6Ov7wD/8w1q5dG1VVVbF///645pprzsmDoOfrztPl0Jd4a4n+KJdlWZbPAVOnTo3JkyfHv/zLv0RERGtra5SVlcU3vvGNWLZs2QfGz5s3L5qammLz5s1t+6ZNmxYTJkyI1atXd+o+Gxsbo6ioKBoaGqKwsDCf6dJDiBFIT4xwvnX273deZ0befffd2LdvXyxfvrxt34ABA2L27NlRW1t7xmNqa2tj6dKl7fZVVlbGxo0bO7yf5ubmaG5ubvu5oaEhIt57UPROrc2/Sj0F6Pcu++b61FPI2+E7K1NPgY/g/b/bH3beI68YeeONN6KlpSWKi4vb7S8uLo5f/OIXZzymrq7ujOPr6uo6vJ/q6uq48847P7C/rKwsn+kC0MsVrUw9A86Fd955J4qKijr8fY/8aO/y5cvbnU1pbW2NN998Mz7+8Y9HLpc7Z/fT2NgYZWVlceLECW//fAhrlR/r1XnWqvOsVedZq87rzrXKsizeeeedGDFixFnH5RUjl156aQwcODDq6+vb7a+vr4+SkpIzHlNSUpLX+IiIgoKCKCgoaLfv4osvzmeqeSksLPRk7SRrlR/r1XnWqvOsVedZq87rrrU62xmR9+X1PSODBw+OSZMmxfbt29v2tba2xvbt26OiouKMx1RUVLQbHxGxbdu2DscDAP1L3m/TLF26NBYsWBCf/exnY8qUKbFy5cpoamqKL3/5yxERMX/+/Bg5cmRUV1dHRMRtt90WM2bMiBUrVsTcuXNj3bp1sXfv3njwwQfP7SMBAHqlvGNk3rx58frrr8ftt98edXV1MWHChHjiiSfaLlI9fvx4DBjwmxMu06dPj7Vr18bf//3fx9/+7d/Gpz71qdi4cWOP+I6RgoKCuOOOOz7wlhAfZK3yY706z1p1nrXqPGvVeT1hrfL+nhEAgHPJ/6YBAJISIwBAUmIEAEhKjAAASfX5GFm1alWMGTMmhgwZElOnTo2f/exnZx2/fv36+PSnPx1DhgyJa6+9Nn70ox+dp5mml89a1dTURC6Xa7cNGTLkPM42nd27d8eNN94YI0aMiFwud9b/s/S+nTt3xmc+85koKCiIcePGRU1NTbfPsyfId6127tz5gedVLpc767+P6Cuqq6tj8uTJMXTo0Bg+fHhUVVXFkSNHPvS4/via1ZW16q+vWQ888EBcd911bV9oVlFRET/+8Y/PekyK51SfjpHHHnssli5dGnfccUfs378/ysvLo7KyMl577bUzjt+zZ0/ccsst8dWvfjUOHDgQVVVVUVVVFYcPHz7PMz//8l2riPe+re/kyZNt27Fjx87jjNNpamqK8vLyWLVqVafGv/jiizF37tyYNWtWHDx4MJYsWRILFy6MJ598sptnml6+a/W+I0eOtHtuDR8+vJtm2HPs2rUrFi1aFM8880xs27YtTp8+HTfccEM0NTV1eEx/fc3qylpF9M/XrFGjRsVdd90V+/bti71798bv/d7vxU033RTPPvvsGccne05lfdiUKVOyRYsWtf3c0tKSjRgxIquurj7j+D/5kz/J5s6d227f1KlTs7/8y7/s1nn2BPmu1Zo1a7KioqLzNLueKyKyDRs2nHXMt771rezqq69ut2/evHlZZWVlN86s5+nMWv30pz/NIiJ76623zsucerLXXnsti4hs165dHY7pz69Z/19n1spr1m9ccskl2UMPPXTG36V6TvXZMyPvvvtu7Nu3L2bPnt22b8CAATF79uyora094zG1tbXtxkdEVFZWdji+r+jKWkVEnDp1KkaPHh1lZWVnLe3+rr8+rz6KCRMmRGlpaXzhC1+Ip59+OvV0kmhoaIiIiGHDhnU4xnPrPZ1ZqwivWS0tLbFu3bpoamrq8F+ypHpO9dkYeeONN6KlpaXtm2HfV1xc3OH7z3V1dXmN7yu6slbjx4+Phx9+ODZt2hSPPvpotLa2xvTp0+Pll18+H1PuVTp6XjU2Nsavf/3rRLPqmUpLS2P16tXx+OOPx+OPPx5lZWUxc+bM2L9/f+qpnVetra2xZMmSuP7668/6bdX99TXr/+vsWvXn16xDhw7FRRddFAUFBfG1r30tNmzYEFddddUZx6Z6TuX9dfAQ8d4/QPz/ZT19+vS48sor43vf+1585zvfSTgzerPx48fH+PHj236ePn16PP/883H//ffH97///YQzO78WLVoUhw8fjqeeeir1VHq8zq5Vf37NGj9+fBw8eDAaGhrihz/8YSxYsCB27drVYZCk0GfPjFx66aUxcODAqK+vb7e/vr4+SkpKznhMSUlJXuP7iq6s1W+74IILYuLEiXH06NHumGKv1tHzqrCwMC688MJEs+o9pkyZ0q+eV4sXL47NmzfHT3/60xg1atRZx/bX16z35bNWv60/vWYNHjw4xo0bF5MmTYrq6uooLy+Pf/7nfz7j2FTPqT4bI4MHD45JkybF9u3b2/a1trbG9u3bO3yvrKKiot34iIht27Z1OL6v6Mpa/baWlpY4dOhQlJaWdtc0e63++rw6Vw4ePNgvnldZlsXixYtjw4YNsWPHjhg7duyHHtNfn1tdWavf1p9fs1pbW6O5ufmMv0v2nOrWy2MTW7duXVZQUJDV1NRk//Vf/5X9xV/8RXbxxRdndXV1WZZl2Ze+9KVs2bJlbeOffvrpbNCgQdm9996bPffcc9kdd9yRXXDBBdmhQ4dSPYTzJt+1uvPOO7Mnn3wye/7557N9+/Zlf/qnf5oNGTIke/bZZ1M9hPPmnXfeyQ4cOJAdOHAgi4jsvvvuyw4cOJAdO3Ysy7IsW7ZsWfalL32pbfwLL7yQfexjH8v+5m/+JnvuueeyVatWZQMHDsyeeOKJVA/hvMl3re6///5s48aN2S9/+cvs0KFD2W233ZYNGDAg+8lPfpLqIZw3X//617OioqJs586d2cmTJ9u2X/3qV21jvGa9pytr1V9fs5YtW5bt2rUre/HFF7Of//zn2bJly7JcLpdt3bo1y7Ke85zq0zGSZVn23e9+N7vsssuywYMHZ1OmTMmeeeaZtt/NmDEjW7BgQbvxP/jBD7IrrrgiGzx4cHb11VdnW7ZsOc8zTieftVqyZEnb2OLi4uyLX/xitn///gSzPv/e//jpb2/vr8+CBQuyGTNmfOCYCRMmZIMHD84++clPZmvWrDnv804h37W6++67s8svvzwbMmRINmzYsGzmzJnZjh070kz+PDvTOkVEu+eK16z3dGWt+utr1le+8pVs9OjR2eDBg7NPfOIT2e///u+3hUiW9ZznVC7Lsqx7z70AAHSsz14zAgD0DmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgqf8DgyUK1RABovAAAAAASUVORK5CYII=\n", "text/plain": [ - "<Figure size 432x288 with 1 Axes>" + "<Figure size 640x480 with 1 Axes>" ] }, - "metadata": { - "needs_background": "light" - }, + "metadata": {}, "output_type": "display_data" } ], @@ -206,7 +205,7 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 10, "id": "8e051b62-1e63-467d-9f43-6a48813a0e2b", "metadata": {}, "outputs": [], @@ -218,7 +217,7 @@ }, { "cell_type": "code", - "execution_count": 20, + "execution_count": 11, "id": "4ea29864-913c-4dbe-8751-4be371b6d396", "metadata": {}, "outputs": [], @@ -228,7 +227,7 @@ }, { "cell_type": "code", - "execution_count": 21, + "execution_count": 12, "id": "4b09a818-f03a-45ef-996e-bcd431e9a8f4", "metadata": {}, "outputs": [], @@ -239,7 +238,28 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 16, + "id": "5e951704", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "<coffea.jetmet_tools.CorrectedJetsFactory.CorrectedJetsFactory at 0x7f16a65fc190>" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "jmestuff['fatjet_factory']['2018']" + ] + }, + { + "cell_type": "code", + "execution_count": 17, "id": "6b9387c1-28d0-4e28-bf35-bebce40cfea7", "metadata": {}, "outputs": [], @@ -254,7 +274,7 @@ }, { "cell_type": "code", - "execution_count": 23, + "execution_count": 18, "id": "10b431d3-6bab-4a0e-a474-2f5257e8623e", "metadata": {}, "outputs": [], @@ -265,7 +285,7 @@ }, { "cell_type": "code", - "execution_count": 27, + "execution_count": 19, "id": "0f0255f9-238c-461c-a003-1113e353996c", "metadata": {}, "outputs": [ @@ -273,41 +293,87 @@ "name": "stdout", "output_type": "stream", "text": [ - "CPU times: user 3.87 s, sys: 616 ms, total: 4.49 s\n", - "Wall time: 4.48 s\n" + "CPU times: user 2.3 s, sys: 165 ms, total: 2.46 s\n", + "Wall time: 7.43 s\n" ] } ], "source": [ "%%time\n", - "j = jmestuff[\"jet_factory\"][\"2018\"].build(add_jec_variables(events.Jet,events.fixedGridRhoFastjetAll), jec_cache)" + "#j = jmestuff[\"jet_factory\"][\"2018\"].build(add_jec_variables(events.Jet,events.fixedGridRhoFastjetAll), jec_cache)\n", + "j = jmestuff[\"fatjet_factory\"][\"2018\"].build(add_jec_variables(events.FatJet,events.fixedGridRhoFastjetAll), jec_cache)" ] }, { "cell_type": "code", - "execution_count": 28, + "execution_count": 20, "id": "0f86f278-ab3e-4e3f-9c18-f2820cb4f840", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "['JES_FlavorQCD',\n", + "['JES_AbsoluteStat',\n", + " 'JES_AbsoluteScale',\n", + " 'JES_AbsoluteSample',\n", + " 'JES_AbsoluteFlavMap',\n", + " 'JES_AbsoluteMPFBias',\n", + " 'JES_Fragmentation',\n", + " 'JES_SinglePionECAL',\n", + " 'JES_SinglePionHCAL',\n", + " 'JES_FlavorQCD',\n", + " 'JES_TimePtEta',\n", + " 'JES_RelativeJEREC1',\n", + " 'JES_RelativeJEREC2',\n", + " 'JES_RelativeJERHF',\n", + " 'JES_RelativePtBB',\n", + " 'JES_RelativePtEC1',\n", + " 'JES_RelativePtEC2',\n", + " 'JES_RelativePtHF',\n", " 'JES_RelativeBal',\n", - " 'JES_HF',\n", - " 'JES_BBEC1',\n", - " 'JES_EC2',\n", - " 'JES_Absolute',\n", - " 'JES_Absolute_2018',\n", - " 'JES_HF_2018',\n", - " 'JES_EC2_2018',\n", - " 'JES_RelativeSample_2018',\n", - " 'JES_BBEC1_2018',\n", + " 'JES_RelativeSample',\n", + " 'JES_RelativeFSR',\n", + " 'JES_RelativeStatFSR',\n", + " 'JES_RelativeStatEC',\n", + " 'JES_RelativeStatHF',\n", + " 'JES_PileUpDataMC',\n", + " 'JES_PileUpPtRef',\n", + " 'JES_PileUpPtBB',\n", + " 'JES_PileUpPtEC1',\n", + " 'JES_PileUpPtEC2',\n", + " 'JES_PileUpPtHF',\n", + " 'JES_PileUpMuZero',\n", + " 'JES_PileUpEnvelope',\n", + " 'JES_SubTotalPileUp',\n", + " 'JES_SubTotalRelative',\n", + " 'JES_SubTotalPt',\n", + " 'JES_SubTotalScale',\n", + " 'JES_SubTotalAbsolute',\n", + " 'JES_SubTotalMC',\n", " 'JES_Total',\n", + " 'JES_TotalNoFlavor',\n", + " 'JES_TotalNoTime',\n", + " 'JES_TotalNoFlavorNoTime',\n", + " 'JES_FlavorZJet',\n", + " 'JES_FlavorPhotonJet',\n", + " 'JES_FlavorPureGluon',\n", + " 'JES_FlavorPureQuark',\n", + " 'JES_FlavorPureCharm',\n", + " 'JES_FlavorPureBottom',\n", + " 'JES_TimeRunA',\n", + " 'JES_TimeRunB',\n", + " 'JES_TimeRunC',\n", + " 'JES_TimeRunD',\n", + " 'JES_CorrelationGroupMPFInSitu',\n", + " 'JES_CorrelationGroupIntercalibration',\n", + " 'jet_energy_uncertainty_CorrelationGroupbJES',\n", + " 'JES_CorrelationGroupbJES',\n", + " 'JES_CorrelationGroupFlavor',\n", + " 'JES_CorrelationGroupUncorrelated',\n", " 'JES_jes']" ] }, - "execution_count": 28, + "execution_count": 20, "metadata": {}, "output_type": "execute_result" } @@ -662,7 +728,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.15" + "version": "3.8.15" } }, "nbformat": 4, -- GitLab From a6356bf358e2c472011b5723ae068a82481f8501 Mon Sep 17 00:00:00 2001 From: mmarchegiani <matmarcheg@gmail.com> Date: Fri, 10 Feb 2023 16:32:17 +0100 Subject: [PATCH 02/11] Temporary solution for plotting ggH --- pocket_coffea/utils/plot_utils.py | 78 ++++++++++++++++++++++++++++++- scripts/plot/make_plots.py | 4 +- 2 files changed, 79 insertions(+), 3 deletions(-) diff --git a/pocket_coffea/utils/plot_utils.py b/pocket_coffea/utils/plot_utils.py index ae7dc0fa..df6004f5 100644 --- a/pocket_coffea/utils/plot_utils.py +++ b/pocket_coffea/utils/plot_utils.py @@ -277,6 +277,9 @@ def plot_data_mc_hist1D( samples = h.keys() samples_data = list(filter(lambda d: 'DATA' in d, samples)) samples_mc = list(filter(lambda d: 'DATA' not in d, samples)) + + print("********************") + print(samples_mc) h_mc = h[samples_mc[0]] @@ -292,6 +295,8 @@ def plot_data_mc_hist1D( for year in years: for cat in categories: + if 'noreweight' in cat: + continue if only_cat: if isinstance(only_cat, list): if not cat in only_cat: @@ -370,7 +375,7 @@ def plot_data_mc_hist1D( [ h[d][slicing_mc] for d in samples_mc - if d.endswith(f'_{f}') + if (d.endswith(f'_{f}') & ('GluGluH' not in d)) ] ) ) @@ -382,7 +387,7 @@ def plot_data_mc_hist1D( [ h[d][slicing_mc_nominal] for d in samples_mc - if d.endswith(f'_{f}') + if (d.endswith(f'_{f}') & ('GluGluH' not in d)) ] ) ) @@ -392,6 +397,57 @@ def plot_data_mc_hist1D( nevents = { f: round(sum(dict_mc_nominal[f].values()), 1) for f in flavors } + print("************") + print(histname) + print(samples_mc) + dict_gghbb = { + f: stack_sum( + hist.Stack.from_iter( + [ + h[d][slicing_mc] + for d in samples_mc + if (d.endswith(f'_{f}') & ('GluGluHToBB' in d)) + ] + ) + ) + for f in flavors + } + dict_gghbb_nominal = { + f: stack_sum( + hist.Stack.from_iter( + [ + h[d][slicing_mc_nominal] + for d in samples_mc + if (d.endswith(f'_{f}') & ('GluGluHToBB' in d)) + ] + ) + ) + for f in flavors + } + dict_gghcc = { + f: stack_sum( + hist.Stack.from_iter( + [ + h[d][slicing_mc] + for d in samples_mc + if (d.endswith(f'_{f}') & ('GluGluHToCC' in d)) + ] + ) + ) + for f in flavors + } + dict_gghcc_nominal = { + f: stack_sum( + hist.Stack.from_iter( + [ + h[d][slicing_mc_nominal] + for d in samples_mc + if (d.endswith(f'_{f}') & ('GluGluHToCC' in d)) + ] + ) + ) + for f in flavors + } else: dict_mc = {d: h[d][slicing_mc] for d in samples_mc} dict_mc_nominal = {d: h[d][slicing_mc_nominal] for d in samples_mc} @@ -406,6 +462,15 @@ def plot_data_mc_hist1D( dict_mc_nominal[sample] = histo_reweighted stack_mc = hist.Stack.from_dict(dict_mc) stack_mc_nominal = hist.Stack.from_dict(dict_mc_nominal) + if flavorsplit == '5f': + stack_gghbb = hist.Stack.from_dict(dict_gghbb) + stack_gghbb_nominal = hist.Stack.from_dict(dict_gghbb_nominal) + stack_gghcc = hist.Stack.from_dict(dict_gghcc) + stack_gghcc_nominal = hist.Stack.from_dict(dict_gghcc_nominal) + hist_gghbb = stack_sum(stack_gghbb) + hist_gghbb_nominal = stack_sum(stack_gghbb_nominal) + hist_gghcc = stack_sum(stack_gghcc) + hist_gghcc_nominal = stack_sum(stack_gghcc_nominal) if not is_mc_only: # Sum over eras if era axis exists in data histogram @@ -463,6 +528,12 @@ def plot_data_mc_hist1D( x, h_data.values(), yerr=np.sqrt(h_data.values()), **opts_data ) # stack_data.plot(stack=True, color='black', ax=ax) + if flavorsplit == '5f': + sf = 100 + hist_gghbb_nominal = sf * hist_gghbb_nominal + hist_gghcc_nominal = sf * hist_gghcc_nominal + hist_gghbb_nominal.plot(stack=False, histtype='step', ax=ax, color='red', label=f'ggHbb x {sf}') + hist_gghcc_nominal.plot(stack=False, histtype='step', ax=ax, color='green', label=f'ggHcc x {sf}') if rebinning: syst_err_up, syst_err_down = get_systematic_uncertainty( stack_mc, @@ -515,6 +586,7 @@ def plot_data_mc_hist1D( rax.tick_params(axis='y', labelsize=fontsize) if log: + ax.set_yscale("log") ax.set_xlim(0, 1) exp = math.floor( math.log(max(stack_sum(stack_mc_nominal).values()), 10) @@ -529,6 +601,7 @@ def plot_data_mc_hist1D( if histname in config.plot_options.keys(): if 'xlim' in config.plot_options[histname].keys(): ax.set_xlim(*config.plot_options[histname]['xlim']) + """ if 'scale' in config.plot_options[histname].keys(): ax.set_yscale(config.plot_options[histname]['scale']) if config.plot_options[histname]['scale'] == 'log': @@ -539,6 +612,7 @@ def plot_data_mc_hist1D( ) ax.set_ylim((0.01, 10 ** (exp + 2))) # ax.legend(handles, labels, loc="upper right", fontsize=fontsize, ncols=2) + """ if 'ylim' in config.plot_options[histname].keys(): if isinstance(config.plot_options[histname]['ylim'], tuple): ax.set_ylim(*config.plot_options[histname]['ylim']) diff --git a/scripts/plot/make_plots.py b/scripts/plot/make_plots.py index e998b553..76d1edef 100644 --- a/scripts/plot/make_plots.py +++ b/scripts/plot/make_plots.py @@ -112,7 +112,9 @@ if not os.path.exists(config.plots): def make_plots(entrystart, entrystop): _accumulator = slice_accumulator(accumulator, entrystart, entrystop) for (histname, h) in _accumulator['variables'].items(): - plot_data_mc_hist1D(h, histname, config, flavorsplit=None, mcstat=True, stat_only=False) + #if ('FatJetGood_DDX' in histname) | ('FatJetGood_particleNetMD_Xqq' in histname) | ('FatJetGood_particleNet_H4qvsQCD' in histname): + # continue + plot_data_mc_hist1D(h, histname, config, flavorsplit='5f', mcstat=True, stat_only=False, log=True) # Filter dictionary of histograms with `args.only` accumulator['variables'] = { k : v for k,v in accumulator['variables'].items() if args.only in k } -- GitLab From b63fe558a07113ab6c698ea379f11b4e0706b546 Mon Sep 17 00:00:00 2001 From: mmarchegiani <matmarcheg@gmail.com> Date: Thu, 16 Mar 2023 11:04:59 +0100 Subject: [PATCH 03/11] Nicer log plots --- pocket_coffea/utils/plot_utils.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pocket_coffea/utils/plot_utils.py b/pocket_coffea/utils/plot_utils.py index 00c0d4ab..09d688ae 100644 --- a/pocket_coffea/utils/plot_utils.py +++ b/pocket_coffea/utils/plot_utils.py @@ -728,7 +728,7 @@ def plot_data_mc_hist1D( exp = math.floor( math.log(max(stack_sum(stack_mc_nominal).values()), 10) ) - ax.set_ylim((0.01, 10 ** (exp + 2))) + ax.set_ylim((0.01, 10 ** (exp + 3))) # if flavorsplit == '5f': # ax.legend(handles, labels, loc="upper right", fontsize=fontsize, ncols=2) # else: @@ -747,7 +747,7 @@ def plot_data_mc_hist1D( handles_new.append(handles[i]) labels = labels_new handles = handles_new - ax.legend(handles, labels, fontsize=fontsize, ncols=2) + ax.legend(handles, labels, fontsize=fontsize, ncols=2, loc="upper right") if ("variables" in config.plot_options) & (histname in config.plot_options["variables"].keys()): cfg_plot = config.plot_options["variables"][histname] if 'xlim' in cfg_plot.keys(): @@ -760,7 +760,7 @@ def plot_data_mc_hist1D( max(stack_sum(stack_mc_nominal).values()), 10 ) ) - ax.set_ylim((0.01, 10 ** (exp + 2))) + ax.set_ylim((0.01, 10 ** (exp + 3))) # ax.legend(handles, labels, loc="upper right", fontsize=fontsize, ncols=2) if 'ylim' in cfg_plot.keys(): if isinstance(cfg_plot['ylim'], tuple): -- GitLab From 05a3c61438e2cc6a828fe0026456ae3d20005ca6 Mon Sep 17 00:00:00 2001 From: mmarchegiani <matmarcheg@gmail.com> Date: Mon, 20 Mar 2023 23:31:20 +0100 Subject: [PATCH 04/11] Update genWeightsProcessor for skim --- pocket_coffea/workflows/genweights.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/pocket_coffea/workflows/genweights.py b/pocket_coffea/workflows/genweights.py index 1a404377..00fe5d08 100644 --- a/pocket_coffea/workflows/genweights.py +++ b/pocket_coffea/workflows/genweights.py @@ -48,4 +48,13 @@ class genWeightsProcessor(BaseProcessorABC): if self._hasSubsamples: raise Exception("This processor cannot compute the sum of genweights of subsamples.") + self.process_extra_before_skim() + self.skim_events() + if not self.has_events: + return self.output + + if self.cfg.save_skimmed_files: + self.export_skimmed_chunk() + return self.output + return self.output -- GitLab From 35de55d76ed059c8e85c13e07699bc554acfe966 Mon Sep 17 00:00:00 2001 From: mmarchegiani <matmarcheg@gmail.com> Date: Tue, 21 Mar 2023 01:13:54 +0100 Subject: [PATCH 05/11] Updates for skim --- pocket_coffea/workflows/genweights.py | 48 ++++++++++++++++++++++++++- 1 file changed, 47 insertions(+), 1 deletion(-) diff --git a/pocket_coffea/workflows/genweights.py b/pocket_coffea/workflows/genweights.py index 00fe5d08..dc0ffab7 100644 --- a/pocket_coffea/workflows/genweights.py +++ b/pocket_coffea/workflows/genweights.py @@ -1,9 +1,16 @@ +import numpy as np import awkward as ak import copy +from coffea.lumi_tools import LumiMask +from coffea.analysis_tools import PackedSelection + from .base import BaseProcessorABC from ..utils.configurator import Configurator +from ..parameters.event_flags import event_flags, event_flags_data +from ..parameters.lumi import goldenJSON +from ..parameters.btag import btag class genWeightsProcessor(BaseProcessorABC): def __init__(self, cfg: Configurator) -> None: @@ -11,7 +18,8 @@ class genWeightsProcessor(BaseProcessorABC): self.output_format = { "sum_genweights": {}, "cutflow": { - "initial": {s: 0 for s in self.cfg.datasets} + "initial": {s: 0 for s in self.cfg.datasets}, + "skim": {s: 0 for s in self.cfg.datasets} } } @@ -19,9 +27,14 @@ class genWeightsProcessor(BaseProcessorABC): self._dataset = self.events.metadata["dataset"] self._sample = self.events.metadata["sample"] self._year = self.events.metadata["year"] + self._btag = btag[self._year] self._isMC = self.events.metadata["isMC"] == "True" if self._isMC: + self._era = "MC" self._xsec = self.events.metadata["xsec"] + else: + self._era = self.events.metadata["era"] + self._goldenJSON = goldenJSON[self._year] # Check if the user specified any subsamples without performing any operation if self._sample in self._subsamplesCfg: @@ -29,6 +42,39 @@ class genWeightsProcessor(BaseProcessorABC): else: self._hasSubsamples = False + def skim_events(self): + self._skim_masks = PackedSelection() + mask_flags = np.ones(self.nEvents_initial, dtype=np.bool) + flags = event_flags[self._year] + if not self._isMC: + flags += event_flags_data[self._year] + for flag in flags: + mask_flags = getattr(self.events.Flag, flag) + self._skim_masks.add("event_flags", mask_flags) + + self._skim_masks.add("PVgood", self.events.PV.npvsGood > 0) + + if not self._isMC: + mask_lumi = LumiMask(self._goldenJSON)( + self.events.run, self.events.luminosityBlock + ) + self._skim_masks.add("lumi_golden", mask_lumi) + + for skim_func in self._skim: + mask = skim_func.get_mask( + self.events, + year=self._year, + sample=self._sample, + btag=self._btag, + isMC=self._isMC, + ) + self._skim_masks.add(skim_func.id, mask) + + self.events = self.events[self._skim_masks.all(*self._skim_masks.names)] + self.nEvents_after_skim = self.nevents + self.output['cutflow']['skim'][self._dataset] += self.nEvents_after_skim + self.has_events = self.nEvents_after_skim > 0 + def apply_object_preselection(self, variation): pass -- GitLab From 22cc6e329ec02c2e6c5b81b82f35be0211dc9901 Mon Sep 17 00:00:00 2001 From: mmarchegiani <matmarcheg@gmail.com> Date: Tue, 21 Mar 2023 13:28:35 +0100 Subject: [PATCH 06/11] Documenting hadd script --- scripts/README.md | 85 +++++++++++++++++++++++++++++++++++ scripts/hadd_skimmed_files.py | 9 ++-- 2 files changed, 91 insertions(+), 3 deletions(-) create mode 100644 scripts/README.md diff --git a/scripts/README.md b/scripts/README.md new file mode 100644 index 00000000..940cc5a1 --- /dev/null +++ b/scripts/README.md @@ -0,0 +1,85 @@ +# PocketCoffea scripts + +The current guide is a documentation of how to run the scripts of the PocketCoffea package. + +## Merge skimmed files with hadd + +In order to merge the skimmed n-tuples, the `hadd` command from ROOT needs to be used. +To source ROOT inside the `pocket_coffea` environment, run the following command: + +``` +source /cvmfs/sft.cern.ch/lcg/views/LCG_102/x86_64-centos7-gcc11-opt/setup.sh +``` + +To run the `hadd` script run the following command: + +``` +cd /path/to/PocketCoffea +python scripts/hadd_skimmed_files.py -fl /path/to/output/output_all.coffea -o /path/to/outputdir +``` + +## Build dataset + +The first step is to build the json files containing the list of files of the UL datasets to process. In our case, we need to include the `TTToSemiLeptonic` and `TTTo2L2Nu` MC datasets and `SingleMuon` as data dataset. + +A step-by-step detailed guide on the creation of the dataset file can be found at this [link](https://pocketcoffea.readthedocs.io/en/latest/examples.html). + +### Append sum of genweights to datasets definitions + +In order to update the `datasets_definitions.json` and include the `sum_genweights` metadata in the datasets definitions, a dedicated script needs to be run with the following command: + +``` +cd /path/to/PocketCoffea +scripts/dataset/append_genweights.py --cfg dataset/dataset_definitions.json -i output/genweights/genweights_2016_PreVFP/output_all.coffea --overwrite +scripts/dataset/append_genweights.py --cfg dataset/dataset_definitions.json -i output/genweights/genweights_2016_PostVFP/output_all.coffea --overwrite +``` + +N.B.: the argument `-i output/genweights/genweights_201*/output_all.coffea` should be the path to the Coffea output file where the genweights are stored. The option `--overwrite` automatically updates the config file passed as `--cfg` argument with the genweights, in this case `dataset/dataset_definitions.json`. + +### Build dataset with sum of genweights + +Now the steps of the [Build dataset](#build-dataset) section need to be repeated to generate the json datasets with the additional metadata `sum_genweights`. + +## Run the analysis + +In order to run the analysis workflow and produce the output histograms, run the following command: + +``` +cd /path/to/PocketCoffea +runner.py --cfg /path/to/config/config.py --full +``` + +N.B.: the argument `--full` will process all the datasets together at once and save the output in a single output file, `output_all.coffea`. Otherwise, the datasets are processed separately and an output file is saved for each dataset. + +## Accumulate output files + +If the output of the [previous step](#run-the-analysis) has been produced without the argument `--full`, the output files need to be merged in a single output file `output_all.coffea`. If the output has been produced with the argument `--full` and the output file `output_all.coffea` is already existing, skip this step and continue with the [next one](#produce-datamc-plots). + +Once the Coffea output files are produced, one needs to merge the files into a single file by using the script `accumulate_files.py` by running this command: + +``` +cd /path/to/output +python accumulate_files.py -i output1.coffea output2.coffea output3.coffea -o output_all.coffea +``` + +## Produce data/MC plots + +To produce plots from a Coffea output file execute the plotting script: + +``` +cd /path/to/PocketCoffea +make_plots.py --cfg /path/to/config/config.py -i /path/to/output/output_all.coffea +``` + +## Run trigger SF script + +To run the script that computes the single electron trigger SF and produces the validation plots, run: + +``` +cd /path/to/PocketCoffea +python scripts/plot/trigger_efficiency.py --cfg config/semileptonic_triggerSF/semileptonic_triggerSF_2017.py --save_plots +``` + +If the additional argument `--save_plots` is not specified, only the scale factor maps are saved without saving the plots. + +The output plots are saved in `/path/to/output/plots/trigger_efficiency` and `/path/to/output/plots/trigger_scalefactor`, while the 1D and 2D scale factor maps are saved in the folder specified in `workflow_options['output_triggerSF']` in the config file. diff --git a/scripts/hadd_skimmed_files.py b/scripts/hadd_skimmed_files.py index 2280fcca..6d870ba2 100644 --- a/scripts/hadd_skimmed_files.py +++ b/scripts/hadd_skimmed_files.py @@ -18,7 +18,7 @@ parser.add_argument( ) parser.add_argument("-o", "--outputdir", required=False, type=str, help="Output folder") parser.add_argument( - "--only-datasets", nargs="+", type=str, help="Restring list of datasets" + "--only-datasets", nargs="+", type=str, help="Restricting list of datasets" ) parser.add_argument("-f", "--files", type=int, help="Limit number of files") parser.add_argument("-e", "--events", type=int, help="Limit number of files") @@ -43,6 +43,9 @@ groups_metadata = {} for dataset in df["skimmed_files"].keys(): if args.only_datasets and dataset not in args.only_datasets: continue + outputdir_dataset = os.path.join(args.outputdir, dataset) + if not os.path.exists(outputdir_dataset): + os.makedirs(outputdir_dataset) groups_metadata[dataset] = defaultdict(dict) nevents_tot = 0 nfiles = 0 @@ -54,7 +57,7 @@ for dataset in df["skimmed_files"].keys(): if (args.files and (nfiles + 1) > args.files) or ( args.events and (nevents_tot + nevents) > args.events ): - outfile = f"{args.outputdir}/{dataset}/{dataset}_{ngroup}.root" + outfile = os.path.join(outputdir_dataset, f"{dataset}_{ngroup}.root") workload.append((outfile, group[:])) groups_metadata[dataset]["files"][outfile] = group[:] group.clear() @@ -68,7 +71,7 @@ for dataset in df["skimmed_files"].keys(): # add last group if len(group): - outfile = f"{args.outputdir}/{dataset}/{dataset}_{ngroup}.root" + outfile = os.path.join(outputdir_dataset, f"{dataset}_{ngroup}.root") workload.append((outfile, group[:])) groups_metadata[dataset]["files"][outfile] = group[:] -- GitLab From 2cdc97f3411628950176da20d345a999ab807b1f Mon Sep 17 00:00:00 2001 From: mmarchegiani <matmarcheg@gmail.com> Date: Wed, 22 Mar 2023 03:06:57 +0100 Subject: [PATCH 07/11] Update object preselection for mutag --- pocket_coffea/parameters/object_preselection.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/pocket_coffea/parameters/object_preselection.py b/pocket_coffea/parameters/object_preselection.py index e93c0453..291f816d 100644 --- a/pocket_coffea/parameters/object_preselection.py +++ b/pocket_coffea/parameters/object_preselection.py @@ -89,9 +89,6 @@ object_preselection = { "eta": 2.4, "msd": 40, "jetId": 2, - "nsubjet": 2, - "nmusj": 1, - "dimuon_pt_ratio": 0.6, }, }, } -- GitLab From ff7cb05bfd6a8472aa4849c1965e29b32eaa695e Mon Sep 17 00:00:00 2001 From: mmarchegiani <matmarcheg@gmail.com> Date: Fri, 24 Mar 2023 02:02:36 +0100 Subject: [PATCH 08/11] Define fatjet selection --- pocket_coffea/lib/jets.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/pocket_coffea/lib/jets.py b/pocket_coffea/lib/jets.py index 79051fde..07e37cee 100644 --- a/pocket_coffea/lib/jets.py +++ b/pocket_coffea/lib/jets.py @@ -226,7 +226,7 @@ def jet_selection(events, Jet, finalstate): # Only jets that are more distant than dr to ALL leptons are tagged as good jets leptons = events["LeptonGood"] # Mask for jets not passing the preselection - presel_mask = ( + mask_presel = ( (jets.pt > cuts["pt"]) & (np.abs(jets.eta) < cuts["eta"]) & (jets.jetId >= cuts["jetId"]) @@ -234,18 +234,20 @@ def jet_selection(events, Jet, finalstate): # Lepton cleaning if "dr" in cuts.keys(): dR_jets_lep = jets.metric_table(leptons) - lepton_cleaning_mask = ak.prod(dR_jets_lep > cuts["dr"], axis=2) == 1 + mask_lepton_cleaning = ak.prod(dR_jets_lep > cuts["dr"], axis=2) == 1 if Jet == "Jet": - jetpuid_mask = (jets.puId >= cuts["puId"]["value"]) | ( + mask_jetpuid = (jets.puId >= cuts["puId"]["value"]) | ( jets.pt >= cuts["puId"]["maxpt"] ) - good_jets_mask = presel_mask & lepton_cleaning_mask & jetpuid_mask + mask_good_jets = mask_presel & mask_lepton_cleaning & mask_jetpuid elif Jet == "FatJet": - raise NotImplementedError() + # Apply the msd and preselection cuts + mask_msd = (events.FatJet.msoftdrop > cuts["msd"]) + mask_good_jets = mask_presel & mask_msd - return jets[good_jets_mask], good_jets_mask + return jets[mask_good_jets], mask_good_jets def btagging(Jet, btag): -- GitLab From 6fb0e4bd7757fc5ffd47ed7a94516c4f75ecb5ca Mon Sep 17 00:00:00 2001 From: mmarchegiani <matmarcheg@gmail.com> Date: Fri, 24 Mar 2023 02:43:50 +0100 Subject: [PATCH 09/11] Define FatJetGood --- pocket_coffea/utils/configurator.py | 12 ++++++++++++ pocket_coffea/workflows/tthbb_base_processor.py | 4 ++++ 2 files changed, 16 insertions(+) diff --git a/pocket_coffea/utils/configurator.py b/pocket_coffea/utils/configurator.py index a3c626a1..8b5bac46 100644 --- a/pocket_coffea/utils/configurator.py +++ b/pocket_coffea/utils/configurator.py @@ -209,6 +209,18 @@ class Configurator: self.subsamples_list += sam self.total_samples_list = list(set(self.samples + self.subsamples_list)) + for key, subscfg in self.subsamples_cuts.items(): + if isinstance(subscfg, dict): + # Convert it to StandardSelection + self.subsamples[key] = StandardSelection(subscfg) + elif isinstance(subscfg, StandardSelection): + self.subsamples[key] = subscfg + elif isinstance(subscfg, CartesianSelection): + self.subsamples[key] = subscfg + # Unique set of cuts + logging.info("Subsamples:") + logging.info(self.subsamples) + def filter_dataset(self, nfiles): filtered_dataset = {} for sample, ds in self.fileset.items(): diff --git a/pocket_coffea/workflows/tthbb_base_processor.py b/pocket_coffea/workflows/tthbb_base_processor.py index 8602188b..5c035380 100644 --- a/pocket_coffea/workflows/tthbb_base_processor.py +++ b/pocket_coffea/workflows/tthbb_base_processor.py @@ -66,6 +66,10 @@ class ttHbbBaseProcessor(BaseProcessorABC): self.events.ElectronGood, self.events.MuonGood ) + self.events["FatJetGood"], self.jetGoodMask = jet_selection( + self.events, "FatJet", self.cfg.finalstate + ) + def count_objects(self, variation): self.events["nMuonGood"] = ak.num(self.events.MuonGood) self.events["nElectronGood"] = ak.num(self.events.ElectronGood) -- GitLab From d90c609c837ff6465cf19906ac133b25e27a4055 Mon Sep 17 00:00:00 2001 From: mmarchegiani <matmarcheg@gmail.com> Date: Fri, 24 Mar 2023 02:52:28 +0100 Subject: [PATCH 10/11] Treat subsamples as StandardSelection or CartesianSelection --- config/test_subsamples.py | 107 +++++++++++++++++----------- pocket_coffea/lib/hist_manager.py | 4 +- pocket_coffea/utils/configurator.py | 39 +++++----- pocket_coffea/workflows/base.py | 66 +++++------------ 4 files changed, 105 insertions(+), 111 deletions(-) diff --git a/config/test_subsamples.py b/config/test_subsamples.py index a75afa41..cdc64b6f 100644 --- a/config/test_subsamples.py +++ b/config/test_subsamples.py @@ -1,12 +1,65 @@ from pocket_coffea.parameters.cuts.preselection_cuts import * from pocket_coffea.workflows.tthbb_base_processor import ttHbbBaseProcessor -from pocket_coffea.lib.cut_functions import get_nObj_min, get_nObj_eq, get_nBtag, get_HLTsel, get_nObj_less +from pocket_coffea.lib.cut_functions import get_nObj_min, get_nObj_eq, get_HLTsel, get_nObj_less from pocket_coffea.parameters.histograms import * from pocket_coffea.lib.columns_manager import ColOut from pocket_coffea.parameters.btag import btag_variations -from pocket_coffea.lib.weights_manager import WeightCustom -import numpy as np +def ptmsd(events, params, **kwargs): + # Mask to select events with a fatjet with minimum softdrop mass and maximum tau21 + #return (events.FatJetGood[:,0].pt > params["pt"]) & (events.FatJetGood[:,0].msoftdrop > params["msd"]) + mask = (events.FatJetGood.pt > params["pt"]) & (events.FatJetGood.msoftdrop > params["msd"]) + + assert not ak.any(ak.is_none(mask, axis=1)), f"None in ptmsd\n{events.FatJetGood.pt[ak.is_none(mask, axis=1)]}" + + return ak.where(~ak.is_none(mask, axis=1), mask, False) + +def get_ptmsd(pt, msd, name=None): + if name == None: + name = f"pt{pt}msd{msd}" + return Cut( + name=name, + params= {"pt" : pt, "msd" : msd}, + function=ptmsd, + collection="FatJetGood" + ) + +def flavor_mask(events, params, **kwargs): + mask = { + "l" : events.FatJetGood.hadronFlavour < 4, + "c" : events.FatJetGood.hadronFlavour == 4, + "b" : events.FatJetGood.hadronFlavour == 5, + "cc" : abs(events.FatJetGood.hadronFlavour == 4) & (events.FatJetGood.nBHadrons == 0) & (events.FatJetGood.nCHadrons >= 2), + "bb" : abs(events.FatJetGood.hadronFlavour == 5) & (events.FatJetGood.nBHadrons >= 2) + } + + if params["flavor"] in ["bb", "cc"]: + return mask[params["flavor"]] + elif params["flavor"] == "b": + return mask[params["flavor"]] & ~mask["bb"] + elif params["flavor"] == "c": + return mask[params["flavor"]] & ~mask["cc"] + elif params["flavor"] == "l": + return mask[params["flavor"]] & ~mask["bb"] & ~mask["cc"] & ~mask["b"] & ~mask["c"] + else: + raise NotImplementedError + +def get_flavor(flavor): + return Cut( + name=flavor, + params={"flavor": flavor}, + function=flavor_mask, + collection="FatJetGood" + ) + +samples = ["TTToSemiLeptonic", + "ttHTobb", + "DATA_SingleMuon" + ] + +subsamples = {} +for s in filter(lambda x: 'DATA' not in x, samples): + subsamples[s] = {f"{s}_{f}" : [get_flavor(f)] for f in ['l', 'c', 'b', 'cc', 'bb']} cfg = { "dataset" : { @@ -16,35 +69,17 @@ cfg = { "datasets/DATA_SingleMuon_local.json", "datasets/DATA_SingleEle.json"], "filter" : { - "samples": [ - # "TTToSemiLeptonic", - "TTbbSemiLeptonic", - "ttHTobb", - "DATA_SingleMu", - # "DATA_SingleEle" - ], + "samples": samples, "samples_exclude" : [], "year": ['2018'] }, - "subsamples": { - "TTbbSemiLeptonic": - { - "tt<3b": [get_nObj_less(3, coll="BJetGood")], - "tt+3b": [get_nObj_eq(3, coll="BJetGood")], - "tt+4b": [get_nObj_eq(4, coll="BJetGood")], - "tt>4b": [get_nObj_min(5, coll="BJetGood")], - }, - "ttHTobb":{ - "ttH<4b": [get_nObj_less(4, coll="BJetGood")], - "ttH+4b": [get_nObj_min(4, coll="BJetGood")], - } - } + "subsamples": subsamples }, # Input and output files "workflow" : ttHbbBaseProcessor, - "output" : "output/test_subsamples", + "output" : "output/test/test_subsamples", "worflow_options" : {}, "run_options" : { @@ -67,16 +102,15 @@ cfg = { }, # Cuts and plots settings - "finalstate" : "semileptonic", - "skim": [get_nObj_min(4, 15., "Jet"), + "finalstate" : "mutag", + "skim": [get_nObj_min(1, 400., "FatJet"), get_HLTsel("semileptonic")], - "preselections" : [semileptonic_presel], + "preselections" : [semileptonic_presel, get_nObj_min(1, 450., "FatJetGood")], "categories": { "baseline": [passthrough], + "pt450msd40" : [get_ptmsd(450., 40.)], }, - - "weights": { "common": { "inclusive": ["genWeight","lumi","XS", @@ -111,23 +145,16 @@ cfg = { "variables": { - + **fatjet_hists(coll="FatJetGood", fields=["pt", "eta", "phi",]), + **fatjet_hists(coll="FatJetGood", fields=["pt", "eta", "phi",], pos=0), **ele_hists(coll="ElectronGood", pos=0), **muon_hists(coll="MuonGood", pos=0), }, "columns": { "common": { - "inclusive": [ ColOut("JetGood", ["pt", "eta","phi"]), - ColOut("ElectronGood", ["pt","eta","phi"])] + "inclusive": [ ColOut("FatJetGood", ["pt", "eta", "phi"])] }, - "bysample":{ - "ttHTobb": { - "inclusive": [ColOut("MuonGood", ["pt", "eta", "phi"])] - }, - "ttH+4b":{ - "inclusive": [ColOut("BJetGood", ["pt", "eta", "phi"])] - } - } + "bysample" : { subs : { "inclusive" : [ColOut("FatJetGood", ["hadronFlavour", "nBHadrons", "nCHadrons"])] } for subs in subsamples } } } diff --git a/pocket_coffea/lib/hist_manager.py b/pocket_coffea/lib/hist_manager.py index 5e1dd11d..4f65ad18 100644 --- a/pocket_coffea/lib/hist_manager.py +++ b/pocket_coffea/lib/hist_manager.py @@ -258,7 +258,7 @@ class HistManager: weights_manager, categories, shape_variation="nominal", - subsamples_masks=None, # This is a dictionary with name:ak.Array(bool) + subsamples=None, # This is a dictionary with name:ak.Array(bool) custom_fields=None, ): ''' @@ -367,7 +367,7 @@ class HistManager: # Mask the events, the weights and then flatten and remove the None correctly for category, cat_mask in categories.get_masks(): # loop directly on subsamples - for subsample, subs_mask in subsamples_masks.items(): + for subsample, subs_mask in subsamples.get_masks(): # logging.info(f"\t\tcategory {category}, subsample {subsample}") mask = cat_mask & subs_mask # Skip empty categories and subsamples diff --git a/pocket_coffea/utils/configurator.py b/pocket_coffea/utils/configurator.py index 8b5bac46..454b0741 100644 --- a/pocket_coffea/utils/configurator.py +++ b/pocket_coffea/utils/configurator.py @@ -186,37 +186,38 @@ class Configurator: def load_subsamples(self): # subsamples configuration - # Dictionary with subsample_name:[list of Cut ids] - self.subsamples_cuts = { - key: subscfg - for key, subscfg in self.dataset.get("subsamples", {}).items() - if key in self.samples - } - # Map of subsamples names + subsamples_dict = self.dataset.get("subsamples", {}) self.subsamples = {} + self.subsamples_names = {} self.has_subsamples = {} + # Save list of subsamples for each sample for sample in self.samples: - if sample in self.subsamples_cuts: - self.subsamples[sample] = list(self.subsamples_cuts[sample].keys()) + if sample in subsamples_dict.keys(): + self.subsamples_names[sample] = list(subsamples_dict[sample].keys()) self.has_subsamples[sample] = True else: - self.subsamples[sample] = [] self.has_subsamples[sample] = False # Complete list of samples and subsamples self.subsamples_list = [] - for sam in self.subsamples.values(): + for sam in self.subsamples_names.values(): self.subsamples_list += sam self.total_samples_list = list(set(self.samples + self.subsamples_list)) - for key, subscfg in self.subsamples_cuts.items(): - if isinstance(subscfg, dict): - # Convert it to StandardSelection - self.subsamples[key] = StandardSelection(subscfg) - elif isinstance(subscfg, StandardSelection): - self.subsamples[key] = subscfg - elif isinstance(subscfg, CartesianSelection): - self.subsamples[key] = subscfg + # Now saving the subsamples + for sample in self.samples: + if sample in subsamples_dict: + subscfg = subsamples_dict[sample] + if isinstance(subscfg, dict): + # Convert it to StandardSelection + self.subsamples[sample] = StandardSelection(subscfg) + elif isinstance(subscfg, StandardSelection): + self.subsamples[sample] = subscfg + elif isinstance(subscfg, CartesianSelection): + self.subsamples[sample] = subscfg + else: + # if there is no configured subsample, the full sample becomes its subsample + self.subsamples[sample] = StandardSelection({sample: [passthrough]}) # Unique set of cuts logging.info("Subsamples:") logging.info(self.subsamples) diff --git a/pocket_coffea/workflows/base.py b/pocket_coffea/workflows/base.py index bbcbc9fd..3d6f4d36 100644 --- a/pocket_coffea/workflows/base.py +++ b/pocket_coffea/workflows/base.py @@ -55,7 +55,7 @@ class BaseProcessorABC(processor.ProcessorABC, ABC): self._categories = self.cfg.categories # Subsamples configurations: special cuts to split a sample in subsamples - self._subsamplesCfg = self.cfg.subsamples_cuts + self._subsamples = self.cfg.subsamples # Weights configuration self.weights_config_allsamples = self.cfg.weights_config @@ -116,15 +116,7 @@ class BaseProcessorABC(processor.ProcessorABC, ABC): self._era = self.events.metadata["era"] self._goldenJSON = goldenJSON[self._year] # Loading metadata for subsamples - if self._sample in self._subsamplesCfg: - self._subsamples = self._subsamplesCfg[self._sample] - self._subsamples_names = list(self._subsamples.keys()) - self._hasSubsamples = True - else: - # if there is no configured subsample, the full sample becomes its subsample - self._subsamples = {self._sample: [passthrough]} - self._subsamples_names = [self._sample] - self._hasSubsamples = False + self._hasSubsamples = self.cfg.has_subsamples[self._sample] def load_metadata_extra(self): ''' @@ -276,28 +268,12 @@ class BaseProcessorABC(processor.ProcessorABC, ABC): isMC=self._isMC, ) - # Defining the subsamples cut - # saving all the cuts in a single selector - self._subsamples_masks = PackedSelection() - self._subsamples_cuts_ids = [] - # saving the map of cut ids for each subsample - self._subsamples_map = defaultdict(list) - self._subsamples_masks = PackedSelection() - - for subs, cuts in self._subsamples.items(): - for cut in cuts: - if cut.id not in self._subsamples_cuts_ids: - self._subsamples_masks.add( - cut.id, - cut.get_mask( - self.events, - year=self._year, - sample=self._sample, - isMC=self._isMC, - ), - ) - self._subsamples_cuts_ids.append(cut.id) - self._subsamples_map[subs].append(cut.id) + self._subsamples[self._sample].prepare( + events=self.events, + year=self._year, + sample=self._sample, + isMC=self._isMC, + ) def define_common_variables_before_presel(self, variation): ''' @@ -377,11 +353,7 @@ class BaseProcessorABC(processor.ProcessorABC, ABC): # If subsamples are defined we also save their metadata if self._hasSubsamples: - for subs in self._subsamples_names: - # Get the mask - subsam_mask = self._subsamples_masks.all( - *self._subsamples_map[subs] - ) + for subs, subsam_mask in self._categories.get_masks(): mask_withsub = mask_on_events & subsam_mask self.output["cutflow"][category][subs] = ak.sum(mask_withsub) if self._isMC: @@ -410,7 +382,7 @@ class BaseProcessorABC(processor.ProcessorABC, ABC): self.hists_managers = HistManager( self.cfg.variables, self._sample, - self._subsamples_names, + self._subsamples[self._sample].keys(), self._categories, variations_config=self.cfg.variations_config[self._sample], custom_axes=self.custom_axes, @@ -434,21 +406,17 @@ class BaseProcessorABC(processor.ProcessorABC, ABC): throught the HistManager. ''' # Filling the autofill=True histogram automatically - subs_masks = {} - for subs in self._subsamples_names: - # Get the mask - subs_masks[subs] = self._subsamples_masks.all(*self._subsamples_map[subs]) # Calling hist manager with the subsample masks self.hists_managers.fill_histograms( self.events, self.weights_manager, self._categories, - subsamples_masks=subs_masks, + subsamples=self._subsamples[self._sample], shape_variation=variation, custom_fields=self.custom_histogram_fields, ) # Saving the output - for subs in self._subsamples_names: + for subs in self._subsamples[self._sample].keys(): # Saving the output for each subsample for var, H in self.hists_managers.get_histograms(subs).items(): self.output["variables"][var][subs] = H @@ -466,7 +434,7 @@ class BaseProcessorABC(processor.ProcessorABC, ABC): If Subsamples are defined a columnsmager is created for each of them. ''' self.column_managers = {} - for subs in self._subsamples_names: + for subs in self._subsamples[self._sample].keys(): self.column_managers[subs] = ColumnsManager( self.cfg.columns[subs], subs, self._categories ) @@ -483,14 +451,12 @@ class BaseProcessorABC(processor.ProcessorABC, ABC): # TODO Fill column accumulator for different variations if self._hasSubsamples: # call the filling for each - for subs in self._subsamples_names: - # Get the mask - subsam_mask = self._subsamples_masks.all(*self._subsamples_map[subs]) + for subs in self._subsamples[self._sample].keys(): # Calling hist manager with a subsample mask self.output["columns"][subs] = self.column_managers[subs].fill_columns( self.events, self._categories, - subsample_mask=subsam_mask, + subsample_mask=self._subsamples[self._sample].get_mask(subs), ) else: self.output["columns"][self._sample] = self.column_managers[ @@ -658,7 +624,7 @@ class BaseProcessorABC(processor.ProcessorABC, ABC): # This is computed before any preselection self.output['sum_genweights'][self._sample] = ak.sum(self.events.genWeight) if self._hasSubsamples: - for subs in self._subsamples_names: + for subs in self._subsamples[self._sample].keys(): self.output['sum_genweights'][subs] = self.output['sum_genweights'][ self._sample ] -- GitLab From 4fe74a09ba85ad2a14bec6913e203e2f4c07c3a5 Mon Sep 17 00:00:00 2001 From: mmarchegiani <matmarcheg@gmail.com> Date: Fri, 24 Mar 2023 03:24:02 +0100 Subject: [PATCH 11/11] Undo changes inherited from other branches --- .../parameters/object_preselection.py | 3 + pocket_coffea/utils/plot_utils.py | 81 ++---------------- pocket_coffea/workflows/genweights.py | 57 +------------ scripts/README.md | 85 ------------------- scripts/hadd_skimmed_files.py | 9 +- 5 files changed, 12 insertions(+), 223 deletions(-) delete mode 100644 scripts/README.md diff --git a/pocket_coffea/parameters/object_preselection.py b/pocket_coffea/parameters/object_preselection.py index 291f816d..e93c0453 100644 --- a/pocket_coffea/parameters/object_preselection.py +++ b/pocket_coffea/parameters/object_preselection.py @@ -89,6 +89,9 @@ object_preselection = { "eta": 2.4, "msd": 40, "jetId": 2, + "nsubjet": 2, + "nmusj": 1, + "dimuon_pt_ratio": 0.6, }, }, } diff --git a/pocket_coffea/utils/plot_utils.py b/pocket_coffea/utils/plot_utils.py index 46d3081d..00c0d4ab 100644 --- a/pocket_coffea/utils/plot_utils.py +++ b/pocket_coffea/utils/plot_utils.py @@ -468,9 +468,6 @@ def plot_data_mc_hist1D( samples = h.keys() samples_data = list(filter(lambda d: 'DATA' in d, samples)) samples_mc = list(filter(lambda d: 'DATA' not in d, samples)) - - print("********************") - print(samples_mc) h_mc = h[samples_mc[0]] @@ -490,8 +487,6 @@ def plot_data_mc_hist1D( for year in years: for cat in categories: - if 'noreweight' in cat: - continue if only_cat: if isinstance(only_cat, list): if not cat in only_cat: @@ -571,7 +566,7 @@ def plot_data_mc_hist1D( [ h[d][slicing_mc] for d in samples_mc - if (d.endswith(f'_{f}') & ('GluGluH' not in d)) + if d.endswith(f'_{f}') ] ) ) @@ -583,7 +578,7 @@ def plot_data_mc_hist1D( [ h[d][slicing_mc_nominal] for d in samples_mc - if (d.endswith(f'_{f}') & ('GluGluH' not in d)) + if d.endswith(f'_{f}') ] ) ) @@ -593,57 +588,6 @@ def plot_data_mc_hist1D( nevents = { f: round(sum(dict_mc_nominal[f].values()), 1) for f in flavors } - print("************") - print(histname) - print(samples_mc) - dict_gghbb = { - f: stack_sum( - hist.Stack.from_iter( - [ - h[d][slicing_mc] - for d in samples_mc - if (d.endswith(f'_{f}') & ('GluGluHToBB' in d)) - ] - ) - ) - for f in flavors - } - dict_gghbb_nominal = { - f: stack_sum( - hist.Stack.from_iter( - [ - h[d][slicing_mc_nominal] - for d in samples_mc - if (d.endswith(f'_{f}') & ('GluGluHToBB' in d)) - ] - ) - ) - for f in flavors - } - dict_gghcc = { - f: stack_sum( - hist.Stack.from_iter( - [ - h[d][slicing_mc] - for d in samples_mc - if (d.endswith(f'_{f}') & ('GluGluHToCC' in d)) - ] - ) - ) - for f in flavors - } - dict_gghcc_nominal = { - f: stack_sum( - hist.Stack.from_iter( - [ - h[d][slicing_mc_nominal] - for d in samples_mc - if (d.endswith(f'_{f}') & ('GluGluHToCC' in d)) - ] - ) - ) - for f in flavors - } else: dict_mc = {d: h[d][slicing_mc] for d in samples_mc} dict_mc_nominal = {d: h[d][slicing_mc_nominal] for d in samples_mc} @@ -672,15 +616,6 @@ def plot_data_mc_hist1D( dict_mc_nominal[sample] = histo_reweighted stack_mc = hist.Stack.from_dict(dict_mc) stack_mc_nominal = hist.Stack.from_dict(dict_mc_nominal) - if flavorsplit == '5f': - stack_gghbb = hist.Stack.from_dict(dict_gghbb) - stack_gghbb_nominal = hist.Stack.from_dict(dict_gghbb_nominal) - stack_gghcc = hist.Stack.from_dict(dict_gghcc) - stack_gghcc_nominal = hist.Stack.from_dict(dict_gghcc_nominal) - hist_gghbb = stack_sum(stack_gghbb) - hist_gghbb_nominal = stack_sum(stack_gghbb_nominal) - hist_gghcc = stack_sum(stack_gghcc) - hist_gghcc_nominal = stack_sum(stack_gghcc_nominal) if not is_mc_only: # Sum over eras if era axis exists in data histogram @@ -731,12 +666,6 @@ def plot_data_mc_hist1D( x, h_data.values(), yerr=np.sqrt(h_data.values()), **opts_data ) # stack_data.plot(stack=True, color='black', ax=ax) - if flavorsplit == '5f': - sf = 100 - hist_gghbb_nominal = sf * hist_gghbb_nominal - hist_gghcc_nominal = sf * hist_gghcc_nominal - hist_gghbb_nominal.plot(stack=False, histtype='step', ax=ax, color='red', label=f'ggHbb x {sf}') - hist_gghcc_nominal.plot(stack=False, histtype='step', ax=ax, color='green', label=f'ggHcc x {sf}') if rebinning: syst_err2_up, syst_err2_down = get_systematic_uncertainty( stack_mc, @@ -799,7 +728,7 @@ def plot_data_mc_hist1D( exp = math.floor( math.log(max(stack_sum(stack_mc_nominal).values()), 10) ) - ax.set_ylim((0.01, 10 ** (exp + 3))) + ax.set_ylim((0.01, 10 ** (exp + 2))) # if flavorsplit == '5f': # ax.legend(handles, labels, loc="upper right", fontsize=fontsize, ncols=2) # else: @@ -818,7 +747,7 @@ def plot_data_mc_hist1D( handles_new.append(handles[i]) labels = labels_new handles = handles_new - ax.legend(handles, labels, fontsize=fontsize, ncols=2, loc="upper right") + ax.legend(handles, labels, fontsize=fontsize, ncols=2) if ("variables" in config.plot_options) & (histname in config.plot_options["variables"].keys()): cfg_plot = config.plot_options["variables"][histname] if 'xlim' in cfg_plot.keys(): @@ -831,7 +760,7 @@ def plot_data_mc_hist1D( max(stack_sum(stack_mc_nominal).values()), 10 ) ) - ax.set_ylim((0.01, 10 ** (exp + 3))) + ax.set_ylim((0.01, 10 ** (exp + 2))) # ax.legend(handles, labels, loc="upper right", fontsize=fontsize, ncols=2) if 'ylim' in cfg_plot.keys(): if isinstance(cfg_plot['ylim'], tuple): diff --git a/pocket_coffea/workflows/genweights.py b/pocket_coffea/workflows/genweights.py index dc0ffab7..1a404377 100644 --- a/pocket_coffea/workflows/genweights.py +++ b/pocket_coffea/workflows/genweights.py @@ -1,16 +1,9 @@ -import numpy as np import awkward as ak import copy -from coffea.lumi_tools import LumiMask -from coffea.analysis_tools import PackedSelection - from .base import BaseProcessorABC from ..utils.configurator import Configurator -from ..parameters.event_flags import event_flags, event_flags_data -from ..parameters.lumi import goldenJSON -from ..parameters.btag import btag class genWeightsProcessor(BaseProcessorABC): def __init__(self, cfg: Configurator) -> None: @@ -18,8 +11,7 @@ class genWeightsProcessor(BaseProcessorABC): self.output_format = { "sum_genweights": {}, "cutflow": { - "initial": {s: 0 for s in self.cfg.datasets}, - "skim": {s: 0 for s in self.cfg.datasets} + "initial": {s: 0 for s in self.cfg.datasets} } } @@ -27,14 +19,9 @@ class genWeightsProcessor(BaseProcessorABC): self._dataset = self.events.metadata["dataset"] self._sample = self.events.metadata["sample"] self._year = self.events.metadata["year"] - self._btag = btag[self._year] self._isMC = self.events.metadata["isMC"] == "True" if self._isMC: - self._era = "MC" self._xsec = self.events.metadata["xsec"] - else: - self._era = self.events.metadata["era"] - self._goldenJSON = goldenJSON[self._year] # Check if the user specified any subsamples without performing any operation if self._sample in self._subsamplesCfg: @@ -42,39 +29,6 @@ class genWeightsProcessor(BaseProcessorABC): else: self._hasSubsamples = False - def skim_events(self): - self._skim_masks = PackedSelection() - mask_flags = np.ones(self.nEvents_initial, dtype=np.bool) - flags = event_flags[self._year] - if not self._isMC: - flags += event_flags_data[self._year] - for flag in flags: - mask_flags = getattr(self.events.Flag, flag) - self._skim_masks.add("event_flags", mask_flags) - - self._skim_masks.add("PVgood", self.events.PV.npvsGood > 0) - - if not self._isMC: - mask_lumi = LumiMask(self._goldenJSON)( - self.events.run, self.events.luminosityBlock - ) - self._skim_masks.add("lumi_golden", mask_lumi) - - for skim_func in self._skim: - mask = skim_func.get_mask( - self.events, - year=self._year, - sample=self._sample, - btag=self._btag, - isMC=self._isMC, - ) - self._skim_masks.add(skim_func.id, mask) - - self.events = self.events[self._skim_masks.all(*self._skim_masks.names)] - self.nEvents_after_skim = self.nevents - self.output['cutflow']['skim'][self._dataset] += self.nEvents_after_skim - self.has_events = self.nEvents_after_skim > 0 - def apply_object_preselection(self, variation): pass @@ -94,13 +48,4 @@ class genWeightsProcessor(BaseProcessorABC): if self._hasSubsamples: raise Exception("This processor cannot compute the sum of genweights of subsamples.") - self.process_extra_before_skim() - self.skim_events() - if not self.has_events: - return self.output - - if self.cfg.save_skimmed_files: - self.export_skimmed_chunk() - return self.output - return self.output diff --git a/scripts/README.md b/scripts/README.md deleted file mode 100644 index 940cc5a1..00000000 --- a/scripts/README.md +++ /dev/null @@ -1,85 +0,0 @@ -# PocketCoffea scripts - -The current guide is a documentation of how to run the scripts of the PocketCoffea package. - -## Merge skimmed files with hadd - -In order to merge the skimmed n-tuples, the `hadd` command from ROOT needs to be used. -To source ROOT inside the `pocket_coffea` environment, run the following command: - -``` -source /cvmfs/sft.cern.ch/lcg/views/LCG_102/x86_64-centos7-gcc11-opt/setup.sh -``` - -To run the `hadd` script run the following command: - -``` -cd /path/to/PocketCoffea -python scripts/hadd_skimmed_files.py -fl /path/to/output/output_all.coffea -o /path/to/outputdir -``` - -## Build dataset - -The first step is to build the json files containing the list of files of the UL datasets to process. In our case, we need to include the `TTToSemiLeptonic` and `TTTo2L2Nu` MC datasets and `SingleMuon` as data dataset. - -A step-by-step detailed guide on the creation of the dataset file can be found at this [link](https://pocketcoffea.readthedocs.io/en/latest/examples.html). - -### Append sum of genweights to datasets definitions - -In order to update the `datasets_definitions.json` and include the `sum_genweights` metadata in the datasets definitions, a dedicated script needs to be run with the following command: - -``` -cd /path/to/PocketCoffea -scripts/dataset/append_genweights.py --cfg dataset/dataset_definitions.json -i output/genweights/genweights_2016_PreVFP/output_all.coffea --overwrite -scripts/dataset/append_genweights.py --cfg dataset/dataset_definitions.json -i output/genweights/genweights_2016_PostVFP/output_all.coffea --overwrite -``` - -N.B.: the argument `-i output/genweights/genweights_201*/output_all.coffea` should be the path to the Coffea output file where the genweights are stored. The option `--overwrite` automatically updates the config file passed as `--cfg` argument with the genweights, in this case `dataset/dataset_definitions.json`. - -### Build dataset with sum of genweights - -Now the steps of the [Build dataset](#build-dataset) section need to be repeated to generate the json datasets with the additional metadata `sum_genweights`. - -## Run the analysis - -In order to run the analysis workflow and produce the output histograms, run the following command: - -``` -cd /path/to/PocketCoffea -runner.py --cfg /path/to/config/config.py --full -``` - -N.B.: the argument `--full` will process all the datasets together at once and save the output in a single output file, `output_all.coffea`. Otherwise, the datasets are processed separately and an output file is saved for each dataset. - -## Accumulate output files - -If the output of the [previous step](#run-the-analysis) has been produced without the argument `--full`, the output files need to be merged in a single output file `output_all.coffea`. If the output has been produced with the argument `--full` and the output file `output_all.coffea` is already existing, skip this step and continue with the [next one](#produce-datamc-plots). - -Once the Coffea output files are produced, one needs to merge the files into a single file by using the script `accumulate_files.py` by running this command: - -``` -cd /path/to/output -python accumulate_files.py -i output1.coffea output2.coffea output3.coffea -o output_all.coffea -``` - -## Produce data/MC plots - -To produce plots from a Coffea output file execute the plotting script: - -``` -cd /path/to/PocketCoffea -make_plots.py --cfg /path/to/config/config.py -i /path/to/output/output_all.coffea -``` - -## Run trigger SF script - -To run the script that computes the single electron trigger SF and produces the validation plots, run: - -``` -cd /path/to/PocketCoffea -python scripts/plot/trigger_efficiency.py --cfg config/semileptonic_triggerSF/semileptonic_triggerSF_2017.py --save_plots -``` - -If the additional argument `--save_plots` is not specified, only the scale factor maps are saved without saving the plots. - -The output plots are saved in `/path/to/output/plots/trigger_efficiency` and `/path/to/output/plots/trigger_scalefactor`, while the 1D and 2D scale factor maps are saved in the folder specified in `workflow_options['output_triggerSF']` in the config file. diff --git a/scripts/hadd_skimmed_files.py b/scripts/hadd_skimmed_files.py index 6d870ba2..2280fcca 100644 --- a/scripts/hadd_skimmed_files.py +++ b/scripts/hadd_skimmed_files.py @@ -18,7 +18,7 @@ parser.add_argument( ) parser.add_argument("-o", "--outputdir", required=False, type=str, help="Output folder") parser.add_argument( - "--only-datasets", nargs="+", type=str, help="Restricting list of datasets" + "--only-datasets", nargs="+", type=str, help="Restring list of datasets" ) parser.add_argument("-f", "--files", type=int, help="Limit number of files") parser.add_argument("-e", "--events", type=int, help="Limit number of files") @@ -43,9 +43,6 @@ groups_metadata = {} for dataset in df["skimmed_files"].keys(): if args.only_datasets and dataset not in args.only_datasets: continue - outputdir_dataset = os.path.join(args.outputdir, dataset) - if not os.path.exists(outputdir_dataset): - os.makedirs(outputdir_dataset) groups_metadata[dataset] = defaultdict(dict) nevents_tot = 0 nfiles = 0 @@ -57,7 +54,7 @@ for dataset in df["skimmed_files"].keys(): if (args.files and (nfiles + 1) > args.files) or ( args.events and (nevents_tot + nevents) > args.events ): - outfile = os.path.join(outputdir_dataset, f"{dataset}_{ngroup}.root") + outfile = f"{args.outputdir}/{dataset}/{dataset}_{ngroup}.root" workload.append((outfile, group[:])) groups_metadata[dataset]["files"][outfile] = group[:] group.clear() @@ -71,7 +68,7 @@ for dataset in df["skimmed_files"].keys(): # add last group if len(group): - outfile = os.path.join(outputdir_dataset, f"{dataset}_{ngroup}.root") + outfile = f"{args.outputdir}/{dataset}/{dataset}_{ngroup}.root" workload.append((outfile, group[:])) groups_metadata[dataset]["files"][outfile] = group[:] -- GitLab