You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

242 lines
179 KiB

10 months ago
  1. {
  2. "cells": [
  3. {
  4. "cell_type": "code",
  5. "execution_count": 1,
  6. "metadata": {},
  7. "outputs": [],
  8. "source": [
  9. "import uproot\n",
  10. "import awkward as ak\n",
  11. "import matplotlib\n",
  12. "import matplotlib.pyplot as plt\n",
  13. "input_tree = uproot.open({\"/work/guenther/reco_tuner/data/param_data_selected_all_p.root\": \"Selected\"})\n",
  14. "array = input_tree.arrays()\n",
  15. "array = array[[field for field in ak.fields(array) if \"scifi_hit\" not in field]]\n",
  16. "df = ak.to_pandas(array)"
  17. ]
  18. },
  19. {
  20. "cell_type": "code",
  21. "execution_count": 2,
  22. "metadata": {},
  23. "outputs": [
  24. {
  25. "data": {
  26. "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuEAAAJQCAYAAAAg1rpCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB0JklEQVR4nO3dd3xUVf7/8fdAKgmEQOiEIhAINhCCIrAYRBBRQVFUumJbo6LiqqCC+LWhiKIbKy5FlKIuqKCsrLQgLSJioQhSAoJCaCGkQJLz+4Nf7maSmWSSTEvyej4eeTxm5p577uee3Jl8cubcc2zGGCMAAAAAXlPN1wEAAAAAVQ1JOAAAAOBlJOEAAACAl5GEAwAAAF5GEg4AAAB4GUk4AAAA4GUk4QAAAICXkYQDAAAAXkYSDgAAAHgZSThQxc2cOVM2m002m02jRo1yW7179+616m3RooXb6gXc4dNPP9X111+vJk2aKDg42LpWr7jiCl+HhkqmRYsW1vW1d+9eX4cDPxLg6wAAAPCmESNG6MMPP/R1GACqOJJwACXau3evWrZsKUlq3rw5vTmosObOnWuXgHfp0kXt27dXWFiYJKlNmza+Cg1AFUMSDgCoMmbNmmU9njRpkiZMmODDaABUZTZjjPF1EAD8Gz3hqCzq16+vI0eOSJIOHDigJk2a+DgiVHYtWrTQvn37JEl79uzhHhlYuDETAFBlHD9+3HrcqFEjH0YCoKojCQcAVBk5OTnW42rV+BMIwHf4BIJf2r17t2rVqmVN6zR16tQS97n33nvtpsQ7efJkmY8/atQoq66ZM2dKko4dO6aXX35ZXbp0Ub169RQaGqrzzjtPd955pzZt2lSq+s+ePasZM2Zo4MCBat68uUJDQ1WrVi21bdtWo0eP1rJly1yuKzU1VVOmTFHv3r3VuHFjhYSEqEaNGmrevLk6deqk2267TTNmzNAff/zhcP/ipijM35Y/FEWS9u3bZ5Uv/FNQWaYoXL9+ve6//36df/75ioyMVEhIiJo2baqrr75a//znP3X69OkS63jmmWes4z7zzDOSziVes2fPVu/eva0p6Ro1aqSBAwdq8eLFLsUmScuXL9fo0aN14YUXqnbt2goMDFRUVJRiY2N1xRVX6Mknn9Tq1at19uxZl+ssibuuFUfTpB04cEBPP/20Lr74YtWuXVthYWFq166dHnjgAevrc3dwdi2sWbNGd955p9q1a6eIiAjZbDY99NBDDuvYv3+//u///k89evRQ48aNFRwcrDp16qhjx4569NFH9dtvv7l07gUVvn6Lu05Pnz6tt99+W9ddd52aN2+uGjVqqGbNmmrTpo3uuOMOLV++vMR2cPRey83N1bx58zRgwACdd955Cg0Nlc1m06JFixzWkZycrIcfflgdOnRQvXr1FBQUpIYNG6pnz56aPHmyXU+/K+3hiWshNzdXCxYs0IgRI9S2bVtFRkYqMDBQdevW1WWXXaYxY8bo22+/lSujYd1xvq5auXKlw6kqFy1apAEDBqhFixYKCQlRw4YNddVVV2nmzJnKy8tz2/FPnjypuXPn6p577tGll16qqKgoBQUFqVatWmrdurWGDBmiTz75xKVjOvtcX7hwoa677jo1a9ZMwcHBql+/vvr06aM5c+a49PsoaNu2bRo/fry6dOmiBg0aKCgoSPXq1dOll16qCRMm6ODBg6VtgqrHAH5q9uzZRpKRZIKCgszmzZudll20aJFVtlq1amb16tXlOvbIkSOt+mbMmGHWrVtnmjRpYr1W+KdatWpm3LhxLtW9fv1606pVK6d15f9cddVV5vDhw8XWtWjRIhMZGVliXZJMkyZNHNYxY8YMq8zIkSOdbnPlp6A9e/ZYrzdv3rzY80hPTze33HJLifU3atTIfPXVV8XWNXHiRKv8xIkTzYEDB8zll19ebL233367yc3NdVrnqVOnzPXXX+9yO7z//vvFxugqd14rzZs3t8rv2bPHLFy40ERERDitMzQ01CxevNgt51H4WsjOzjb33nuvw+OOGTPGbt/c3Fzz9NNPm5CQkGLbICAgwIwfP97k5eUVe+7F/Ti7ThcsWGAaNmxY4v7XXnutOXHihNN2KPxe++OPP0yPHj0c1rVw4UK7fY8dO2YGDRpUYgy1a9c2n3zySbG/D09eC6tXrzYxMTEutffjjz/utB53nq+rVqxYYdXbs2dPk5aWZgYOHFjs8ePi4swff/xRbL2F29uRzz77zAQHB7vUbh06dHBaT77C19qJEydK/Ay7+uqrTUZGRontlJWVZe69915TvXr1YusLDQ01b775Zon1VWXMjgK/NXz4cH399deaO3euzpw5oyFDhmjTpk0KDQ21K3fw4EGNHj3aej5+/Hj16NHDbXHs27dPjzzyiI4fP66wsDD16tVLDRo00J9//qkVK1bo9OnTysvL04svvqicnBy9/PLLTutavXq1+vXrp4yMDOu1Ll266Pzzz9eZM2e0fv16/f7775KkZcuWqXv37lqzZo3q1atXpK7vv/9eN910k/X1emhoqC677DK1aNFCwcHBSktL0++//66ff/7Z7nilERsbq4SEBJ06dUqzZ8+WJNWsWVMjRowoU32OZGRkqFevXtq4caP1WuPGjdWjRw+Fh4dr165dWrNmjXJzc3Xo0CFdf/31mjt3rm666aYS605PT9fVV1+tX375RTVq1FCPHj0UHR2tU6dOacWKFTp8+LAkacaMGWrbtq0ef/xxh/UMHz5cX3zxhfW8devW6tixo+rUqaOzZ8/qyJEj+vnnn916w6o7r5XCvv32W91zzz3Kzc1Vs2bN1LVrV9WqVUt79uzRypUrlZOTo8zMTA0ePFi//PKL3Tch7vDwww/rnXfekSRdeOGFuvjiixUYGKjffvvNbohIbm6ubrnlFn322WfWa40aNdKll16q+vXrKz09XRs2bNDvv/+unJwcvfDCCzpy5Ijee+89u+ONHDlSR48elSQlJiZaryckJNiVq1u3bpFYX3vtNY0dO9bqJaxZs6a6du2q6Oho5ebmauvWrUpOTpYxRosXL1bPnj21du1a1ahRo9g2yM7O1vXXX69NmzYpICBAl19+uVq3bq2srCz98MMPdmX//PNP9erVS9u2bbNei42NVYcOHVSzZk0dPnxYa9asUWpqqk6cOKHBgwfrww8/1NChQ4uNQXLvtTBv3jyNGDHC7pugmJgYXXLJJYqIiNDJkyf166+/6tdff1VeXp6ysrIc1uPJ8y2N22+/3fpGIv+9l52drfXr12v37t2SzvXU9+rVS2vXrlWdOnXKfKzDhw8rOztbktS0aVO1b99eDRs2VI0aNZSenq5t27bphx9+kDFGP/74o3r06KEff/zR4TVbWG5urgYNGqRvv/1WQUFBuvzyy9WqVStlZWUpKSlJKSkpkqSlS5fqkUce0dtvv+20rtOnT6tv37767rvvrNdatmypzp07KzIyUsePH9fatWv1xx9/KDMzUw888IDS0tI0fvz4MrdNpebjfwKAYp04ccKuF+Hee++1256Xl2d69+5tbb/00kvN2bNny33cgj3hQUFBRpIZMmRIkV6uEydOmFtvvdXuv//ly5c7rPPYsWN2vemtWrUyycnJRcrNmTPHhIaGWuWuu+46h/UNGDDAKjNo0CBz7Ngxh+WysrLMkiVLzD333ONwe3E94flK06td2n3+/ve/W+WqV69upk6dWqRX+rfffjOdOnWyytWqVcvs3r3bYX0Fe8Lze5ZGjhxpjh49alfu9OnT5rbbbrPKhoeHm/T09CL1bd682a5McT3xv//+u3nuuefMF198UUzLlMzd14ox9r1xwcHBJiwszHz44YdFeo5/+eUXu2Pffvvt5ToXY+yvhfzes+joaIffWGVlZVmPn376aWu/+vXrm/nz5zv8xuKTTz6x68mdP3++01gKvldL8t///tdUq1bNSDKBgYHmueeec3qNtG/f3qr373//u8P6Cr7XAgICrB5XR72a+e2Qm5tr4uPjrf0uueQSh9dCZmameeaZZ4zNZjOSTFhYmNP3iCeuhR9++MHu24qOHTua9evXOyx76NA
  27. "text/plain": [
  28. "<Figure size 800x600 with 1 Axes>"
  29. ]
  30. },
  31. "metadata": {},
  32. "output_type": "display_data"
  33. }
  34. ],
  35. "source": [
  36. "urax = df.loc[(df[\"x_ref\"] > -3200) & (df[\"x_ref\"] < 3200), \"x_ref\"].plot.hist(bins=100, logy=True,title='x positions on reference plane', figsize=(8,6), fontsize=20,edgecolor='black')\n",
  37. "urax.set_xlabel('x [mm]', fontsize=24)\n",
  38. "urax.set_ylabel('Count', fontsize=24)\n",
  39. "urax.title.set_size(24)"
  40. ]
  41. },
  42. {
  43. "cell_type": "code",
  44. "execution_count": 3,
  45. "metadata": {},
  46. "outputs": [],
  47. "source": [
  48. "from scipy.optimize import curve_fit\n",
  49. "import numpy as np\n",
  50. "def fitFastSigmoid(x, p0, p1, p2 ):\n",
  51. " return p0 + p1 * x / (1 + abs(x * p2) ) \n",
  52. "\n",
  53. "def sqrtSigmoid( x, p0, p1, p2 ):\n",
  54. " return p0 + p1 * x / np.sqrt( p2 + x**2 )"
  55. ]
  56. },
  57. {
  58. "cell_type": "code",
  59. "execution_count": 14,
  60. "metadata": {},
  61. "outputs": [
  62. {
  63. "name": "stdout",
  64. "output_type": "stream",
  65. "text": [
  66. "-1152.3191321144784 1144.5112688972854\n",
  67. "BINS=900\n",
  68. "shift=126.0\n",
  69. "[ 4.50871394e+02 5.78026621e-01 -6.72848459e-04]\n",
  70. "at -3000 = 2\n",
  71. "at 0 = 576\n",
  72. "at 3000 = 1151\n",
  73. "at -3000 = 2.3960602381498575\n",
  74. "at 0 = 576.8713937732111\n",
  75. "at 3000 = 1151.3467273082724\n",
  76. "\n",
  77. "-1152.3191321144784 1145.2265343593158\n",
  78. "BINS=1000\n",
  79. "shift=140.0\n",
  80. "[5.00965281e+02 6.42244329e-01 6.72824292e-04]\n",
  81. "at -3000 = 2\n",
  82. "at 0 = 640\n",
  83. "at 3000 = 1279\n",
  84. "at -3000 = 2.6514499904290005\n",
  85. "at 0 = 640.965281389201\n",
  86. "at 3000 = 1279.279112787973\n",
  87. "\n"
  88. ]
  89. },
  90. {
  91. "data": {
  92. "text/html": [
  93. "<div>\n",
  94. "<style scoped>\n",
  95. " .dataframe tbody tr th:only-of-type {\n",
  96. " vertical-align: middle;\n",
  97. " }\n",
  98. "\n",
  99. " .dataframe tbody tr th {\n",
  100. " vertical-align: top;\n",
  101. " }\n",
  102. "\n",
  103. " .dataframe thead th {\n",
  104. " text-align: right;\n",
  105. " }\n",
  106. "</style>\n",
  107. "<table border=\"1\" class=\"dataframe\">\n",
  108. " <thead>\n",
  109. " <tr style=\"text-align: right;\">\n",
  110. " <th></th>\n",
  111. " <th>NBINS</th>\n",
  112. " <th>nBinsReal</th>\n",
  113. " <th>minBinWidth</th>\n",
  114. " <th>maxBinWidth</th>\n",
  115. " </tr>\n",
  116. " </thead>\n",
  117. " <tbody>\n",
  118. " <tr>\n",
  119. " <th>0</th>\n",
  120. " <td>900</td>\n",
  121. " <td>1152</td>\n",
  122. " <td>1.520524</td>\n",
  123. " <td>411.883522</td>\n",
  124. " </tr>\n",
  125. " <tr>\n",
  126. " <th>1</th>\n",
  127. " <td>1000</td>\n",
  128. " <td>1280</td>\n",
  129. " <td>1.378500</td>\n",
  130. " <td>406.784548</td>\n",
  131. " </tr>\n",
  132. " </tbody>\n",
  133. "</table>\n",
  134. "</div>"
  135. ],
  136. "text/plain": [
  137. " NBINS nBinsReal minBinWidth maxBinWidth\n",
  138. "0 900 1152 1.520524 411.883522\n",
  139. "1 1000 1280 1.378500 406.784548"
  140. ]
  141. },
  142. "execution_count": 14,
  143. "metadata": {},
  144. "output_type": "execute_result"
  145. },
  146. {
  147. "data": {
  148. "image/png": "iVBORw0KGgoAAAANSUhEUgAAB5AAAAdnCAYAAAADTeAUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd3hUVdf38e+kF5JAgDR670VAIRRDD0izgqARRFEeQUUsgOANYqGIFEUUARFEikoRLPQuLVKkVwMIBClCQgmp+/ljzMiQQoDAJPD7XFfe9z777NlnnTPjPKxZ++xjMcYYRERERERERERERERERETknufk6ABERERERERERERERERERCRnUAFZREREREREREREREREREQAFZBFRERERERERERERERERORfKiCLiIiIiIiIiIiIiIiIiAigArKIiIiIiIiIiIiIiIiIiPxLBWQREREREREREREREREREQFUQBYRERERERERERERERERkX+pgCwiIiIiIiIiIiIiIiIiIoAKyCIiIiIiIiIiIiIiIiIi8i8VkEVEREREREREREREREREBFABWUREREREROSesnr1atq0aUNISAgWi4V58+bZ9iUmJtKnTx+qVKmCt7c3ISEhPPPMM5w4ccJujPj4eF5++WUKFCiAt7c3bdu25dixY3Z9zp07R0REBH5+fvj5+REREcH58+fvwBmKiIiIiIjIrXBxdAC3S0pKCidOnMDHxweLxeLocERERERERHIcYwwXLlwgJCQEJyfNL75XXLp0iWrVqvHss8/y2GOP2e27fPkyW7Zs4Z133qFatWqcO3eOXr160bZtW37//Xdbv169erFgwQJmzpxJ/vz5ef3112ndujWbN2/G2dkZgE6dOnHs2DEWLlwIwAsvvEBERAQLFizIcqzK7UVERERERDJ2u/J6izHGZNtoOcixY8coUqSIo8MQERERERHJ8f766y8KFy7s6DDEASwWC3PnzuXhhx/OsE9kZCQPPPAAR44coWjRosTExFCwYEG++eYbOnToAMCJEycoUqQIv/zyC+Hh4ezZs4eKFSuyYcMGateuDcCGDRsIDQ1l7969lCtXLkvxKbcXERERERG5vuzO6+/aO5B9fHwA6wXz9fV1cDQiIiIiIiI5T2xsLEWKFLHlTyLpiYmJwWKxkDdvXgA2b95MYmIizZs3t/UJCQmhcuXKrFu3jvDwcNavX4+fn5+teAxQp04d/Pz8WLduXYYF5Pj4eOLj423bqXPelduLiIiIiIikdbvy+ru2gJy6tJWvr6+STBERERERkUxoaWDJyJUrV+jbty+dOnWy5dYnT57Ezc2NfPny2fUNDAzk5MmTtj4BAQFpxgsICLD1Sc+QIUN4991307QrtxcREREREclYduf1esiViIiIiIiIiKSRmJjIk08+SUpKCuPGjbtuf2OM3Y8W6f2AcW2fa/Xr14+YmBjb319//XVzwYuIiIiIiMhNUwFZREREREREROwkJibSvn17oqKiWLJkid3dv0FBQSQkJHDu3Dm715w6dYrAwEBbn7///jvNuKdPn7b1SY+7u7vtbmPddSwiIiIiIuIYKiCLiIiIiIiIiE1q8fjAgQMsXbqU/Pnz2+2vWbMmrq6uLFmyxNYWHR3Nzp07qVu3LgChoaHExMSwadMmW5+NGzcSExNj6yMiIiIiIiI50137DOSsSk5OJjEx0dFhiEgO4+bmhpOT5tiIiIiIyN3n4sWLHDx40LYdFRXFtm3b8Pf3JyQkhMcff5wtW7bw008/kZycbHtmsb+/P25ubvj5+fHcc8/x+uuvkz9/fvz9/XnjjTeoUqUKTZs2BaBChQq0aNGCbt26MX78eABeeOEFWrduTbly5bL1fFJSUkhISMjWMUUkd3F1dcXZ2dnRYYiIiIjcNe7ZArIxhpMnT3L+/HlHhyIiOZCTkxMlSpTAzc3N0aGIiIiIiGSr33//nUaNGtm2e/fuDUDnzp0ZNGgQ8+fPB6B69ep2r1uxYgUNGzYEYNSoUbi4uNC+fXvi4uJo0qQJX3/9tV0B59tvv+WVV16hefPmALRt25axY8dm67kkJCQQFRVFSkpKto4rIrlP3rx5CQoKyvQ56yIiIiKSNRZjjHF0ELdDbGwsfn5+xMTEpPvMpOjoaM6fP09AQABeXl76x6WI2KSkpHDixAlcXV0pWrSovh9ERETkrnW9vEnE0TL7jBpjOHr0KImJiYSEhGgFIZF7lDGGy5cvc+rUKfLmzUtwcLCjQxIRERG5Y25XXn9P3oGcnJxsKx5f+ywnERGAggULcuLECZKSknB1dXV0OCIiIiIico2kpCQuX75MSEgIXl5ejg5HRBzI09MTgFOnThEQEKDlrEVERERu0T05PTf1mcdKMEUkI6lLVycnJzs4EhERERERSU/qv9X12BkRgf9+50v93U9EREREbt49WUBOpWVpRSQj+n4QEREREckd9G93EQF9F4iIiIhkp3u6gCwiIiIiIiIiIiIiIiIiIv9RAfku0LBhQ3r16uXoMERERERERETkJim3FxERERGRnEIF5HvMypUrsVgsnD9/3tGhiIiIiIiIiMhNUG4vIiIiIiK3kwrIIiIiIiIiIiIiIiIiIiICqICcLaJj4lh36AzRMXG3/ViXLl3imWeeIU+ePAQHB/Pxxx/b7Z82bRq1atXCx8eHoKAgOnXqxKlTpwA4fPgwjRo1AiBfvnxYLBa6dOkCwMKFC6lfvz558+Ylf/78tG7dmkOHDt328xERERERERFxtDuZ14NyexERERERydlUQL5FsyKPUm/ocjpN2Ei9ocuZFXn0th7vzTffZMWKFcydO5fFixezcuVKNm/ebNufkJDAe++9xx9//MG8efOIioqyJZJFihRh9uzZAOzbt4/o6GjGjBkDWJPX3r17ExkZybJly3BycuKRRx4hJSXltp6PiIiIiIiIiCPd6bwelNuLiIiIiEjO5uLoAHKz6Jg4+s3ZQYqxbqcYeHvOTh4sW5BgP89sP97FixeZNGkSU6dOpVmzZgBMmTKFwoUL2/p07drV9r9LlizJJ598wgMPPMDFixfJkycP/v7+AAQEBJA3b15b38cee8zuWJMmTSIgIIDdu3dTuXLlbD8XEREREREREUe703k9KLcXEREREZGcT3cg34KoM5dsSWaqZGM4fObybTneoUOHSEhIIDQ01Nbm7+9PuXLlbNtbt26lXbt2FCtWDB8fHxo2bAjA0aOZz6A+dOgQnTp1omTJkvj6+lKiRIksvU5EREREREQkt7rTeT0otxcRERERkZxPBeRbUKKAN04W+zZni4XiBbxuy/GMMZnuv3TpEs2bNydPnjxMmzaNyMhI5s6dC1iXv8pMmzZtOHv2LBMmTGDjxo1s3LgxS68TERERERERya3udF4Pyu1FRERERCTnUwH5FgT7eTLk0So4W6zZprPFwoePVr5ty1yVLl0aV1dXNmzYYGs7d+4c+/fvB2Dv3r2cOXOGoUOH0qBBA8qXL8+pU6fsxnBzcwMgOTnZ1nb27Fn27NnDgAEDaNKkCRUqVODcuXO35RxEREREREREcoo7ndeDcnsREREREcn59AzkW9Th/qI8WLYgh89cpngBr9uaZObJk4fnnnuON998k/z58xMYGEj//v1xcrLOAyhatChubm58+umndO/enZ07d/Lee+/ZjVGsWDEsFgs//fQTDz30EJ6enuTLl4/8+fPz5ZdfEhwczNGjR+nbt+9tOw8RERERERGRnOJO5vWg3F5ERERERHI+3YGcDYL9PAktlf+2J5kAH330EQ8++CBt27aladOm1K9fn5o1awJQsGBBvv76a77//nsqVqzI0KFDGTFihN3rCxUqxLvvvkvfvn0JDAykZ8+eODk5MXPmTDZv3kzlypV57bXX+Oijj277uYiIiIiIiIjkBHcyrwfl9iIiIiIikrNZzPUevpNLxcbG4ufnR0xMDL6+vnb7rly5QlRUFCVKlMDDw8NBEYpITqbvCREREbkXZJY3ieQEyu1FJKv0nSAiIiL3otuV1+sOZBERERERERERERERERERAVRAFhERERE
  149. "text/plain": [
  150. "<Figure size 2400x2400 with 4 Axes>"
  151. ]
  152. },
  153. "metadata": {},
  154. "output_type": "display_data"
  155. }
  156. ],
  157. "source": [
  158. "import pandas as pd\n",
  159. "from math import ceil\n",
  160. "selection = (df[\"x_ref\"] > -3000) & (df[\"x_ref\"] < 3000) \n",
  161. "data = df.loc[selection, \"x_ref\"]\n",
  162. "#func = sqrtSigmoid\n",
  163. "func = fitFastSigmoid\n",
  164. "#NBINS=[800, 900, 1000, 1100]\n",
  165. "NBINS=[900, 1000]\n",
  166. "#SHIFTS = [48, 50, 53, 52]\n",
  167. "SHIFTS=[126., 140.]\n",
  168. "minBinWidths = []\n",
  169. "maxBinWidths = []\n",
  170. "nBinsReal = []\n",
  171. "fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(24, 24))\n",
  172. "for i, BINS in enumerate(NBINS):\n",
  173. " out, equal_bins = pd.qcut(data, q=BINS, retbins=True)\n",
  174. " binWidths = np.array([ abs(equal_bins[j] - equal_bins[j+1]) for j in range(len(equal_bins)-1)])\n",
  175. " minBinWidths.append( binWidths.min() )\n",
  176. " maxBinWidths.append( binWidths.max() )\n",
  177. " bin_numbering = np.arange(0, BINS+1)\n",
  178. " equalbins_center = equal_bins[int(BINS/10):int(9*BINS/10)]\n",
  179. " bin_numbering_center = bin_numbering[int(BINS/10):int(9*BINS/10)]\n",
  180. " print(equalbins_center[0], equalbins_center[-1])\n",
  181. " popt, pcov = curve_fit( func, xdata=equalbins_center, ydata=bin_numbering_center )\n",
  182. " if i<2:\n",
  183. " axes[0,i].plot(equal_bins, bin_numbering, '.', label='data')\n",
  184. " axes[0,i].plot(equal_bins, func(equal_bins, *popt), label=f'p0={popt[0]}\\np1={popt[1]}\\np2={popt[2]}\\n')\n",
  185. " axes[0,i].legend()\n",
  186. " else:\n",
  187. " axes[1,i-2].plot(equal_bins, bin_numbering, '.', label='data')\n",
  188. " axes[1,i-2].plot(equal_bins, func(equal_bins, *popt), label=f'p0={popt[0]}\\np1={popt[1]}\\np2={popt[2]}\\n')\n",
  189. " axes[1,i-2].legend()\n",
  190. " shift = SHIFTS[i]\n",
  191. " nBinsReal.append(ceil(func(3000, popt[0]+shift, popt[1], popt[2])))\n",
  192. " print(f'BINS={BINS}')\n",
  193. " print(f'shift={shift}')\n",
  194. " print(popt)\n",
  195. " print(f'at -3000 = {int(func(-3000., popt[0]+shift, popt[1], popt[2]))}')\n",
  196. " print(f'at 0 = {int(func(0., popt[0]+shift, popt[1], popt[2]))}')\n",
  197. " print(f'at 3000 = {int(func(3000., popt[0]+shift, popt[1], popt[2]))}')\n",
  198. " print(f'at -3000 = {func(-3000., popt[0]+shift, popt[1], popt[2])}')\n",
  199. " print(f'at 0 = {func(0., popt[0]+shift, popt[1], popt[2])}')\n",
  200. " print(f'at 3000 = {func(3000., popt[0]+shift, popt[1], popt[2])}')\n",
  201. " print()\n",
  202. "binData = pd.DataFrame({'NBINS':NBINS, 'nBinsReal': nBinsReal,'minBinWidth':minBinWidths, 'maxBinWidth':maxBinWidths})\n",
  203. "#ax = binData.plot(x='NBINS', y=['minBinWidth', 'maxBinWidth'], logy=True)\n",
  204. "binData"
  205. ]
  206. },
  207. {
  208. "cell_type": "code",
  209. "execution_count": null,
  210. "metadata": {},
  211. "outputs": [],
  212. "source": []
  213. }
  214. ],
  215. "metadata": {
  216. "kernelspec": {
  217. "display_name": "Python 3.10.6 (conda)",
  218. "language": "python",
  219. "name": "python3"
  220. },
  221. "language_info": {
  222. "codemirror_mode": {
  223. "name": "ipython",
  224. "version": 3
  225. },
  226. "file_extension": ".py",
  227. "mimetype": "text/x-python",
  228. "name": "python",
  229. "nbconvert_exporter": "python",
  230. "pygments_lexer": "ipython3",
  231. "version": "3.10.6"
  232. },
  233. "orig_nbformat": 4,
  234. "vscode": {
  235. "interpreter": {
  236. "hash": "a2eff8b4da8b8eebf5ee2e5f811f31a557e0a202b4d2f04f849b065340a6eda6"
  237. }
  238. }
  239. },
  240. "nbformat": 4,
  241. "nbformat_minor": 2
  242. }