cetin/tracking-parametrisation-tuner
1
1
Fork 0
Code Issues Pull Requests repo.project_board Releases Wiki Activity
aa2649befc
tracking-parametrisation-tuner/parameterisations/notebooks/HougHistogram_old.ipynb

435 lines
217 KiB
Plaintext
Raw Normal View History

first commit
2023-12-19 13:00:59 +01:00
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import pandas as pd\n",
"from pathlib import Path\n",
"import numpy as np\n",
"\n",
"dataDir = Path('/scratch/guenther/MCParticleData')\n",
"file = 'bsphiphi500.h5'\n",
"zReference = 8520.\n",
"zMagnetParams = [5212.38, 406.609, -1102.35, -498.039]\n",
"xParams = [18.6195, -5.55793]\n",
"xLayers = [0,3,4,7,8,11]\n",
"uvLayers = [1,2,5,6,9,10]\n",
"\n",
"with pd.HDFStore(dataDir/file, mode='r') as eventData:\n",
" tmp = []\n",
" # event loop\n",
" for iEvent in range(0, 500):\n",
" hits = eventData[str(iEvent)]\n",
" tmp.append(hits.copy())\n",
" allHits = pd.concat(tmp, ignore_index=True)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Hough Histogram Studies\n",
"## Occupancy of the Detector\n",
"* data contains all x hits of 500 events in the range [-3000,3000]mm\n",
"* the first plot shows the coarse distribution of these hits along x\n",
"* the occupancy is higher in the central region than at the edges as expected\n",
"\n",
"## Hough Histogram Binning\n",
"* the idea is to adjust the histogram's bin width along x such that every bin contains the same number of hits on average\n",
"* two approaches have been tested, one of which is shown below\n",
"1. in a **first attempt**, I fitted the (right) flank of the occupancy distribution with a 5th order polynomial to describe the occupancy for most of the range\n",
" - to find regions of equal content, i.e. the desired bins, the culumative distribution must be calculated from which the quantiles (i.e. 500 quantiles for 500 bins) can (in theory) be calculated by $x_{Q}=F^{-1}(p)$ where $p$ is the \"probability\" given by $i_{bin}/N_{Bins}$\n",
" - for polynomials of order higher than 3 this is, however, not analytically doable\n",
" - other functions might be possible, but since there is an easier approach, I didn't look further into this\n",
"2. the **second approach** numerically calculates the quantiles for the full occupancy distribution using `pandas.qcut` function\n",
" - this function directly gives the bin edges containing an equal amount of x hits\n",
" - the bin edges are then enumerated and this enumeration is plotted against the edges, which is shown for different number of bins below\n",
" - the resulting distribution is fitted with $f(x)=p_0+\\frac{p_1\\cdot x}{1+|p_2\\cdot x|}$, which is chosen because it's computionally faster than e.g. a Sigmoid\n",
" - the fitted function gives a direct mapping $x\\to i_{Bin}$ when downcasted to integer values\n",
" - the resulting minimal and maximal width of the bins is shown in the table and the last plot below for different total numbers of bins\n",
" - 1000 bins have minimal bin width of 1.625mm which is of the order of the Hough Projection's resolution for high momentum tracks\n",
" - the parameters for this scenario are $p_0=505.675291$, $p_1=0.453316164$ and $p_2=5.57027462\\times 10^{-4}$\n",
" - to avoid negative bin numbers arising from the imperfection of the parametrisation, the parameter $p_0$ can be shifted by +2.5 for 1000 bins\n",
" - also because of the imperfection of the parametrisation, the real number of bins later on will be slightly higher, which is also shown in the table\n",
" - the inverse. i.e. the mapping $i_{bin}\\to x$ is $x=\\frac{p_0-i_{Bin}}{p_0\\cdot p_2-p_1-p_2\\cdot i_{Bin}}$ for $i<p_0$ and $x=\\frac{p_0+i_{Bin}}{p_0\\cdot p_2+p_1-p_2\\cdot i_{Bin}}$ else"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"BINS=600\n",
"shift=1.1\n",
"[ 3.03414161e+02 2.71945159e-01 -5.56801117e-04]\n",
"at -3000 = 0\n",
"at 3000 = 610\n",
"at -3000 = -0.9960452190164233\n",
"at 3000 = 610.0243669660781\n",
"\n",
"BINS=700\n",
"shift=1.5\n",
"[ 3.53980141e+02 3.17291325e-01 -5.56893995e-04]\n",
"at -3000 = 0\n",
"at 3000 = 711\n",
"at -3000 = -0.9359283090630584\n",
"at 3000 = 711.8962100461829\n",
"\n",
"BINS=800\n",
"shift=1.8\n",
"[ 4.04541076e+02 3.62631197e-01 -5.56941658e-04]\n",
"at -3000 = 0\n",
"at 3000 = 813\n",
"at -3000 = -0.9838529169481944\n",
"at 3000 = 813.6660040895719\n",
"\n",
"BINS=900\n",
"shift=2.2\n",
"[4.55109246e+02 4.07977502e-01 5.56998216e-04]\n",
"at -3000 = 0\n",
"at 3000 = 915\n",
"at -3000 = -0.9217380096175134\n",
"at 3000 = 915.5402300772885\n",
"\n",
"BINS=1000\n",
"shift=5.5\n",
"[5.05675291e+02 4.53316164e-01 5.57027462e-04]\n",
"at -3000 = 2\n",
"at 3000 = 1020\n",
"at -3000 = 2.037684447748063\n",
"at 3000 = 1020.3128978501425\n",
"\n",
"BINS=1100\n",
"shift=4.9\n",
"[5.56237804e+02 4.98644231e-01 5.57020439e-04]\n",
"at -3000 = 1\n",
"at 3000 = 1121\n",
"at -3000 = 1.0860058517370135\n",
"at 3000 = 1121.1896027202529\n",
"\n",
"BINS=1200\n",
"shift=3.2\n",
"[6.06805431e+02 5.43977948e-01 5.57024839e-04]\n",
"at -3000 = 0\n",
"at 3000 = 1220\n",
"at -3000 = -0.95986992994915\n",
"at 3000 = 1220.9707309946493\n",
"\n",
"BINS=1300\n",
"shift=3.6\n",
"[ 6.57370222e+02 5.89309937e-01 -5.57035442e-04]\n",
"at -3000 = 0\n",
"at 3000 = 1322\n",
"at -3000 = -0.9015253398330287\n",
"at 3000 = 1322.841969580662\n",
"\n",
"BINS=1400\n",
"shift=3.9\n",
"[7.07935490e+02 6.34651601e-01 5.57051233e-04]\n",
"at -3000 = 0\n",
"at 3000 = 1424\n",
"at -3000 = -0.9482046457474098\n",
"at 3000 = 1424.6191845963053\n",
"\n",
"BINS=1500\n",
"shift=4.2\n",
"[ 7.58500798e+02 6.80003618e-01 -5.57081680e-04]\n",
"at -3000 = 0\n",
"at 3000 = 1526\n",
"at -3000 = -0.992094855437017\n",
"at 3000 = 1526.3936918183863\n",
"\n"
]
},
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>NBINS</th>\n",
" <th>nBinsReal</th>\n",
" <th>minBinWidth</th>\n",
" <th>maxBinWidth</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>600</td>\n",
" <td>610</td>\n",
" <td>2.750</td>\n",
" <td>72.399902</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>700</td>\n",
" <td>712</td>\n",
" <td>2.375</td>\n",
" <td>61.944386</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>800</td>\n",
" <td>813</td>\n",
" <td>2.000</td>\n",
" <td>54.544922</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>900</td>\n",
" <td>916</td>\n",
" <td>1.750</td>\n",
" <td>49.346866</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>1000</td>\n",
" <td>1017</td>\n",
" <td>1.625</td>\n",
" <td>44.815297</td>\n",
" </tr>\n",
" <tr>\n",
" <th>5</th>\n",
" <td>1100</td>\n",
" <td>1119</td>\n",
" <td>1.375</td>\n",
" <td>40.286656</td>\n",
" </tr>\n",
" <tr>\n",
" <th>6</th>\n",
" <td>1200</td>\n",
" <td>1221</td>\n",
" <td>1.375</td>\n",
" <td>36.824951</td>\n",
" </tr>\n",
" <tr>\n",
" <th>7</th>\n",
" <td>1300</td>\n",
" <td>1323</td>\n",
" <td>1.250</td>\n",
" <td>34.324951</td>\n",
" </tr>\n",
" <tr>\n",
" <th>8</th>\n",
" <td>1400</td>\n",
" <td>1424</td>\n",
" <td>1.125</td>\n",
" <td>31.924148</td>\n",
" </tr>\n",
" <tr>\n",
" <th>9</th>\n",
" <td>1500</td>\n",
" <td>1527</td>\n",
" <td>1.000</td>\n",
" <td>30.625244</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" NBINS nBinsReal minBinWidth maxBinWidth\n",
"0 600 610 2.750 72.399902\n",
"1 700 712 2.375 61.944386\n",
"2 800 813 2.000 54.544922\n",
"3 900 916 1.750 49.346866\n",
"4 1000 1017 1.625 44.815297\n",
"5 1100 1119 1.375 40.286656\n",
"6 1200 1221 1.375 36.824951\n",
"7 1300 1323 1.250 34.324951\n",
"8 1400 1424 1.125 31.924148\n",
"9 1500 1527 1.000 30.625244"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAxsAAAJ9CAYAAABdBG86AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABc8UlEQVR4nO3deZxkVX3//9dbVgVZFNSJOGwRotH8SILBYBREQ1xQUfErMSpgxLjgbsSIChiJosQF1y9GREMSQL4RI1GDCogBl6AxJGIAkUFEkFV2RpbP749zyymKqp7unr7TPT2v5+NRj9t17zn3nnu6uro+dbZUFZIkSZI01+4z3wWQJEmStDgZbEiSJEnqhcGGJEmSpF4YbEiSJEnqhcGGJEmSpF4YbEiSJEnqhcGGpLVakv2TVJIz57ssK5PkuK6sh83BuZZ159p9lQs2u+sf1l3/uPm4/nQl2b0r57L5LsuoJGd2Zdt/ZP+CLTMs/PJNJclvJbkzyTfmuyyaviR/1r3m3jffZVkbGWxI0zD0Ie/MOT7v7t2Hrr3n8rxriiQ7dfe//3yXZbaS7Jjk6CT/neSmJMuTXJbku0k+nuT5SR6wGsoxCB5W9thpFa8zCBKWzWXakXyv6/JuM8tiDv/NDj/uSHJtkh8nOSXJW5NsO9trzKJMm3X3ddjquubq1gXvh63q62wB+xtgHeCvh3cmeerQ6+xJU50gyclduh8mWX+6F+7r/9Ba4gTgx8BBSbaa78KsbQw2pPm1O3AosPf8FmPe7ES7//3ntxizk+RlwHnAq4FHARsBtwBbAo8BXk77J/fiObrkFcAFwDVTpLkF+MUUjzu6dBd357p1jso2U9d0179izLHX0V4X28zBde5gxb1fB9wP2B54FnAEcHGSzyXZckL+W7tyXjwHZdmMdl+HzsG5AH5KK9sNc3S+ubA/7f52miLNXNbpapNkF+DZwHeq6mvDx6rqy8Dx3dNPJrnfhHM8F3gucDfw51X1qx6LrE5V3QW8B9iQufv70zQZbEjSLCR5HPAJYH3ga8BuwIZV9QDgvsAOwEHAt4Cai2tW1V9V1W9V1UemSHZUVT1kiscPu3M9qTvXd+eibDNVVR/prv9XPV/qnKF7f3BV3RfYHHgqcCLtd7MP8J9JHjqmnN/tyjnlt9Xzoape3JXt8/NdlplYyHW6Em/otsdMOP464CpgW1ogew9dC+dHu6dHV9W357qAmtI/0b6MeXGSLea7MGsTgw1Jmp1XA6G1bDylqs4afEtZzUVV9dGq2pXJH040D6rql1X1laraF3g6cDvwUOD/zW/JtFAleSCtNexXwD+PS1NV19LeFwBek+SxI0k+CDwY+AlwSD8l1SRVdSvwL7QviF44z8VZqxhsSKtoeJBmkvt2/ZUvSHJbkquSnJDk4SN5tklSrGjO3W9M//JtxlzrGUm+kOTKJL/qzv/FJH+ykjI+MsmJXfrbkvxvksOTbJgJA3UHZezKSZLHdn2Nr0hyV5IPDqV9VJK3J/lmkp924xau7ermpUnWGVOmAj7dPd1tzP3vPibPH3X1+bOha3wtyZ8myRT3/xtJjklyeZLbk/wkyfuTbDZVva3Eo7vtl7sm+omq6rYpyrZLks+mjbm4Pck1Sb6f5N1JdhxJu6gHiA/2AVt3u84YeU2cOdflqKqvAG/qnu6S5Bkj5Zw4mDnJfbq/+zO61+IdSa7u+uIfm+QpQ2nPBC4Zej76ej9sOO3Qe8pmSY7s/mZvTfLLcemmusfufeOMJNcnuTnJt5K8YELae/zdT0hzrzrpylq0Fj6AT4/c37Kp8o+5xhOT/HNWvNddmeTzSfaYIs+v3zuTLE3yyaH3ikuSHJVkk4kVNbU/AzYAvlpVv5yUqKpOAk6hfb76VLoxGUmeBryoS3Zg98G3d0nW6eryQ0m+l+QXXX3+fFJ9JtkoyY1dXe41xbnT1WuldSkdPb5x2rio/0hyQ/f+dlHaGLeHTTjndF/76yd5bZJzkvyy+9v7RZL/SvLRJH84odgnddsDpqw4zal157sA0iKyCXA28LvAclqf3C2B5wN/nOQPqmrQR/kuWh/yjWn9/G/n3v2uf/0BNsl6tA/mfzZ0/Mbu/HsBeyV5X1W9ebRQSZ4MfJHWV3WQb1vgHcCewJkru7Ek/wf4B9p7xg3DZeucCTxwqNw3Aw+gffDYDXh2kmdV1Z1DeX5B6260Ca1f/XUj57xHX+YkRwLD93cTrQ/8k7rHM5P8WVXdPZLvEcA3aHUFrRn9IcDrgWcAH5/y5lfuXl1vpiNJaH2Ih+/pRuD+tNfQ7wJLWEPHs8zSzbTXxZa0D2vXc8/XwehrZK58Eng77VvnF9D+Xqbj77v0AzfQXs9bAI/sHl/pjl1HG6cy6L7xi5Fz3Tzm/FsC3wO2o72nzLh/f5LX0r5Rr6589wUeCzw2yR9W1aunyD4Tt9Hu6QHAerTX8nCQffUMyvwuVnzzPyj3g2hj2/ZO8p6VdL/7/4Bju7LcRHstbQO8kfbFxq5Vdcfk7GPt2W3PnkbaV9LG4z0SeFuSo4D/2x37u6o6fYbXXhWPAIavN3gdLWFFfR5SVX8zSFBVtyQ5ATiQ9qH81Ann3oNWr7fSxqb9Wve++2VWfHFwZ3ft36S1/rwwyTOqalJ9TnztJ1kXOI0Vge3gNfJA2uvkd7qfvzXmvIPr/U6Sh1TVlROur7lUVT58+FjJAziO9oZ25phjZ3bHrqd9c/kntNlK7gM8HrisO37SmLyHdceOW8n1P9CluwT4U2Djbv/GwMtob7QF/OlIvi1oH3AK+A7wqG7/erQPSTd15b5XGWj/RKp73AScDGzTHVt38HP3/J+BlwJLgXW7fRvRmqqv6M7xl2Pua/9J9TqS7rVduquAVwCbdfs3BJ4H/Lw7/lcj+dYDftgduxh4Qrf/PrRA4yrgl9Mpw5gyfabL9yvgObN4Tf3lUP1+DNh6qGxbA38BHDLhdXjYmPMtm3RswvUH6XefYbkHr9llq5J2qtf+bMs23b/ZCen/sUv/s5H9u4+7B+AJ3f67aH3179/tD+2D3H608TNj/6ZWUpYzWfF391PgKcB9umO/OSbd/hPKfEv3+vwM8ODu2ObAUUOvvRfMtIyT6mSqMs0g/75DZfswsEW3/4HA0UPHXjgm7+DY9cDXWfF+twHwEtqXOgW8coavpQDXdnn3nGaeP2fF+8OXup8vBzZdXa/pLs8OtG/z96IF0+n2Pwh4Gy0IuBvYZSTfHwyVf4sJ5z6+S/PZkf2b0v5XFfB52hcng/8L2wCf7Y5dSfdePpPXPm3CjcHr+4W0sXLQ/u8uBV7FyP+CkWtc2uXfZ7a/Cx8zfO3OdwF8+FgTHlO9yQ+9Od7K0AeBoePP7Y7fDqw/cuwwVhJsAA+nfaC5HthuQpr/053nf0b2H97t/8Xom/pIvpUFG/8+eMOfRd09vjvHJWOO7b+yf5601oubaK0ffzAhzWO7f5jXDdcxrdtC0b4Z23GKss0m2Pjt7p/dIP8yWuvTK4DfB9aZIu8Dh/L+zSxeh4eNObasO3Yz7Z/4uMeBY9LvPsP7Hrxm75riOoPHzYO6meI893rtz7ZsE+pqWr9X4K+GfpfrDe3ffdw90FqkitaNbrpl+vXf1ErSncmKD3qPmka6/Uf27z50L6fRfcCcUD8XDR+fThkn1clUZZpOftqH+ou6Y/80Ie8gKFzGyHvS0D3/D7DBmLwf7o6fPsPX0sOHzv2QGeT76lC+Ap4x29fzbF7T0zzn27tzfnrMsf/qjr12zLFNaf/z7vV3Cryr23/KuNdel+ZfuzRvmulrn/blTAEfn+U9/0uX/31zVY8+pn44ZkOaOydX1Y/H7B+8sW1Aa0KeqRfTvu0+pap+MiHNP9M+UP92kiVD+5/TbY+pMf2Mq/UvnnTOYX9bI92TpquqvklrPdgmyW/M4hTPpbXg/HtNmDmp2qwuP6F9a/v7Q4f26bb/XFUXTCjbWbMoE9VmdXoyreUEWmvE/rR/hOcC1yb5xIS+yc+jTcF6PSPz9c+BjWjfYI57bDSH17nPFNfp43p9u37o5+msi3Jjt31Qkr7+l365qv5
"text/plain": [
"<Figure size 864x720 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABWoAAAI/CAYAAADuquC4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAEAAElEQVR4nOydeVxV1fr/34tRURQnFEFFRJkP4KyVY5rZpGapaVnZeBu++SvLexvvLa9Ntznt3gbTMq2slNSyUiszzRHNEQeQUQQEmYfDeX5/HM6OIxxERAFd79eLF7DPWnuvtc/en/2sZ6/1PEpE0Gg0Go1Go9FoNBqNRqPRaDQaTcPh1NAN0Gg0Go1Go9FoNBqNRqPRaDSaSx3tqNVoNBqNRqPRaDQajUaj0Wg0mgZGO2o1Go1Go9FoNBqNRqPRaDQajaaB0Y5ajUaj0Wg0Go1Go9FoNBqNRqNpYLSjVqPRaDQajUaj0Wg0Go1Go9FoGhjtqNVoNBqNRqPRaDQajUaj0Wg0mgbGpaEbcCbat28v/v7+Dd0MjUZzGtu3b88UkQ4N3Y6zRWuKRtM40Zqi0WjqE60pGo2mPtGaotFo6pOaNKXRO2r9/f3Ztm1bQzdDo9GchlLqWEO3oS5oTdFoGidaUzQaTX2iNUWj0dQnWlM0Gk19UpOm6NAHGo1Go9FoNBqNRqPRaDQajUbTwGhHrUaj0Wg0Go1Go9FoNBqNRqPRNDDaUavRaDQajUaj0Wg0Go1Go9FoNA1Mo49RWx1lZWUkJydTXFzc0E3RaOqdZs2a4efnh6ura0M35ZJBa4rmUkRrzflDa4rmYkZrx4VHa4rmUkHry4VBa4rmUqGpakqTdNQmJyfj6emJv78/SqmGbo5GU2+ICFlZWSQnJ9O9e/eGbs4lg9YUzaWG1przi9YUzcWK1o6GQWuK5lJA68uFQ2uK5lKgKWtKkwx9UFxcTLt27bSoaC46lFK0a9dOv928wGhN0VxqaK05v2hN0VysaO1oGLSmaC4FtL5cOLSmaC4FmrKmNElHLaBFRXPRoq/thkGfd82lhr7mzy/6/GouVvS13TDo8665FNDX+YVDn2vNpUBTvc6brKO2MfHcc8/x6quvOvx8+fLl7Nu37wK2SKPRNGW0pmg0mvpEa4pGo6lPtKZoNJr6RuuKRvMX2lF7AdCiomkSZB2BuB8auhWaWqA1RdMkyEmC/SsbuhWaWqA1RdMkyEuHPV83dCs0tUBriqZJUJAFu79o6FZoaonWFU2jp/gUxH4GIue8K+2orSNz5swhKCiIK6+8koMHDwLw/vvv069fPyIjI7nxxhspLCzk999/JyYmhlmzZhEVFcWRI0eqLafRNAjlZbBvBSy6Ad7uDTEPgaW8oVt1SaI1RXNRYLHA4Z9gyRR40wTf3AdlRQ3dqksSrSmaiwIRiP8VvpgOr4fC13dbnSuaC47WFM1FgQgkboav74HXQqyacjK+oVt1yaJ1RdPkEYHk7bDiAfhPMCy/H1J3nPNuz8lRq5TyUkotU0odUErtV0oNUkq1VUr9qJQ6VPG7TaXyf1dKHVZKHVRKXXXOrT8Lth/L5t31h9l+LPvc97V9O0uXLmXnzp18/fXXbN26FYAJEyawdetWdu3aRUhICB9++CGDBw/m+uuv55VXXiE2NpYePXpUW06juaDkpsK6F+D1MPjiNuts2hFPwb2/gJNzgzRJKRWklIqt9JOrlHpEa4rWFE0ToPAk/P629YXPpzdC0ha4fCb87Xdwbd5gzVJKzVRK7VVK7VFKLVFKNdOaojVF0wQoyoZN8+CdfrDwOoj/BQbcBw9sgRbtGrp1TQKtKRpNJYpPwZb3Yf5g+OgqOPgd9JkO92+Ctuc/G7xS6iOl1Aml1J5K287aHlFK9VFK/Vnx2VvqAgbgrE9NAa0rmiZOaSFs/xj+ewV8MAL2fAMRN8E9P4Nvn3Pevcs51n8T+F5EJiql3AAP4B/AWhF5USk1G5gNPKGUCgUmA2FAZ+AnpVQvETnv0/e2H8tm6gebKTVbcHNxYvFdA+nTrc2ZKzpgw4YNjB8/Hg8PDwCuv/56APbs2cNTTz1FTk4O+fn5XHVV9WO82pbTaOqd5O2weR7sW26dOdvrKuh7JwRe2WAOWhsichCIAlBKOQMpwDdYNURritYUTWMkZTts/RD2fAXmYug6yPrSJ+Q6cHFv0KYppXyBh4FQESlSSn2BVTNC0ZqiNUXTODn+J2x+r0JTisCvH4x7D8LGNehLn6aG1hSNpoLMQ7B5PuxaAmWF4BMJ170F4TeCe8sL2ZKPgXeARZW21WWMMx+4B9gMrAbGAN+d78bXt6aA1hVNEyUnEbZ+ANsXQnEOdIyAa16zOmmbtaq3w9TZUauUagUMAW4HEJFSoFQpdQMwrKLYQuBn4AngBmCpiJQA8Uqpw0B/YFNd21BbNh/NotRswSJQZraw+WjWOQtLdS+vbr/9dpYvX05kZCQff/wxP//8c7V1a1tOo6kXystgf4zVSEneCu6toP+90P9uuzfI249ls/loFgMD2p3z/VEPjASOiMgxrSlaUzSNDEu5dSbK729D0mZwawlRU6HfDOgYZhRrJJriAjRXSpVhfZmcCvwdrSlaUzSNB4sFDv8Im96xhjlw9YDISdB3BviYGrp1TRKtKZpLGlvIlE3vwqE14OxmdaL0uwt8ezdQk+RXpZT/aZvPaoyjlEoAWonIJgCl1CJgHBfAUXs+NAW0rmiaCCKQ8Bv88R4cXA0o66SUAfdB14FwHia2n8uM2gAgA1iglIoEtgP/B3QUkTQAEUlTSnlXlPfF+ubHRnLFtvPOwIB2uLk4UWa24OrixMCAc1syNWTIEG6//XZmz56N2Wzm22+/5d577yUvLw8fHx/KyspYvHgxvr7W7nl6epKXl2fUd1ROo6lXSvJg20fWmSl5qdCmO1z9MkTdAu6eVifKrsO08XBjT+oplm1Pxlxef29Jz5HJwJKKv7WmaE3RNAbKiqwB8je9CyePgFdXGPMSO9uN5ffkUtokuJG9r/FoioikKKVeBRKBIuAHEflBKaU1RWuKpjFQWgi7l1pDHGQdAs/OcOU/oc90tp+AzQeyaJOYyJ7UUyggrHNrsgtLG/rlT5NAa4rmksRcAn8us64eTN8DHu1h6Gzri+SW3mw/ls1X3/xp6IlNWyb09msoTTlbe6Ss4u/Tt5936ltTQOuKpglQboa9X8PGN62a0rwtXPaIVVNa+1n9KT8fMcY+mXkldPB0rxdNORdHrQvQG3hIRP5QSr2Jdbq+I6pzM1ebDk0pdQ/WKf107dr1HJpopU+3Niy+a2C9zezp3bs3kyZNIioqim7dunHFFVcA8PzzzzNgwAC6detGRESEISSTJ0/m7rvv5q233mLZsmUOy2k09UJBpvVtz5b/QfEpkr36crjPP9jbYgBtaEb27+m08cjmXyv3UlJmQbDenLabsT7fktaFijAq12Od9VZj0Wq2aU3RmqKpbwoyrXHdtr4PhVmke4YRF/Uqf7a6gtyTwgcxeyi3iJ2WNAZNqYj1dgPQHcgBvlRKTaupSjXbtKZoTdHUN4UnrTbKH/+FopMkNw9iY9dnsYTcwMliIW/9cT74Ld7Qlco4KRrLC+VGjdYUzSVFca51KfLm+VBwgiKvXqz2m83uNqMI8vBmz4/pZOYlsu7gCczlVR/rX25PZsndjUpTHNkjF42dAlpXNI2Y0kLY+al19eCpRIq8erLabzbrXIbgldeKsAMW1h/cxroDJ7BUY6vUh6YokWrv7TNXVKoTsFlE/Cv+vwKrozYQGFbxVsgH+FlEgpRSfwcQkbkV5dcAz9mm7juib9++sm3bNrtt+/fvJyQkpE7t1miaAnW5xnfv3Uv++tfom/UtblLC3tZDeCbrKnaYu9s5TpwUOCmFRQTLabe/AtxdazcAUkptF5G+Z9XIWlAR6uABERld8f9BtKZoNOeFmq793fsPIL+9QVja17hYStjlMYiXckezydwLQdk5Yx3RkJqilLoJGCMiMyr+vw0YiDW0itYUjeYcqMs1vmv/Qco3vkNE2jJcywvZ5TGIubmj2Wzuhc3/UBt
"text/plain": [
"<Figure size 1728x720 with 10 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEGCAYAAACevtWaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAlbElEQVR4nO3deXQc1YHv8e9tba21tVmyLdmWsYyxwdgYEcBmCAyZifOIHSbhHBKGsGQhhCFgZrIwk8yDvDPJYUggMIEZcAIYXtgCBAIzDMkw4OcxAowMJt7xJtnyot1arf2+P6q61S21bMmSaXXp9zmnTndXVVffKlu/un3r9i1jrUVERLzFF+sCiIjI+FO4i4h4kMJdRMSDFO4iIh6kcBcR8aDEWBcAID8/35aUlMS6GCIicWXjxo311top0ZZNiHAvKSmhoqIi1sUQEYkrxpiq4ZapWUZExIMU7iIiHqRwFxHxoAnR5i4iE0NPTw/V1dV0dnbGuigSxu/3U1xcTFJS0ojfE9NwN8asAFaUlpbGshgi4qquriYzM5OSkhKMMbEujgDWWhoaGqiurmb27Nkjfl9Mm2Wsta9aa28MBAKxLIaIuDo7O8nLy1OwTyDGGPLy8kb9bUpt7iISQcE+8ZzMv0l8t7lXrocD70FRGUw/B/xZsS6RiMiEEN81933/A//9f+DJlXD3THjofHj5Znj/UTi0Cfp6Yl1CETkFXnnlFe6+++7jrlNZWUlqaiqLFy9m0aJFLF26lJ07dwJQUVHBrbfeetz333777dx///2h15/97Gf5xje+EXr9d3/3d9x3333HLUtGRkaoLE8//XRo/po1a7jllluO+/ljFd8190v/Hs7/Fhz6EA5uhOoK+PgPsOkpZ3miH6Ytcmr2RUuguAyyZ4G+dorEtZUrV7Jy5coTrjdnzhw2bdoEwCOPPMJPf/pTnnjiCcrKyigrKzvue5cuXcrzzz/PqlWr6O/vp76+npaWltDy8vJy7r//fs4///wTliUY7ldfffWJd26cxH9vmbRcKL3MmQCshaNVbthvdB4rHoV3H3LXz4eic52p+FyYvsTZhohMCJWVlSxfvpyLLrqId999l0WLFnHDDTdw5513Ultby1NPPcW2bduoqKjgwQcf5PrrrycrK4uKigqOHDnCPffcw5VXXjlkuy0tLeTk5ACwdu1afv7zn/Pv//7v3HXXXezfv5+9e/eyf/9+Vq1axa233sqyZcu4/fbbAdi6dStnnXUWhw8fpqmpibS0NLZv384555zDmjVrQmXZt28fV199Nb29vSxfvjz02XfccQfbt29n8eLFXHfddeTk5HDo0CGWL1/Onj17+Ku/+ivuueeecT2OMQ13a+2rwKtlZWXfHLeNGgM5Jc501peceX09ULvNqdkf/AAOVsCuPwLuLQZz57hhX+Y8Tl0IiSnjViSRePTjV7ey7VDLiVcchQXTs7hzxZknXG/37t08//zzrF69mvPOO4+nn36a9evX88orr/DTn/6UK664ImL9w4cPs379enbs2MHKlStD4b5nzx4WL15Ma2srHR0dvPfee1E/b8eOHbz11lu0trYyb948vv3tbzN9+nQSExPZv38/5eXlXHjhhRw8eJB33nmHQCDA2WefTXJycsR2brvtNr797W9z7bXX8tBDD4Xm33333aGTCTjNMps2beLDDz8kJSWFefPm8Z3vfIcZM2aM5nAeV3w3y4xUQpLTPDNtEZz3dWdeZ8tAc87BjbBvHWz+rbPMl+QEfDDsi8og9zTwxfclCpF4MXv2bBYuXAjAmWeeyWWXXYYxhoULF1JZWTlk/SuuuAKfz8eCBQuoqakJzQ9vlnnuuee48cYbef3114e8//LLLyclJYWUlBQKCgqoqamhuLiYZcuWUV5eTnl5OX/7t3/LwYMHKS8vJxAIsHTp0iHbefvtt3nxxRcB+OpXv8oPfvCDYffxsssuI9gNfMGCBVRVVSncx4U/C077tDMFNR90w96t4X/4FGxY7a4fGGjOKXJDPyPqSJsinjCSGvapkpIy8M3Z5/OFXvt8Pnp7e4+7vrU26jZXrlzJDTfccMLPS0hICH3G0qVLKS8vZ/PmzZx11lnMmDGDe++9l6ysLL72ta9F3dZIuy0O95njZfKGezSBImda4F4c6e+Dup1u2Ltt+P9zL9h+Z3n2zIGwn74YChao/V5kglq/fj1z5swZ1XuWLVvGvffey2mnnUZCQgK5ubkcPXqUrVu38qtf/Srq+s8++yzXXHMNTz31VGh+ZmYmra2tY96H0VC4H48vAQoXONOSa5153e1w+KOB3jnVG2HrSwPvyZjqrF8QnObDlDMgOS02+yAyiQXb3K21JCcn8+tf/3pU71+4cCH19fURvVwWLlxIW1sb+fn5Q9Z/4IEHuPrqq3nggQf40pe+FJp/9tlnk5iYyKJFi7j++utDF3ZPJTPcV5hPUllZmY3rm3W01kDNZqjZBrXboXarU+PvDf5c2EDu7MjALzzTuZCboPOrTBzbt29n/vz5sS6GRBHt38YYs9FaG7VPZ/x3hZwIMgudqfQzA/P6+6Bxn9NLJzjVbIOdrw006yQkQ/48N+zDgj8wQ33xRWRMvNcVcqLwJUB+qTMtCPuBQ08n1O90avg1W53HqrcHeuoAJGcOCnx3Ss/75PdDROKS2gQ+aUn+gW6Z4Y4dhbodA4Ffuw22vgwb1wysk1HohH544BecAcnpn+AOiEg8ULhPFKnZMPMCZwqyFlqPhDXtuLX9iseh95i7koGcWVBwphP8+ac77fs5syE9X807IpOUwn0iMwaypjlTcHgFcNrzmyojA792O3z8Oti+gfWSM5yQzy1xH2cPPGYV62KuiIfprzse+RIgb44zzV8xML+3ywn9xn3QtG/gsXaHM6BaX3fYNhKdfvqDQz9ntjN0g7puisQ1hbuXJKbAlHnONFh/H7Qcigz94GN1BXQ1R66fMXVo6Acf03LV3CNxZc2aNXzve9+jqKiInp4e5s+fz5NPPklaWhoPP/wwaWlpXHvttcO+/5xzzuHxxx9n8eLF9Pb2EggEeOSRR7jmmmsAOPfcc/nVr37Fyy+/zMUXX8xnPvOZiPeHD1S2du1akpOTQ8MXXH/99Xz+85+POtjZWCjcJwtfAmTPcKbZF0cusxaONQ0N/cZ9sPct+OjpyPVTspzafbTwzypyPktkgrnqqqt48MEHAbj66qt57rnnuOGGG7jppptO+N7gMASLFy/mo48+Yt68eZSXl3PNNdfQ3t7O3r17WbRoEUuWLDnhttauXUtGRkbUsWnGk/q5i1MLT8t1puJzhy7vORa9uefIFtjxGvSH3RTFl+Q09+S6zTvZs9xROt1Hv+6XK8c3kiF/AVatWsWxY8dITU3l8ccfZ968edx3331s2bKFxx57jM2bN/OVr3yFDRs2RGy/t7eX9vb20K9E77rrLjIyMvjud7/LJZdcwvnnn89bb73F0aNHefTRR/mzP/szli1bxmuvvcbNN99MeXk5N910E2vWrAFgw4YNLFmyhISEhIha+Ouvv86qVavIz88PhX5lZSUPP/wwCQkJ/OY3v+GXv/wlAOvWreO+++477pDFo6V+7nJiSaluF8wov1zs74Pmaifsw08ATZVQ/T50DmruSc0JC/ySgdDPKXF+vJWQdKr3RkbqP++AI5vHd5tTF8Lnjn8HJTjxkL9PPvkk69atIzExkTfeeIN/+Id/4MUXX2TVqlVccsklvPTSS/zkJz/hkUceIS3NuX703HPPsX79eg4fPszpp5/OihUron52b28vGzZs4LXXXuPHP/4xb7zxBkuXLuVHP/oR4Nyk48477+SZZ56htbWV8vJyli1bFrGNzs5OvvnNb/Lmm29SWlrKVVddBUBJSQk33XRT6GQC8Oijjw47ZPFYqFlGxsaX4Ab0rOjLjzVBU5UT9kfdx6ZKJzR2/Edkrd/4nF48we3llEB2yUD4q2vnpHGiIX+bm5u57rrr2LVrF8YYenqc/0c+n481a9Zw9tln861vfSsidIPNMtZa/uZv/oaf/exn3HHHHUM++4tf/CLgtKMHhxcuKSm
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"from scipy.optimize import curve_fit\n",
"import matplotlib.pyplot as plt\n",
"\n",
"\n",
"data = allHits.loc[(allHits['pc'].isin(xLayers)) & (allHits['x']<=3200) & (allHits['x']>=-3200)]['x']\n",
"\n",
"urax = data.plot.hist(bins=25, title='Integrated SciFi Hit Distribution (X Layers)', figsize=(12,10), fontsize=20, edgecolor='black')\n",
"urax.set_xlabel('x [mm]', fontsize=24)\n",
"urax.set_ylabel('Count', fontsize=24)\n",
"urax.title.set_size(24)\n",
"#out_0,equal_bins_0 = pd.qcut(data, q=25, retbins=True)\n",
"#urax0 = data.plot.hist(bins=equal_bins_0, title='Integrated SciFi Hit Distribution (X Layers)', figsize=(12,10), fontsize=20, edgecolor='black')\n",
"#urax0.set_xlabel('x [mm]', fontsize=24)\n",
"#urax0.set_ylabel('Count', fontsize=24)\n",
"#urax0.title.set_size(24)\n",
"def fit3(x, p0, p1, p2, p3, p4, p5, p6, p7):\n",
" return p0 + p1 * x + p2 * x**2 + p3 * x**3 + p4 * x**4 + p5 * x**5 + p6 * x**6 + p7 * x**7\n",
"\n",
"def fitFastSigmoid(x, p0, p1, p2 ):\n",
" return p0 + p1 * x / (1 + abs(x * p2) ) \n",
"\n",
"def FastSigmoidInverse(y, p0, p1, p2):\n",
" return (y - p0) / (p0 * p2 + p1 - p2 * y)#(p0 - y) / (p0 * p2 - p1 - p2 * y)\n",
"\n",
"NBINS=[600, 700, 800, 900, 1000, 1100, 1200, 1300, 1400, 1500]\n",
"SHIFTS = [1.1, 1.5, 1.8, 2.2, 5.5, 4.9, 3.2, 3.6, 3.9, 4.2]\n",
"minBinWidths = []\n",
"maxBinWidths = []\n",
"nBinsReal = []\n",
"fig, axes = plt.subplots(nrows=2, ncols=5, figsize=(24, 10))\n",
"for i, BINS in enumerate(NBINS):\n",
" out,equal_bins = pd.qcut(data, q=BINS, retbins=True)\n",
" binWidths = np.array([ abs(equal_bins[i] - equal_bins[i+1]) for i in range(len(equal_bins)-1)])\n",
" minBinWidths.append( binWidths.min() )\n",
" maxBinWidths.append( binWidths.max() )\n",
" #equal_bincenters = np.array([0.5 * ( equal_bins[i] + equal_bins[i+1] ) for i in range(len(equal_bins)-1)] ) \n",
" popt, pcov = curve_fit( fitFastSigmoid, xdata=equal_bins[:-1], ydata=np.arange(0, BINS) )\n",
" if i<5:\n",
" axes[0,i].plot(equal_bins[:-1], np.arange(0, BINS), '.', label='data')\n",
" axes[0,i].plot(equal_bins[:-1], fitFastSigmoid(equal_bins[:-1], *popt), label=f'p0={popt[0]}\\np1={popt[1]}\\np2={popt[2]}\\n')\n",
" axes[0,i].legend()\n",
" else:\n",
" axes[1,i-5].plot(equal_bins[:-1], np.arange(0, BINS), '.', label='data')\n",
" axes[1,i-5].plot(equal_bins[:-1], fitFastSigmoid(equal_bins[:-1], *popt), label=f'p0={popt[0]}\\np1={popt[1]}\\np2={popt[2]}\\n')\n",
" axes[1,i-5].legend()\n",
" nBinsReal.append(abs(int(fitFastSigmoid(-3000, popt[0], popt[1], popt[2])) - int(fitFastSigmoid(3000, popt[0], popt[1], popt[2]))))\n",
" ###############33\n",
" # SHIFT A BIT UPWARDS TO AVOID NEGATIVE INDICES p0+2.5\n",
" ################################3#\n",
" shift = SHIFTS[i]\n",
" print(f'BINS={BINS}')\n",
" print(f'shift={shift}')\n",
" print(popt)\n",
" print(f'at -3000 = {int(fitFastSigmoid(-3000., popt[0]+shift, popt[1], popt[2]))}')\n",
" print(f'at 3000 = {int(fitFastSigmoid(3000., popt[0]+shift, popt[1], popt[2]))}')\n",
" print(f'at -3000 = {fitFastSigmoid(-3000., popt[0]+shift, popt[1], popt[2])}')\n",
" print(f'at 3000 = {fitFastSigmoid(3000., popt[0]+shift, popt[1], popt[2])}')\n",
" #print(f'inv at 0 = {FastSigmoidInverse(506, popt[0]+shift, popt[1], popt[2])}')\n",
" #print(f'inv at 1016 = {FastSigmoidInverse(1016, popt[0]+shift, popt[1], popt[2])}')\n",
" print()\n",
"binData = pd.DataFrame({'NBINS':NBINS, 'nBinsReal': nBinsReal,'minBinWidth':minBinWidths, 'maxBinWidth':maxBinWidths})\n",
"ax = binData.plot(x='NBINS', y=['minBinWidth', 'maxBinWidth'], logy=True)\n",
"binData\n",
"#out,equal_bins = pd.qcut(data, q=1000, retbins=True)\n",
"#fig, axis = plt.subplots(figsize=(12,10))\n",
"#axis.set_title('Mapping Bin Number to x Position', fontsize=24)\n",
"#axis.plot(equal_bins[:-1], np.arange(0,1000), '.' )\n",
"#axis.set_xlabel('x [mm]', fontsize=24)\n",
"#axis.set_ylabel('Bin Number', fontsize=24)\n",
"#axis.tick_params(axis='both', which='major', labelsize=20)\n",
"#axis.tick_params(axis='both', which='minor', labelsize=16)\n",
"#popt, pcov = curve_fit( fitFastSigmoid, xdata=equal_bins[:-1], ydata=np.arange(0, 1000) )\n",
"#axis.plot(equal_bins[:-1], fitFastSigmoid(equal_bins[:-1], *popt), label=f'p0={popt[0]}\\np1={popt[1]}\\np2={popt[2]}\\n', color='red', markersize=24)\n",
"#axis.legend(fontsize=20)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3.10.6 (conda)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.6"
},
"vscode": {
"interpreter": {
"hash": "a2eff8b4da8b8eebf5ee2e5f811f31a557e0a202b4d2f04f849b065340a6eda6"
}
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Reference in New Issue Copy Permalink