tracking-parametrisation-tuner/thesis/MatchingParameterisation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "import uproot\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import mplhep\n",
    "import awkward as ak\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "\n",
    "mplhep.style.use([\"LHCbTex2\"])\n",
    "plt.rcParams[\"savefig.dpi\"] = 600\n",
    "# %matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "tree = uproot.open({\n",
    "    \"/work/cetin/LHCb/reco_tuner/data_matching/parameterisations/sample4/param_data_B_default_thesis_selected.root\":\n",
    "    \"Selected\"\n",
    "})\n",
    "allcols = tree.arrays()\n",
    "allcols = allcols[allcols.isElectron == 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "zMagnetParamsEndT = [\n",
    "    5283.043384741342,\n",
    "    28.45437483916609,\n",
    "    0.016308934185688337,\n",
    "    -1668.3602223088105,\n",
    "    301.07580662497,\n",
    "]\n",
    "dSlope_xEndT = allcols.tx_l11 - allcols.tx\n",
    "dSlope_xEndT_abs = abs(dSlope_xEndT)\n",
    "x_EndT_abs = abs(allcols.x_l11 + allcols.tx_l11 * (9410.0 - allcols.z_l11))\n",
    "\n",
    "zMagTerm = zMagnetParamsEndT[0] + zMagnetParamsEndT[3] * allcols.tx * allcols.tx\n",
    "\n",
    "z_mag_xEndT = (\n",
    "    zMagTerm + dSlope_xEndT_abs *\n",
    "    (zMagnetParamsEndT[1] + zMagnetParamsEndT[4] * dSlope_xEndT_abs) +\n",
    "    zMagnetParamsEndT[2] * x_EndT_abs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "zmag = (allcols[\"x\"] - allcols[\"x_l11\"] - allcols[\"tx\"] * allcols[\"z\"] +\n",
    "        allcols[\"tx_l11\"] * allcols[\"z_l11\"]) / dSlope_xEndT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABMEAAAOWCAYAAADm+IXMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAAC0EklEQVR4nOz9T2wb2b3nf3/KyfwQYHClkgwMIEhAWsV4tmOR8l09+AGxip1s0yal3t1FWqQ9SyUmrV7NqmWq01pOTMq9uLsrkXa2GbPkBn74rdJiubOdNEsKHhsCZmCyzGCAYAatehZ6yIj6S4kUWZTeL4BoUTp1zpeU3AA/+J5TRhAEgQAAAAAAAIBr7NagCwAAAAAAAACuGiEYAAAAAAAArj1CMAAAAAAAAFx7hGAAAAAAAAC49gjBAAAAAAAAcO0RggEAAAAAAODaIwQDAAAAAADAtUcIBgAAAAAAgGuPEAwAAAAAAADXHiEYAAAAAAAArj1CMAAAAAAAAFx7hGAAAAAAAAC49gjBAAAAAAAAcO0RggEAAAAAAODaIwQDAAAAAADAtUcIBgAAAAAAgGuPEAwAAAAAAADXHiEYAAAAAAAArj1CMAAAAAAAAFx7hGAAAAAAAAC49gjBAAAAAAAAcO0RggEAAAAAAODaIwQDAAAAAADAtUcIBgAAAAAAgGuPEAwAAAAAAADXHiEYAAAAAAAArj1CMAAAAAAAAFx7hGAABs73/UGXAAAAAAC45gjBAAyE67rKZrMaGxvT4uLioMsBAAAAAFxzPx50AQDCx3VdbWxsyHVdeZ6nWq0m3/dlmqYsy9L4+Lii0aji8bhs2z52fTKZlOd5qlQqre/5vi/HcbSxsSHHcej+AgAAAAD0lREEQTDoIoBeKxQKSqfTF74uGo22BTeXZRjGha8xTVP1ev3EnyWTSZVKpQvPWS6XTwypTuL7vlZWVrS6unrhdaLRqGzbViQSUblcVqlUOvZeNoMx13WPXZ9IJFQsFi+8LgAAAAAAnWI7JK6lVCqler2uSqWiRCJx5ljbtlUul1vjeyEIAlWrVeXzeZmmeeo4y7KUz+dVrVZPDcAkqVgsql6vdxRqZTIZVSoV1ev1jgOw5rbEowFYM5xqztd8XeVyWalUqvXaXNfV6uqq0ul0K6w72unVnCcIAlmW1VFdAAAAAAD0Cp1guBHS6bQKhcKx7/ejA8nzPEUikRN/Vq/XzwzJThOLxU7sqMrn80qlUheqLR6Py/O8tu+nUinl8/mO5lhdXVU2mz32fcuyVK1WT7zm6O+DTjAAAAAAwFWjEww3wmlbIy+zZfKiLMtSNBo99v1EInGpAEw6ve6LBGCO4ygSibQFYKZpqlKpdByASQedZ9Vq9dhrqdVqp15z2dcNAAAAAMBlEYLhRjht+12/tuWdtE43a5907UWCJdd1FY/Hj12/tbV1YmDXST07OzttNXDwPQAAAAAgTAjBcCOcFhCNj4/3t5AQ8DxPc3Nzx75fLBYvFYA1NbvIrqNSqXRsyyj+wXGcQZcAAAAAAOciBAP64KSw7fbt2z2dr9NAL5lMHuvSSqVSHR+ifxbLspTJZLqep19c11U6nVYkEpFhGBobG1MkElEsFmvdkbNUKmlxcfHcoMdxnNb11z0U8jxPhUJB8XhchmEc6yoctJv0uwAAAADQuR8PugAA/eM4zokH6udyuZ6tsby83LrLpO/7oTz/y/d9LS4utu5kaVlW6y6inufJdV25rtv6efOaszTDRd/3lU6nT70pwDDzPE/JZFKe54V6u+tN+F0AAAAAuDg6wYAb5KQD9bs5oP8kpmle6ID+QWh2eUkHd9SsVqsqFosqFouqVCqqVqvHOuPOC1JuwtZay7JUqVRUr9dD3fF3E34XAAAAAC6OEAy4IRzHOfFcq+Xl5Z6vlUwmJZ19h8izOI6jZDLZtk0xHo+rUCh0Xdvq6mpri1wmkzkxsLMsS+Vyua1D7rwzwfL5vCzLkmVZKhaLXdcZdmHbAnnYTftdAAAAAOgM2yGBG+K0MKCbw/BPY9u2crnchTtyfN/X3NzcsS2bvu/LcRw5jqNcLtfVIf4rKyutr0/qjDssk8no/fv3Wl1dPTfQs237Rm27G1S3VTqdVj6fP3PMTftdAAAAAOgMnWDADXHSAeG9OAz/NJlM5kLbLD3P0/T09Ilnlh0dF4vFLnXg+WXOssrlcopGo9wdMgQKhUJPugEBAAAA3EyEYMAN4HneiSHOVXSBXYbv+63DzBOJhPL5vCqVisrl8qlnT8Xj8QsHU0cDsE6DtOXl5VAfBH8TeJ53buceAAAAAJyF7ZDADXBaWBSJRPpcycmaYVS5XD7WnWbbttLp9ImhVzKZVKVS6Xgdy7LanmezWc3Pz5/bsda8c2RY73Z53fm+H+ozyAAAAAAMBzrBgBvgtBAsTHfROykAa2oeVH+U67oX2hZpmmZbiOX7vmKxWEcdZUEQEIANQPOcOLajAgAAAOgWIRgwINlsVoZhXOoRi8UutNZpW/nCEuqYpnnu+WSWZbXdrbHppO+d5eh4z/MUiUR6ctaU67pKp9Pndtj5vq9CoaBYLNa2bnPL39jYmAzDUCQS0erq6rHrPc9TNptt3T0zEokonU6f+HturnP0b+ioeDzemu/wo5fbQAuFQts6Y2NjisViJ75G6eD9POmcuMP1nfZed/q7aPJ9X6urq4rH4633v1lfNpvtKITzPE+rq6uKRCJtNR/+HUQikdbWXwAAAAB9FgA3hKRjj3q93pe1U6nUiev38mFZ1oXXL5fLfXn9R2UymbY6bNvu+FrTNI+9jouybfvE9yMajQaVSuVCc1UqlSCVSrXVZZrmsXH1ej3I5XJBNBptWzOfzwdBcPw9OfxIpVKtec4aZ1nWqX/TR/8GTpPL5Tr6N1KpVDr+HRSLxdb7UiwW2+awLOvU2uv1elCpVI7VVKlU2h6H5+vkd3FUPp9vjc/lckGlUgmq1WpQLBbbfl+JROLYtfV6PchkMq3XcbjGer1+7Pd9uK5qtXpubQAAAAB6h04wYEBSqZQqlcqlHvl8/kJrndbxFZZulIt0pJ3UMXbeHSWPKpfLrXO+js4Ti8U67tTxPE8bGxsddxtFo9Fjh7tXKhXFYjG5rqtisahqtap6vd52Q4BCoSDXdRWPx+U4jsrlsur1uur1usrlcuv98zxPKysrJ66dTCY7qvGk96UbjuO01l5eXm6bPxqNqlgsSjqofXFxse1a0zQVjUaP3cCh+b3DP7vo76Ipm80qnU7LNE1Vq1VlMhlFo1FZlqVEIqFKpdKquVQqKRKJtP1t1Go1xePxY3+X29vbmp6ebm3lrVaryufzrd+V7/sc9A8AAAD026BTOKBfFLJOsFwud+n5jnbh6JxOsKOdNM1Hswup3452M53UYXOaZldRL17Hae/LZeY9PNd53UdHu4NO68g72rF2uCPssMPvyWlrl8vljjq36vV6TzvBEonEuX/z583Tae1Nnf4uDneAnfe7PtzpdVrn4uEaTdM88fd69O+3X/8PAgAAAEAnGHAjhL0T7CKO3uFRuvzryGQyqtfrp3Y/Ne9K2cn8J9V1msM3JEilUqeeh3b0joindQAert/3/dD+Xk/7Ozz8/V4cgN/p7yKbzbbWT6VSHY2VDrrbTrohw+HXsb6+fuLv9ejf2vb2dke1AgAAAOgeIRhwA5wWClSr1T5X0r2TgpT37993NV+xWFSlUjm27U46CDzm5uZ6Gix1uv3z8O/tIltGw3QnxfX1deVyORWLxXODJql/wWyhUGit1Ulolkql2n4HF70hw2GH1wvT7woAAAC47gjBgBvgtE6jk7pZwu6kwOL27dtdzxuNRk89b8113WPnVfXDZe/eWavVeltIF0zTVCaTOfUMtoue59Yr5XK59XWnnWOH/x0N478dAAAA4KYjBANuiJO6nDzPC+3WuYu4yFbE86RSKVWr1WMBVKlUCnXwcdnArN88z1OhUJDjOCceet/POi7q6N/Zdfi3AwAAANwkhGDADbGwsHDi98Mc7JzmaODTyxCsOV+lUjm2zkXvyol/aN7dMp/Pn3kOWr8cDrA6DcS
      "text/plain": [
       "<Figure size 1200x900 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(\n",
    "    zmag,\n",
    "    bins=50,\n",
    "    alpha=0.5,\n",
    "    density=True,\n",
    "    color=\"#107E7D\",\n",
    "    label=\"True\",\n",
    ")\n",
    "plt.hist(\n",
    "    z_mag_xEndT,\n",
    "    bins=50,\n",
    "    alpha=0.5,\n",
    "    density=True,\n",
    "    color=\"#F05342\",\n",
    "    label=\"Predicted\",\n",
    ")\n",
    "plt.xlabel(r\"$z_{mag}$\")\n",
    "plt.ylabel(\"Number of tracks (normalised)\")\n",
    "mplhep.lhcb.text(\"Simulation\", loc=0)\n",
    "plt.legend(loc=\"upper right\")\n",
    "# plt.savefig(\n",
    "#     \"/work/cetin/LHCb/reco_tuner/thesis/zmag_new_parameterisation.pdf\", format=\"PDF\"\n",
    "# )\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "weird = zmag < 5200\n",
    "True in weird"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABMAAAAOWCAYAAAAJOu7yAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB2sUlEQVR4nOz9UWyjaXofej5U10zaGYND1QAGBHKBDJWDnQDnYpqsWSx8ruKiYuwez7GRkaqz9s4M9iQljQ0jCzCB5IZxkAQLuyLB4c0C4xE7ydkzPrNAW0p2fDwHMEZsXy/QJbZzscAAjti+EEEgC1exCTvutafJvahIllQSRUkUP/Lj7wcII5Yevnw43W+zvr/e9/0yg8FgEAAAAACQUgtJNwAAAAAA90kABgAAAECqCcAAAAAASDUBGAAAAACpJgADAAAAINUEYAAAAACkmgAMAAAAgFQTgAEAAACQagIwAAAAAFJNAAYAAABAqgnAAAAAAEg1ARgAAAAAqSYAAwAAACDVBGAAAAAApJoADAAAAIBUE4ABAAAAkGoCMAAAAABSTQAGAAAAQKoJwAAAAABINQEYAAAAAKkmAAMAAAAg1QRgAAAAAKSaAAwAAACAVBOAAQAAAJBqAjAAAAAAUk0ABgAAAECqCcAAAAAASDUBGAAAAACp9iCpF/7DP/zD2N3djb/39/5e/MN/+A+TaoMJ+9znPheffPJJvPHGG/FTP/VTSbcDAAAAXPCf/tN/ik8//TTefPPN+PM///Ok2xmLzGAwGEz6RT/88MMol8uRyWQiIuLo6Cj+1t/6W5NugwS88cYb0e/3k24DAAAAuMbCwkJ8+umnSbcxFomsANva2oqIiMFgEJlMJlqtlgBsTpwEYAsLC7G0tHSnsTqdztjGSsKs9x8x++9h1vuPmP33MOv9R8z+e5j1/iNm/z3oP3mz/h5mvf+I2X8Ps95/xOy/h1nvP2L238Os9x8x++9hnP2fjPXGG2+MqbvkJRKANRqNyGQyMRgMYn19PX7mZ37m0rrf+q3fikqlEl/+8pcn2yD35qd+6qei3W7H0tJSHB8f32msQqEwtrGSMOv9R8z+e5j1/iNm/z3Mev8Rs/8eZr3/iNl/D/pP3qy/h1nvP2L238Os9x8x++9h1vuPmP33MOv9R8z+exhn/ydjpenookQPwd/d3Y3f/u3fvvLnm5ub8fz58wl2BAAAAEDaJBaA7e7uxtOnTyf2er1eb2KvBQAAAMD0SCQAKxaL8ejRo4m+5tra2kRfDwAAAIDpkMgZYKVSKd5///2Jnu1lK+V06XQ6USgUhtZUq9WoVqsT6ggAAADSr1arRa1WG1rT6XQm1M3kJBKAra+vx9/7e38v/vRP/zTW19eH3gHy448/vtP2xVarFZubm9Htdm89BuPX7/ej3W4PrbFtFQAAAMar1+tdez2eRokEYJVKJVZXV2N7ezu2t7eH1m5ubsbm5uadXm8wGEQmk7nTGIzXKLdlzWazE+oGAAAA5kM2m418Pj+0ptPpRL/fn1BHk5FIAPbxxx9Hs9mMwWBwbe0oNcyeWb2tLAAAAMyyUY4bKhQKqVsllkgA9uzZs2i1Wqersj7/+c/f2xbFTCYjRGMqVavV6PV6M73Sbdbfw6z3H5GO9zDrZv2fwaz3HzH772HW+0+DWf9nMOv9R8z+e5j1/tMgDf8MZv09zHr/Eel4D1wtM0ggHXr06FF8+OGHMRgMolwuxwcffBAff/xxNBqN+OCDDyIiotvtRr1ej1KpFMVi8dav1Wq1otlsRiaTiU8//XRcb4FbOkmR8/n8nVeAjXMsIBnmMcw+8xjSwVyG2ed6e7hEVoC1Wq2IiNjY2Ijf/u3fjohXq8C+9rWvxde+9rXTunq9Ht/61rfiH/2jf3Sn12s0GvGzP/uzdxoDAAAAgNm0kMSLdrvdWF9fPw2/7lulUom33nprIq8FAAAAwHRJJACLeLX6K82vBwAAAMB0SCwA+/KXvzzR13v69OlEXw8AAACA6ZBIAFYsFuM//If/MNHX7PV6E309AAAAAKZDIgHY6upq/OZv/uZEX3NtbW2irwcAAADAdEgkAFtfX4+9vb34N//m30zsNZ8/fz6x1wIAAABgejxI4kWLxWJ87WtfOw3CVldXo1KpxMOHDyObzZ6r/fjjj2+9ffHFixfR7XZjc3Mzut3uGDpnmlSr1ej1eq/9OwPMDvMYZp95DOlgLsPsM4+HywwGg0FSL14ul+PDDz+MTCZz6c8Hg8GVP7uJk3E+/fTTO4/F3RQKhWi325HP5+P4+DjpdgAAAIAL0njtnthdICMi/vAP/zC+9rWvxWAwuPQrIq782U2+AAAAAJhfiQZgn//852Nvby9++MMfxuPHj+/lNcaxggwAAACA2ZXIGWAXVSqVqFQq8fHHH0ej0YhWqxVHR0dRr9ejVCrFo0ePbjzmixcvIiKi1WpFs9kcd8sAAAAAzIhEzwC7zsLCQtTr9fhH/+gf3WmcZrMZX/nKV5wBNgVO9hEvLCzE0tLS0NpqtRrVanVCnQEAAED61Wq1qNVqQ2s6nU70+/1UnQE2FSvA7lupVIq33nor6TY4o9/vR7vdHlpz27t/AgAAAJfr9XrXXo+n0VQHYKurq1EsFscy1sbGxljGYTxGWQHm1q0AAAAwXtlsNvL5/NCakxVgaTLVWyBJnzTeShUAAADSJI3X7oneBXKS/uRP/iTpFgAAAABIwNwEYGtra0m3MLNarVasrKwk3QYAAADArcxNANZqtZJuYep0u93IZDLXfi0vL4/tLDYAAIBRDAYRn/aHfznQBxjVVB+CPw69Xi82Nzej2+0m3crUqdfrI9dubW3dYyfAfRsMIvrX/AVxIRORyUymHwCA6/QHEQ9+fXjNj38j4g1/fwFGkEgA9tFHH8WTJ0+i1WrFxsZG/OZv/uZrNY8ePYqPPvroTq9zEnoNBoPIuKp7zbNnz0aqq1QqVoDBjPMXSAAAYJ4lEoCtra1Fs9mMiIjt7e149OhR/P2///fP1ZTL5dMaxq9er0e3243Nzc1rz/d69OjRhLoCkvS9DyO+UU66CwAAgPFLJABrNpuRyWRi8F82bF92Pte3vvWtePfdd19buTUYcZP32eeN+px5sr29HcViMba3t5NuBZiA7314fc039yPeWIj4pbfuvx8AAIBJSiQAe/r0abz77runj1dXV1+reeuttyKXy8XHH38cg8EgcrlcPHz4MHK53Eiv0e1248WLF87+usT+/n60Wq3Y3d1NuhVgAj7+JOJXfm+02m99P+Krfyci++a9tgQAMBZWsAOjSiQA293djZWVlXj+/HlsbGzE3/pbf+vSuqdPn8Zv/dZvRbPZjC9/+cu3fr2tra34rd/6rVs/P22ePXsWuVwunjx5knQrwAT8TjPiz/9qtNo/+8uI7zYjfvWn77cnAIDrWMEOjFNmMMX7Az/88MN48uRJ/PEf//Gdx3r48GG8ePFiDF3NtmazGeXyX/+KpFgsRqVSiZWVlUtX4o1boVCIdrsd+Xw+jo+P7/31gIif/bcRP7zBf0Z/9r+K+IP//v76AQC4zsefROR/c7Rf4v3kZyPa71jBDuOUxmv3haQbGOatt96Kz3/+82MZy10MX9na2jr3uNVqRb1ej7W1tchkMuduUACkw8u/uN96AIBxu80KdoBhpjoAi4h4//33xzLO2TPH5lWr1YpGozG0Zn9/P8rlcmxsbEyoK+C+Lf7E/dYDAIzb7//oZvU/uGE9MH+mPgAbxwqwX/7lX4633rIpvFgsxu7ubmxvb8f6+vrQVXH1ev3cVklgdn31Szer/7kb1gMAjJsV7MC4TX0ANg71ej3+5E/+JOk2psL6+npsbm7G7u5uHB0dxcuXL2N7e/vSu2s2m81YWVmZfJPAWH29FPG5z4xW+5OfjfhG6X77AQC4jhXswLglchfISXr//fdjMBjEv/t3/y7+yT/5J0m3M3VyuVxsbm7G5uZm7O/vx9OnT6Pb7Z7+vNFoxM7OTmxubo71dTudThQKhTuPU61Wo1qtjqEjSK/Pvxnx7Z9/dZek63znFxwgCwAk76tfutlNfKx
      "text/plain": [
       "<Figure size 1200x900 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "array = {\"z_mag\": zmag, \"z_mag_pred\": z_mag_xEndT, \"diff\": zmag - z_mag_xEndT}\n",
    "df = pd.DataFrame(array)\n",
    "sns.regplot(data=df, x=\"z_mag\", y=\"diff\", x_bins=30, fit_reg=False)\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": []
  }
 ],
 "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
}