Projektpraktikum/trackinglosses_rad_length_beginVelo.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "import uproot\t\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pylab\n",
    "from mpl_toolkits import mplot3d\n",
    "import awkward as ak\n",
    "from scipy.optimize import curve_fit\n",
    "from methods.fit_linear_regression_model import fit_linear_regression_model\n",
    "import sklearn\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "40402 10099\n",
      "50501\n"
     ]
    }
   ],
   "source": [
    "file = uproot.open(\n",
    "    \"tracking_losses_ntuple_B_default_radlength_beginVelo.root:PrDebugTrackingLosses.PrDebugTrackingTool/Tuple;1\"\n",
    ")\n",
    "\n",
    "# selektiere nur elektronen von B->K*ee\n",
    "allcolumns = file.arrays()\n",
    "found = allcolumns[\n",
    "    (allcolumns.isElectron) & (~allcolumns.lost) & (allcolumns.fromB)\n",
    "]  # B: 9056\n",
    "lost = allcolumns[\n",
    "    (allcolumns.isElectron) & (allcolumns.lost) & (allcolumns.fromB)\n",
    "]  # B: 1466\n",
    "\n",
    "electrons = allcolumns[(allcolumns.isElectron) & (allcolumns.fromB)]\n",
    "\n",
    "print(ak.num(found, axis=0), ak.num(lost, axis=0))\n",
    "print(ak.num(electrons, axis=0))\n",
    "# ak.count(found, axis=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "stretch factor:  0.24996287312509283\n"
     ]
    }
   ],
   "source": [
    "rad_length_found = ak.to_numpy(found[\"rad_length_frac\"])\n",
    "eta_found = ak.to_numpy(found[\"eta\"])\n",
    "phi_found = ak.to_numpy(found[\"phi\"])\n",
    "\n",
    "rad_length_lost = ak.to_numpy(lost[\"rad_length_frac\"])\n",
    "eta_lost = ak.to_numpy(lost[\"eta\"])\n",
    "phi_lost = ak.to_numpy(lost[\"phi\"])\n",
    "\n",
    "stretch_factor = ak.num(eta_lost, axis=0) / ak.num(eta_found, axis=0)\n",
    "print(\"stretch factor: \", stretch_factor)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjYAAAHLCAYAAADbUtJvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAxt0lEQVR4nO3dTWwb653n+59u9xzFieBTkrPJAAHiItCbXoxR1NkEjdiAyc4iS5M2cIHJaeBAJALtDDQJ9cbypgkSvbkLDVI0MnMmAwxgszwvwO1FN3kadjAIcMdmtbO5mHtxVW7AiwzQORJjaOLo5IV34VQdUXwRVWSxisXvBxDOUZFVz5/10FV/PfW8rPX7/b4AAABS4H+LOwAAAIB5IbEBAACpQWIDAABSg8QGAACkBokNAABIDRIbAACQGiQ2AAAgNUhsAABAapDYYGm5rqtmsxl3GCuv1+up2WzK87y4Q0kVz/PUbDbV6/XiDgVYKiQ2WDqe56lYLCqbzcq27YHtm5ubc012ojjmtFzXVbVaVTabVTabXXj502g2m7p+/brK5fLMN2D/8xaLRWUyGTUajfkEGVKcdd9oNJTJZFQul3V0dLSwcmf5zJ1OR5lMRmtra1pbW1M2m5XjOEPvcxxH2Ww2eI/ruhOPuwz/DpAsJDZYOqZpqtVqDW3v9Xrq9Xo6PDwMdVzP84ZuzrMecxaWZenevXsXXvgX6fw5KpVKKpVKMx/XdV3dvn1b9XpdrVZLhUJhoec8aXVfqVRUKBQWXu4snzmXy6nb7Qa/l8vlkZ+hUCioXC7LMAx1u11ZljXxuEn8d4Bk++O4AwDmxbIszbL0WbFYVKvVkmEYczvmrC666C/aqHN07dq1mY9bq9W0tbUV/F6v12c+5mUkse7Pno9FmfUzG4ahSqWiRqOhVqs1Nuntdrva29u7VFzAtGixAfT+xsZfhJNFeY7iPPfU/Xz5CUun0xl7Xp88eTKXlj5gFBIbRM5xHOXzeXU6HTWbTW1ubqpcLgevNxoNlcvl4Dn6qOf7vV5P5XI5+BnV/8JxHBWLRRWLxaHXJpXhOE5wAS6XywM3uknH9GOqVqvK5/PBZzz7erPZDPoadDqdoG/BqONdluu6KhaLyufzymQyqlarM5Xtn9tMJjN0s590jnxHR0dBmdN+xmazqWKxKM/zgr5TxWIxOI/z/u749XS2fpNe935Zm5ubQ5/fN+m7cNakOh73mS/7WQzDCB5B1Wq1odcdx1EulxtoHZs2/lEuqgusoD4QoVar1TdNsy+pXyqV+pVKpW9ZVt+yrH6/3+9XKpX+2a9hu93uS+q32+1g2+HhYd8wjIFt9Xq9Lyk4zuHhYbAtl8sNxDBNGf57Dg8PB8odd8xut9s3DKPf7XaDbbZt9yX16/V6sH+hUAj2r1Qq/W632y+VSgPvu8jZz3m2/LMxtVqt4ByHKduyrH6lUgn2ldSX1DdNMzjmqHPU7/cHzlG9Xg/1GU3T7JumObBtXt8d0zQHYjYMo28YRvB7Uuvef28ul+uXSqW+bdt9y7KGYrrou+C7qI7HfeYwn+Xs8Y+Pj4fiOFtH08Y/7t/BRXWB1UNig8j5F0v/onpWLpcbuMkcHx8PvbdQKAzdXPr90Re6UTeiacoYd9Med0zLskbG5N94/OP4F+nzF9lRxxxn1Oe0LGvgYt7vv79hn72RTFu2/76zx/NvWmfPx0WJjW3boT/jqMTm7LHDfncsyxr6/P4x/fOU1Lr366DVag1s95M9Pzm4zHfhojoeF1+Yz5LL5Ybqw/8j5axp4vfLGvXvYJq6wGrhURQi5zc5f/TRR0OvtVqtgZEUL1++lKRghIrnecHjiLAuKuOyPM+T67ojOzT6jwnODkOXNNDs7gs7jNcvv1arBY8Nzj4S8D/ftGW/ePFi6HX/eJfpe7K9vT1U7qxz28z63XFdV7lcbmC/SqWifr8/8rxcJI66P1+WX0673Z76uzCvOr7MZ/EfJ519PFiv1wc6DV/2u3xWmLrAamBUFBZm1EXRMAwZhiHHcfT48eOhG5h/YzRNc6ZyJ5VxWZNuBP7NPcrJ6vzyRw15D8M/H51OZ+gmkZTRKGG+O/55CpPAjBN33UsKEjX/xi5d/F2Io45zuZxM0wwmGiyVSnry5Ilev34dvGeW73IS6gLJRIsNYuV5nrLZrDzPU6vVUqVSGXpdCt+6MU0ZYY1q8fFvolEO1fXPybwu2oVCQblcTrVaTZ1OR71eT/V6XZVKZaaEMmrTfneiuLnFVfdny/GTBunizxhXHfutNvV6XY7jaHt7eyDRnEcdxVkXSCYSG8Qqn89ra2trbLLhX3TPPnKYdxmX5f+FO2rkhX+RzWQycylrFP+cjJrVdVxcF2m1WsrlcsEyFfV6feFzyVzWRfXq19O41oAwN9O46/5sOR999NGlvgtx1HGpVAoeSe7s7AyN6Jrlu5yEukAykdggNv4Q37N/wfkXJL+Fxm9SHrdmzkX9ZKYp4zLHk95fjC3LCo591suXL2UYRqRzdPiPIqrV6lBzfNjp//0J6iqViiqVysTHE+fP0eeffy5ptla1o6OjS+1/2e/O+ZtftVod+mt+Gepeen8jN00zaIWRpvsuXKaO5+lsn5rzMxHP8l1OQl0gmUhsEDn/hnH+xuHfWBzHUbPZVLPZDJquXdcN/orz/yLPZrPqdDryPC94n+d5wZw2oxKWacro9XrBX3a2bQcdlscdU1IwS+3Zv0D95v1Hjx4FN9wo1vnxZ3eV3p+TYrGoRqOhfD6vw8PD4GYxbdn+jd8/TrPZHJjfxTfuHI3T6/Vm/vyzfnf8Fol8Pq9isRjMd5PJZII6Smrd+8c5e9P2PE+2bQetUNN+F6at43GfeZbP4icXo5KMaeMfZ9q6wIqJe1gW0u3sXCSmaQ4NCbZtu28YRt80zWBYaKlU6huGMTBM1Lbt4DiWZQXzk1Qqlf7h4WG/2+0Gc23oD8NS/aGi05ZhWVbfMIxg/oxJx+z33w8v9oeil0qlfqlUGhi22u12g2Gnpmn22+12//j4OBhmqxHDZ886O1fIqPLr9frYc3uZsv1z6W8/+2Oa5kCZ589Rq9UKhub6w3bPlzNqqPa4z1gqlYJhzPP67rRareBcnJ9DZdznirvu/TIqlUo/l8sNlHN+Xph+f/J3od+fro7HfeZ5fBZ/npxxLvouT/p3cFFdYPWs9fsxLoYCIHau6+rx48fa29vT0dHRQEtLq9VSJpOZW/8kxIM6xiphuDewwvyRRcfHx8Hw6bNM02R6+iVHHWPV0McGWGF+/42dnZ2B/hb+3CO2bdMBc8lRx1g1PIoCVlyj0VCtVhvooGtZlur1+oWdN7EcqGOsEhIbAJK+HMWU5En5MBvqGKuAxAYAAKQGfWwAAEBqLPWoqF/84hf6u7/7O33rW9/SlStX4g4HAABM4d27d/qnf/onffe739XXv/71uR57qRObv/3bv9Vf/MVfxB0GAAAI4dNPP9XHH38812MudWLzjW98Q5L06NGjSNc9uXPnjp4+fRrZ8dNSxsnJiW7evKnnz59rY2MjkjLScJ4WUcYi6kKK/nNQF9NLw7miLlanDNd1tbOzE9zH5ymyxOYHP/iBtre39cknn0RVhL7yla9Ikv7kT/4k0sTmypUrkS8Yl4Yy3r59K0m6ceOGrl69GkkZaThPiyhjEXUhRf85qIvppeFcURerU8bJyYmkL+/j8xRJ5+HXr1/Ltm2m6AYAAAsVSYvN9evXZds2cyUAAICFiuxR1M7OTlSHBgAAGCmyeWy++93v6tWrV1EdHgAAYEjoFpt79+6Nfa3X66nT6ejJkye6ceNG2CKm9v3vf19f/epXR762u7ur3d3dyGMAAABfOjg40MHBwcjXfvWrX0VWbujEptVqTfWev/7rvw5bxNR+/OMf6zvf+U5kx19EYpSWMqKWlvOUhrqQov8c1MX00nCuqIt0lTGpYeEnP/mJbt68OdPxxwm9VtTdu3dVr9e1tbU19Nrh4aGazaZ++MMfzhzgJP6Jef78eaSJDabz9u1bffjhh/rlL38Z6VBKXIy6SA7qIjmoi+SI8v4duo9NuVz
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(\n",
    "    rad_length_lost,\n",
    "    bins=100,\n",
    "    density=True,\n",
    "    alpha=0.5,\n",
    "    color=\"darkorange\",\n",
    "    histtype=\"bar\",\n",
    "    label=\"lost\",\n",
    "    range=[0, 1],\n",
    ")\n",
    "plt.hist(\n",
    "    rad_length_found,\n",
    "    bins=100,\n",
    "    density=True,\n",
    "    alpha=0.5,\n",
    "    color=\"blue\",\n",
    "    histtype=\"bar\",\n",
    "    label=\"found\",\n",
    "    range=[0, 1],\n",
    ")\n",
    "plt.xlim(0, 1)\n",
    "# plt.yscale(\"log\")\n",
    "plt.title(\"radiation length fraction beginVelo\")\n",
    "plt.xlabel(f\"$x/X_0$\")\n",
    "plt.ylabel(\"a.u.\")\n",
    "\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABiMAAAL7CAYAAACfnvZxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAADuDUlEQVR4nOz9cXQd133Y+/42KZNBq0gHoO/SdfzU0AdsXuh7zeidQ0d5YRU7T0DNu2qtdSkCtkjcrtSpBaSB2MZuCkRNXJpxHBpMWt1bCq0BuXHy7iJhGaC4VuW+0gHcyHXo2jFxqlKvZl4uAdPljVLdSjjHslpYsMD9/oDOEcCzf79zZnBmMAC+n7WwJM5vZu/f3rNnzsyZMzPOe+8FAAAAAAAAAAAgITs2OgEAAAAAAAAAALC1cTECAAAAAAAAAAAkiosRAAAAAAAAAAAgUVyMAAAAAAAAAAAAieJiBAAAAAAAAAAASBQXIwAAAAAAAAAAQKK4GAEAAAAAAAAAABLFxQgAAAAAAAAAAJAoLkYAAAAAAAAAAIBEcTECAAAgJaVSSc6cOSPFYlG6u7s3Op1NY35+XsbHx7dkn83MzMjAwIC0t7fL1NTURqeDTWqr7VsqlYpMTU1Jb2+vtLe3S6lU2uiUtrStvI8FAADZwsUIAACAFJRKJXn66adleHhY/WJtampK2tvb5cyZM6nklHZ9UZVKJenu7pbOzk4ZGBiQmZmZjU6ppaampmRsbEzGx8elUqkkUseZM2dqX+i2t7fLwMBAYnW1WtbHZ1Y0s2/ZbMbHx2VsbEympqYyPV5bPUanpqakWCyKc27NX3d3t1pHqVSS3t7eNfP39vbK/Px8w/q2+j4WAABkj/Pe+41OAgAAYLvo7e2Vqakp6erqkunp6TWxYrEopVJJ8vm8zM3NtaS+6hd5uVyuLpZEfUkYHx+XgYEBERHZrIeu1npwzomIyOTkpPT09LSsvgcffFC6urpkZGREKpWKFItFmZ+fl5GRERkaGmpJPeu1FcZnVlj7llaz1lurzM/PS2dnp4iIzM7OSqFQSKyuuJIao93d3bULA4VCQWZnZxsuc+bMGRkeHo61fW+FfSwAANgcuDMCAAAgRR0dHWpsZGSk9uVxq1i/kE2iviTk8/mNTmHdmv2lcqtUfyVf/YIxl8vJ3NycTE9PZ+ZChMjWGJ9ZYe1bWi2N8Zxme+JKaoxOTk7W/r/Zu11eeeUVyeVysbbvrbCPBQAAm8MdG50AAAAAVnR1dUlXV1fLyhseHjYfu9Hq+hDWaD0kYXx8XETqv9DN0vpmfG5OGzGesyqpMZrL5aS/v7+2HZ85c6bhRYbx8XHp7+9veS4AAACtxJ0RAAAAW9DU1BTP2s+AjVgPad6BERfjc3NivaVn9d0Wp0+fNuedmZmRSqUijz/+eNJpAQAArAsXIwAAwLY2MzMjvb29tV+gDgwM1F4YutrU1JR0d3dLb2+vdHd3S7FYrC2jGR4els7OTikWi03NPz8/L8PDw1IsFtV5msljampKHn300dq/H3zwwdoLjKPWJ7Lyq9xqXZ2dndLd3a22pVKpyPj4+Jp5xsfHpbOzU5xztWest1qlUpGBgYE19YS+NG1FftX+qK7bgYGB4It2m10PtxseHpb29vbI/TU+Pi7t7e1r1ue73vWuWp2r+yPJcb+6f9rb26W7u3tNG7bq+Gy2r1qRQ9R9SyPNbD9RxnOccdPI/Pz8mpexN3pUVLP7hNWa3bZX56SN0fWu51wuV3uHTKVSkampKXXekZER6enpWfMOjzjtbyTKtgYAABDkAQAAtqGRkRGfz+e9iHgR8WNjY76rq8vncrnatLm5Oe+9911dXV5EfLlcri0/NDTkRcSPjIzUlV0ul32hUPBdXV1rlunp6amV3dXVVZs+Oztbq0NEfC6XC+YcJY9yuVwrb3Z2dk2s2fpmZ2d9Pp/3PT09a8qttuP29k1PT/tCoVArd2RkxPf09PhCoeD7+/vX9Pfq5RqZnp6uLafFC4WCn56e9t57Pzc3V2tfoVBoaX6FQmFNmZOTk7VlVv+t7i9tPVStHoPVcTM0NBS7v+bm5urGcFWS435ubs7n83nf399fW2Z1LpOTk031y2Ybn1H6ar05RN23NKPZ7adaf6PxHHXcaFbXNTIy4kXE5/P5NWNVRGp5x21TVZRtu9EYbdVYW7395PN5c57V/RC1/Y32sVG3NQAAAA0XIwAAwLa1+suufD5f++JmZGTE9/f3e+/f+pLm9i+bql8Ahb4oLRQK6heo1S/SQl8Y9vf3q2VGzaOZLw2t+qq5NrowEmpHtY25XG7NF2TlcrkWi/KlpPVFWbX9t7dxdfuHhoZakl/1C9XVX6p7/9YXwfl83pfL5TVfykW5GJHP5+vmqeZzexss1sWI23Nq9bhf/WVlVXWsrP4ydSuOzyh9FTeHuPsWTdTtp9F6i9MXmtvH6eqxvHqfoNUVZZ8QZ9v2vrkxut6xtvqiR+jCS/XC5Xra3+hiRNxtDQAA4HY8pgkAAGxbuVyu9liL1S8iHRoakrGxsTXz5vP5Nf+uvhi4UqmseYTH+Pi4lEol9dnd1stOGz2OJkoezbDqGx4elkqlor4Qtfo885mZmbrHh1Rz7O/vX9PeXC4nBw8eFBGRubm5SLlaeeZyOSkUCmumr552+2NE4uZXbeft62BgYEBE3npXw+pHpUQxMjJS145qfq18D0QS4354eFgd9yMjI5LL5SK/6Hezjc8o22acHNazb9HE2X6a0cr9lMjK+lxdZldXV22sViqVNY8fitOmuNt2o312K8ba8PBw7f9Xv0eianx8vJZndf5WrtP1bGsAAAC3u2OjEwAAAMgC7Uulrq4umZubq/uSavWXwwsLC7UvqapfzNz+RVBV9Uu5qLGoeTTDqq/6pdJ73/veYLxQKEg+n5f5+Xl5+umna882b6Saf6u+XK/m2dnZWRer1lH9ArSZvrHyq05bWFhYM331l4xR18F68mmFVo376peboXFfKBSkXC5Hzm2zjM9Wb5taDuvZt2havf0ksZ8Sqb9IILLyBX/1i/Jvfetbtelx2hR3247T51XNjrWuri4pFApSKpVkZmZGSqVSbQxU32ex+kJBq9dpUtsaAADYnrgYAQAA0MDqL8Kmpqbk6aefDn45JvLWlz1aPK081quZL2OrX0Al8ULqZlTr7erqkunp6VTr1X6FntT62AjNjrf5+XmpVCqptj1r4zONbbPV+5aktp8091MHDx6sfUEvsv42ZXXbfvzxx6W3t1dERE6fPi2Tk5O1/1/94uok1mnWtjUAALC58ZgmAACAJoyPj0tnZ6csLCzI5ORk8FEpSf16PWoe63X7r5g11S/nrHmSVK03jX4Xkdqvj29/lFH1sTNb8RfBUcZ9WuMgq+MzyW0ziTGe5PaTxn5KpP4CQdw2ZX3b7unpqbV1amqq9sX/7Y/tavU6zeq2BgAANi/ujAAAAGigu7tbZmZmgo8fWW31s9Dn5+db/kvaZvNYr9VlW88zr/4ad6N+MVx9REpaFyPGxsZkfn5eZmZmpLe3V0ZGRqSjo0N6e3ulUCjUfq28VTQ73m5/H0ArH1MVksXxmfS2mcS+JantJ639lMhb67j62KK4bdoM2/bw8HDt3RDVR3bl8/k1j+1q9TrN4rYGAAA2N+6MAAAAMAwPD8vMzIz09/c3/KJl9ZdCrX5cRZQ8WqFax8zMjDpP9Quv6otY03b742A0lUol1kt4Q6anp6VQKEhHR4cMDw/XviCcnZ1tSflZEXfcW/08Pz/f8i9JszA+09g2k9i3JLH9pL2fqq7j6vsM1tOmrG/b/f39a97RMj4+vubl1iLJrNMsbWsAAGDz42IEAACAwXoZaEj1i5vTp09vaB6rxXl0RvUXuKVSSf0Cufql6O1fiKUll8vV+vvRRx9d8+vx1Xp7e9XnwEfV29srk5OTMjY2Vvtvs49w2UyPMIk63qpflp8+fVodLwMDA8EvqDf7+FzPthlFq/ct691+Qustrb6oKpVKksvlao9ZWk+b1rNtp+X2x12tfnG1SDL7xCxtawAAYPPjYgQAAICh+oXb008
      "text/plain": [
       "<Figure size 2000x800 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nbins = 100\n",
    "vmax = 60\n",
    "\n",
    "fig, ((ax0, ax1)) = plt.subplots(nrows=1, ncols=2, figsize=(20, 8))\n",
    "\n",
    "a0 = ax0.hist2d(\n",
    "    rad_length_found,\n",
    "    eta_found,\n",
    "    density=False,\n",
    "    bins=nbins,\n",
    "    cmap=plt.cm.jet,\n",
    "    cmin=1,\n",
    "    vmax=vmax,\n",
    "    range=[[0, 0.6], [2, 5]],\n",
    ")\n",
    "ax0.set_xlabel(f\"$x/X_0$\")\n",
    "ax0.set_ylabel(f\"$\\eta$\")\n",
    "ax0.set_title(f\"found $\\eta$ rad_length_frac\")\n",
    "\n",
    "a1 = ax1.hist2d(\n",
    "    rad_length_lost,\n",
    "    eta_lost,\n",
    "    density=False,\n",
    "    bins=nbins,\n",
    "    cmap=plt.cm.jet,\n",
    "    cmin=1,\n",
    "    vmax=vmax * stretch_factor,\n",
    "    range=[[0, 0.6], [2, 5]],\n",
    ")\n",
    "ax1.set_xlabel(f\"$x/X_0$\")\n",
    "ax1.set_ylabel(f\"$\\eta$\")\n",
    "ax1.set_title(f\"lost $\\eta$ rad_length_frac\")\n",
    "# ax1.set(xlim=(0,4000), ylim=(-1000,1000))\n",
    "\n",
    "plt.suptitle(\"radiation length fraction and eta beginVelo\")\n",
    "plt.colorbar(a0[3], ax=ax1)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "rad_length_frac = found[\"rad_length_frac\"]\n",
    "# @ z = 9400.mm or 770.mm\n",
    "state = 1\n",
    "\n",
    "if state == 1:\n",
    "    slopex = found[\"ideal_state_770_tx\"]\n",
    "    slopey = found[\"ideal_state_770_ty\"]\n",
    "    x = found[\"ideal_state_770_x\"]\n",
    "    y = found[\"ideal_state_770_y\"]\n",
    "    qop = found[\"ideal_state_770_qop\"]\n",
    "elif state == 2:\n",
    "    slopex = found[\"ideal_state_9410_tx\"]\n",
    "    slopey = found[\"ideal_state_9410_ty\"]\n",
    "    x = found[\"ideal_state_9410_x\"]\n",
    "    y = found[\"ideal_state_9410_y\"]\n",
    "    qop = found[\"ideal_state_9410_qop\"]\n",
    "\n",
    "data = ak.zip({\n",
    "    \"rad_length_frac\": rad_length_frac,\n",
    "    \"x\": x,\n",
    "    \"y\": y,\n",
    "    \"tx\": slopex,\n",
    "    \"ty\": slopey,\n",
    "    \"qop\": qop,\n",
    "})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parameterisation for rad_length_frac:\n",
      "intercept= 0.0\n",
      "coef= {'1': 0.4397780582886272, 'x': -0.00033503537608265674, 'y': -4.789321464958182e-05, 'tx': 0.29235639572920724, 'ty': 0.04582618208302177, 'qop': -6.415064049320476, 'x^2': -1.7148780985188313e-05, 'x y': 6.213092137729638e-05, 'x tx': 0.0357792373581695, 'x ty': -0.03383507316198001, 'x qop': 5.411983741236269, 'y^2': -0.0001816058133305051, 'y tx': -0.060157737422332554, 'y ty': 0.3005860193729231, 'y qop': -0.5258168180545822, 'tx^2': -20.581151289333324, 'tx ty': 33.044369657406456, 'tx qop': -4029.9409136973454, 'ty^2': -125.62371943789807, 'ty qop': 353.08874965548245, 'qop^2': -3.455196781164501, 'x^3': 1.5371780097645036e-07, 'x^2 y': -1.2426481494632311e-06, 'x^2 tx': -0.0003102207951011151, 'x^2 ty': 0.0008294673378711963, 'x^2 qop': -0.0013130710648701907, 'x y^2': 1.3249247615395102e-07, 'x y tx': 0.0018733641082853674, 'x y ty': 0.0011459295349262105, 'x y qop': 0.0019703515773188207, 'x tx^2': 0.20678893428535178, 'x tx ty': -1.2036256002683805, 'x tx qop': 0.5976253324202965, 'x ty^2': -1.0782570211818505, 'x ty qop': 102.9459912574048, 'x qop^2': -399.16918894718305, 'y^3': 6.233244675968308e-07, 'y^2 tx': -0.0013994133271424403, 'y^2 ty': -0.001494031769823323, 'y^2 qop': 0.000739172959467016, 'y tx^2': -0.7290175680182052, 'y tx ty': 1.243439836769458, 'y tx qop': -105.55527881418587, 'y ty^2': 1.1837378224511292, 'y ty qop': -0.7657500120637686, 'y qop^2': 337.9089334089693, 'tx^3': -46.002810282557306, 'tx^2 ty': 453.57769058127064, 'tx^2 qop': -38.01772137810428, 'tx ty^2': -37.626648454968965, 'tx ty qop': 1.8221207378503694, 'tx qop^2': -0.40208255096589585, 'ty^3': -309.9803925733283, 'ty^2 qop': 15.736234283455065, 'ty qop^2': 0.6119680938280082, 'qop^3': -0.0011848472786331894, 'x^4': 1.4258219493967772e-08, 'x^3 y': -2.866883619390137e-09, 'x^3 tx': -4.622450717306492e-05, 'x^3 ty': -9.82760022623097e-05, 'x^3 qop': -0.0003703064328215433, 'x^2 y^2': -2.5736142106325133e-09, 'x^2 y tx': 0.00010255319807583874, 'x^2 y ty': 4.727919866809316e-05, 'x^2 y qop': -0.0006108544453127251, 'x^2 tx^2': 0.0554175462089006, 'x^2 tx ty': 0.16472559466997527, 'x^2 tx qop': 0.4783758545936343, 'x^2 ty^2': -0.029491297111002623, 'x^2 ty qop': 0.4072459238308874, 'x^2 qop^2': 1.274082272606741, 'x y^3': 1.019856199491187e-08, 'x y^2 tx': -4.0220230914655986e-05, 'x y^2 ty': -7.807169754414645e-05, 'x y^2 qop': 0.0009660222713137046, 'x y tx^2': -0.1655492341343554, 'x y tx ty': -0.026431127204681104, 'x y tx qop': 0.6220634937641649, 'x y ty^2': 0.09539033507193473, 'x y ty qop': -1.5617368820840163, 'x y qop^2': 9.235370076473936, 'x tx^3': -29.04572443637665, 'x tx^2 ty': -68.28508620169197, 'x tx^2 qop': -138.17259913877646, 'x tx ty^2': 35.45859347871301, 'x tx ty qop': -273.5835489656032, 'x tx qop^2': -13.43661305575325, 'x ty^3': -32.63704432223729, 'x ty^2 qop': 266.4268072225656, 'x ty qop^2': -19.06231959583929, 'x qop^3': 1.7675392248918753, 'y^4': -5.9619367220875574e-09, 'y^3 tx': 5.3697276143793715e-05, 'y^3 ty': 1.6576491020714457e-05, 'y^3 qop': -0.00017261869555795784, 'y^2 tx^2': 0.05014539344681879, 'y^2 tx ty': -0.07584181301500559, 'y^2 tx qop': -0.11650839062004614, 'y^2 ty^2': -0.015097852966895253, 'y^2 ty qop': 0.30517451156427916, 'y^2 qop^2': -6.650919392954738, 'y tx^3': 67.30430825382783, 'y tx^2 ty': -34.801544500318485, 'y tx^2 qop': -142.686910608299, 'y tx ty^2': 27.185653317937266, 'y tx ty qop': 439.1032304393983, 'y tx qop^2': -18.23933135674488, 'y ty^3': 4.265555264813383, 'y ty^2 qop': -131.08781810964697, 'y ty qop^2': 34.20534645433437, 'y qop^3': 0.11852127731033578, 'tx^4': 5649.699838454551, 'tx^3 ty': 329.9615549901753, 'tx^3 qop': -5.54031835759426, 'tx^2 ty^2': 584.460443200635, 'tx^2 ty qop': -0.7411483064667538, 'tx^2 qop^2': -0.03588219456135393, 'tx ty^3': 173.96226240131907, 'tx ty^2 qop': 1.0493642752682577, 'tx ty qop^2': -0.052660249222163585, 'tx qop^3': 0.001973073444663188, 'ty^4': 215.62061929608907, 'ty^3 qop': -0.5416245412488983, 'ty^2 qop^2': 0.09234573370962393, 'ty qop^3': 0.0002451739001199065, 'qop^4': 2.303474287996047e-05}\n",
      "r2 score= 0.07032743800675567\n",
      "RMSE = 0.1530771991719649\n",
      "Type  <class 'numpy.ndarray'>\n"
     ]
    }
   ],
   "source": [
    "lin_reg, features, xx0_test, xx0_predicted = fit_linear_regression_model(\n",
    "    data,\n",
    "    \"rad_length_frac\",\n",
    "    [\"x\", \"y\", \"tx\", \"ty\", \"qop\"],\n",
    "    4,\n",
    "    include_bias=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAikAAAHLCAYAAAAN9wgKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABwfElEQVR4nO3de3xTdZ4//lcpUAqVhoIiigIJKDcRToozOxcZbeJlvNsE9rc7t/2uNLrd7+yuqw3MqoDOENJxZ3dnJjMmdXbdmZ39jjR1xNs4JuDgjI4jzRGRiwJNAUVEIQ2lUFranN8f9WSaNPdzcu3r+XjkAcm5vXN6ybufy/tTJkmSBCIiIqICMybfARARERHFwiSFiIiIChKTFCIiIipITFKIiIioIDFJISIiooLEJIWIiIgKEpMUIiIiKkhMUoiIiKggMUkhKmCiKMLlco3a62dDMBiEy+WC3+/PdyhElASTFBq1RFGE1WqF2WyGTqdDc3NzvkMK8/v9MJvN0Ov1cDqdRXF9+X7q9Xro9fosR5gZl8uFOXPmwGKxIBgMqnbeQv5eIipmY/MdAFE+iKKIuro6dHV1AQCsVis6OjryHNWfabVatLa2oqysrGiuLwgCAKC5uTn8/3zz+/2oqamBRqMBADQ0NKCjo0PVJKLQv5eIihmTFBqVbDYbampqws/tdnseoykdhZKcyMxmM1pbW8NJCgBMnTpV1Wvwe4koe9jdQ6OSKIr5DoGyzGw25+TrzO8louxhkkKjisvlgtlsht/vD4+7MJvN8Hq94X2CwSAsFgusViuMRiOMRmPEdrfbjSlTpqCsrCz8AeX1emE2m1FWVgaz2Rw+j8vlgl6vh9vthtfrhV6vj9hnOPm68iOdLgm32x2O0+VyYcqUKbBYLOHtzc3N4fek1+tjDoZVcv1UiKIIs9kMo9EInU4Hq9Uace107hWAcJw6nW5EQuJ2u8PPLRZLzIQlEAiEr5noOvEk+15S+2sify8yKaJRRSIahbRaraTVake87vP5JI1GI/l8vvBrTqdTAiDZ7fbwaw0NDRKAiP06OjokAJLJZAo/N5lMEgDJYDBITU1Nks/nCx87/HwdHR2SRqORPB5P+DW73S4BkARBSPheWltbJa1WKwGQGhoapKamJkkQhPBxTU1N0vAfdY/HIwGIuJaS60eLdYzP55MMBkNEzHK88vVTvVeSJEmCIEhNTU3hYwFIACStVhs+p/y+Ozo6Io6V35fBYJDsdnvC66Qi1veSWl8TrVYbEb9Go5E0Gk3aMRIVKyYpNCrFS1IEQYj4MB3++vAPPPlDZniS0tXVFZGkSNKfP4yjP/zkD0mZyWSKed1UkwT5g1f+4B7OYDBEfLDJcQ7fV+n1kx0jCELEvZKkoQ9cAFJXV5ckSanfK3m/4eeTk4zhH+jJkhSn05nwOqmK972k9GsiCMKIeyGfU75nRKWO3T1En/H7/RBFMebgT7mZPtPpwMMHbsoCgUD4unLXQKbk8y9fvnzEttbWVvh8vvDz9vZ2AAhPwVXj+onI99Vms4W7RIZ3rcjxyBLdKwDYsWPHiO3y+dLpCqmtrR1xXTVrpyj9moiiCIPBEHFcU1MTJEmKeY+IShFn9xB9JtEHnPyBlo0CYPI5tVqt4nPF+vDSaDTQaDRwu914+umnR3xoqnn9WOT72traqsr55Pi9Xu+IhLLQZhcBmX1N5HvGZIRGO7akEEWJVeRL/rAYPtVULXKSMLy1QO3z6/V6+P1+tLa2oqmpKefXH/6vUiaTCQaDATabDV6vF8FgEHa7HU1NTVlLtNSW6teEVXFptGOSQvQZ+a/w4TN5ZHLiotPpVL+u/ME6vPlfTUajETU1NSM+CHN1ffn8brc75vZY9zuZ1tZWGAyGcNl+u91eVPVJkn1N5O/FeK1PTF5otGCSQqNSIBAY0XKg1WohCEJ4Sulw7e3t0Gg0aGhoAPDngmDD95P/n265dbkryeVyxTxWSfl2+b0M7zaQzye//2xeH0B4XIXVah3RpZbpukBykbampiY0NTUl7OaJjv/kyZMA1Gs5ivW9lEi6X5PoJM5qtWalRY+oEDFJIRpGrk46vJ6F3J3Q0tIS/mCRPxStVmu4DoY8qNbr9YYHoaby4aXRaMJ/Uev1eni9Xvj9/nAdEb/fn7RmifwhF/2BLH+Yud1uuFwuuFyu8HlFUQy3bii9fjrvz2w2o7m5GUajER0dHeEkJtUPevmDWz6Py+WKqIsik1u9nE5neHBwIsFgUNUuL6VfE7llyGg0wmw2h+up6HQ6jlWh0SPf04uIcml4TQx8VsNieG0KSRqaDipPyW1oaJAaGhpGTJ+VpKHpoHLdCnnqqFarDdf48Pl84anLWq1W8ng8UldXV8T1h08xdTqd4doagiCE62Q0NTWNmEY73PCaHFqtdsTUWqfTKWk0mvC5JGloyu7wuNO5vhx7rEf0exs+VdZut8eNM517JccV6/parTbimoIgSBqNJlw7pbW1NTz1WZ4WHX2d4fck0XtN9L2k1tektbU1fF8EQRjxvUpU6sokSZKylwIREalLFEU8/fTTWLt2LQKBQEQLSGtrK3Q6XdyxHkRUXDgFmYiKhjwrpqurKzyNdzitVpvRQFwiKkwll6QEg0HYbDYAqa9GKhea0mq1CAaDMBqNMJlM2QyTiDIgD05evXo11q5dGx4b5Pf74fV60dHRUVSzfIiKnTw+rKurK+J11T5X893fpCaPxxNe/0Pug05GXvdj+JiDWH3IRFQY5LFAGDY+hOM1iPJDq9WOWE9Kzc/VkhyTUlZWhoaGhpRKmMuzMDweT/g1l8sFi8WCErw1RCVDHotSLAXciEqNXFagvb09oiVFzc/VUT0FORgMRkwXlQ2vUUBEhUmj0TBBIcoTr9eLqVOnjqhRpPbn6qhOUuRFvaJ/0ck3fXgWSEREREOcTmfMWXRqf66W3MDZdMiD8OIVRkpWevrEiRN48cUXMWPGDEyYMCHjOMaPH4/x48dnfDwREWVXf38/+vv7Mz7+3LlzOHbsGG655RZMmzZNxcgiHTlyBCdOnFAcLwDMmDED8+bNG/G61WqNO0Bd6edqtKwlKd3d3QgEAvD7/dBqtZg9e3a2LpWxjo4OAPEXjUtWDvzFF1/Et771LZWjIiKiUvXUU0/hm9/8ZlbOfeTIEcydNQvnVTrfuHHjcPDgQVx++eXh10RRxNSpU+N2tSr9XI2mepKybds2bNq0CVu3bgUASJKEsrIyAEPNPc3NzbjuuuvUvmxG5LLZ8UphJ+vvnjFjBgDghz/8Ia6++uqMYvj617+Op59+WlFLSn19Pdra2vJ2fE9PD1asWIHt27ejqqoqLzGocQ7eB3WOL5X7oPQcatwHpTGocbzSc5TKfejv78eqVavwi1/8IuVjzp0P4edvHcXRYC+uqQqi+ZEHw58b2XDixAmcB2ACUANgUMm5AGw5fx4nTpyISFJsNlvchS8B5Z+r0VRNUtra2mCz2bBq1aoRi2AFAgH4fD48+OCDuPfee3HPPfeoeemMyDcrXmaX7GbKXTxXX301rr322oximDRpEj7/+c9ndKyssrIy4QJr2T6+u7sbALB06VJMnjw5LzGocQ7eB3WOL5X7oPQcatwHpTGocbzSc5TSfZg0aVLKv+tPnT2Pdc/txsDUCvzkW4vw6YGdaAYUDQ1I1SUALlV4jkkxXrNarTAajTEXVpX/Vfq5Gk3VJKW9vT08aCaWuro6NDU1Yc2aNWpeNmPyaOPoPjL5uV6vz3lMRERU3D45fQ6PPLsHZ/oHYLt7CeZMm4RPD+Q7KuW8Xm/cxUZ1Oh0EQQj3oqj1uarq7B55+fpksjloKB0ajQaCIIwYbSyX1V65cmU+wiIioiL1YddZWN27MBAKwV4/lKDkWjmGWiCUPMpjnNfn80GSpIhHU1MTNBoNJEmCz+dT/XNV1STl4MGDWLt2LQ4fPhxz+6FDh/D444+HB9ZkQ6JBOX6/HzqdLmJtj5aWlvDS9DK73Q673c7l0ImIKGUHP+nBmrZ3MXH8WGyqX4JLNJV5iWMsgHEKH0q6WdT8XE07jm3btuH666+Pue2JJ56
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nbins = 100\n",
    "vmax = 40\n",
    "\n",
    "a0 = plt.hist2d(\n",
    "    xx0_test,\n",
    "    xx0_predicted,\n",
    "    density=False,\n",
    "    bins=nbins,\n",
    "    cmap=plt.cm.jet,\n",
    "    cmin=1,\n",
    "    vmax=vmax,\n",
    "    range=[[-0.2, 1.0], [-0.2, 1.0]],\n",
    ")\n",
    "plt.plot([-0.2, 1], [-0.2, 1], marker=\"\", alpha=0.8)\n",
    "plt.xlabel(f\"True $x/X_0$\")\n",
    "plt.ylabel(f\"Predicted $x/X_0$\")\n",
    "plt.title(f\"found rad_length_frac\")\n",
    "# ax1.set(xlim=(0,4000), ylim=(-1000,1000))\n",
    "\n",
    "plt.colorbar(a0[3])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test = ak.zip({\n",
    "    \"tx\": slopex,\n",
    "    \"ty\": slopey,\n",
    "    \"x\": x,\n",
    "    \"y\": y,\n",
    "    \"qop\": qop,\n",
    "})\n",
    "test = np.column_stack([ak.to_numpy(test[feat]) for feat in features])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "tuner",
   "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}