tracking-parametrisation-tuner/parameterisations/notebooks/bend_y_params.ipynb

650 lines
365 KiB
Plaintext
Raw Normal View History

2023-12-19 13:00:59 +01:00
{
"cells": [
{
"cell_type": "code",
2024-02-29 15:54:19 +01:00
"execution_count": 13,
2023-12-19 13:00:59 +01:00
"metadata": {},
"outputs": [],
"source": [
"import uproot\n",
"import awkward as ak\n",
"import matplotlib.pyplot as plt\n",
"import seaborn as sns\n",
2023-12-19 13:00:59 +01:00
"import numpy as np\n",
"import mplhep\n",
"\n",
"from sklearn.preprocessing import PolynomialFeatures\n",
"from sklearn.linear_model import LinearRegression, Lasso, Ridge, ElasticNet\n",
"from sklearn.model_selection import train_test_split\n",
"from sklearn.pipeline import Pipeline\n",
"from sklearn.metrics import mean_squared_error\n",
"\n",
"mplhep.style.use([\"LHCbTex2\"])\n",
2024-02-29 15:54:19 +01:00
"input_tree = uproot.open({\n",
" \"/work/cetin/LHCb/reco_tuner/data/tracking_losses_ntuple_B_BJpsi_def_selected.root\":\n",
" \"Selected\"\n",
"})\n",
2023-12-19 13:00:59 +01:00
"array = input_tree.arrays()\n",
"\n",
2024-02-29 15:54:19 +01:00
"array[\"yStraightOut\"] = array[\n",
" \"ideal_state_770_y\"] + array[\"ideal_state_770_ty\"] * (\n",
" array[\"ideal_state_10000_z\"] - array[\"ideal_state_770_z\"])\n",
"array[\"yDiffOut\"] = array[\"ideal_state_10000_y\"] - array[\"yStraightOut\"]\n",
2024-02-29 15:54:19 +01:00
"array[\"yStraightEndT\"] = array[\"ideal_state_770_y\"] + array[\n",
" \"ideal_state_770_ty\"] * (9410.0 - array[\"ideal_state_770_z\"])\n",
"array[\"yDiffEndT\"] = array[\"ideal_state_9410_y\"] - array[\"yStraightEndT\"]\n",
"\n",
2024-02-29 15:54:19 +01:00
"array[\"dSlope_xEndT\"] = array[\"ideal_state_9410_tx\"] - array[\n",
" \"ideal_state_770_tx\"]\n",
"array[\"dSlope_yEndT\"] = array[\"ideal_state_9410_ty\"] - array[\n",
" \"ideal_state_770_ty\"]\n",
2023-12-19 13:00:59 +01:00
"array[\"dSlope_xEndT_abs\"] = abs(array[\"dSlope_xEndT\"])\n",
"array[\"dSlope_yEndT_abs\"] = abs(array[\"dSlope_yEndT\"])\n",
"\n",
2023-12-19 13:00:59 +01:00
"\n",
"def format_array(name, coef):\n",
" coef = [str(c) + \"f\" for c in coef if c != 0.0]\n",
2023-12-19 13:00:59 +01:00
" code = f\"constexpr std::array {name}\"\n",
" code += \"{\" + \", \".join(list(coef)) + \"};\"\n",
2023-12-19 13:00:59 +01:00
" return code"
]
},
{
"cell_type": "code",
2024-02-29 15:54:19 +01:00
"execution_count": null,
"metadata": {},
2024-02-29 15:54:19 +01:00
"outputs": [],
"source": []
},
{
"cell_type": "code",
2024-02-29 15:54:19 +01:00
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
2024-02-29 15:54:19 +01:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAABMIAAAOVCAYAAACLW0xhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAACgnElEQVR4nOz9T2xbd57ne39IO4EzSdOUXRhAEIEuH1UGLsxiRofyLDKraZG3FhNP0h3SLiRjG7jdplJ5gp4HGpQYV5YXiYsalIB75yIdMlUXGLuzsMmqTpA0UB0xmdUsBpaOq3fGVJNOA9JDYICYNDu5MZKY51lwyBb1hzykDnkOyfcLIIqSvvydr/4lpU++v98J2LZtCwAAAAAAAJhwQa8bAAAAAAAAAEaBIAwAAAAAAABTgSAMAAAAAAAAU4EgDAAAAAAAAFOBIAwAAAAAAABTgSAMAAAAAAAAU4EgDAAAAAAAAFOBIAwAAAAAAABTgSAMAAAAAAAAU4EgDAAAAAAAAFOBIAwAAAAAAABTgSAMAAAAAAAAU4EgDAAAAAAAAFOBIAwAAAAAAABTgSAMAAAAAAAAU4EgDAAAAAAAAFOBIAwAAAAAAABTgSAMAAAAAAAAU4EgDAAAAAAAAFOBIAwAAAAAAABTgSAMAAAAAAAAU4EgDAAAAAAAAFOBIAwAAAAAAABTgSAMAAAAAAAAU4EgDAAAAAAAAFOBIAwAAAAAAABTgSAMAAAAAAAAU4EgDAAAAAAAAFOBIAwAAAAAAABTgSBsDNRqNaXTacXjcc3PzysQCGh+fl7JZFK5XG7gddfW1hSPxzUzM9OxZrFYdLF7AAAAAAAAfyAI87m1tTXNzMwol8spHo8rm81qa2tL6XRa5XJZy8vLmp+f7yu8KhaLmpmZUTqdliTl83mVSiVlMhlZlqV4PK54PK5arTakzwoAAAAAAGD0ArZt2143gYMlk0kVCgWZpqmtra0Da5aXl9tTYVtbWzJNs+uaxWJR8XhckpRKpZTNZvfVRKNRWZYlwzC0tbWlcDh8tE8EAAAAAADABwjCfGptba09sVWtVruGUfPz8yqXywqHw6pWq4fW1Wo1nTlzRrVaTYZhqFQqHVhXLpc1Pz8vSYrFYtrY2Bj8EwEAAAAAAPAJtkb6VCsEM02z50RWIpGQ1Ay6CoXCoXXJZLK93bG1/kEMw2ivWSwWj3QOGQAAAAAAgF8QhPnQ7vO+DMPoWX/u3Ln28zt37hxYUy6XO9a9cOFC1zUvXrzYft4tNAMAAAAAABgXBGE+VC6X288ty+rrtadPnz7w/ZlMpv08Fos5njKTek+aAQAAAAAAjAOCMB86depU+3m5XO4Ixg6yewrssAmy3dsbex2of9Bat27dcvQaAAAAAAAAvzrudQPYb2+Ytby83PXA+ta0Vjgc7pjkatk7VbZ7K2U3pmm2Qzi3JsKefvppPXr0SMeOHdM//af/1JU1AQAAAADA+Pif//N/6vHjxzpx4oS++uqrkV6bIMyHTNOUaZrtAKtYLCqZTCqfz++rzeVy7bBq9/bH3XafDSY5O3fsoDrLshxPkx3m66+/lm3bajQa2tnZOdJaAAAAAABgfH399dcjvyZBmE+99957ikaj7bcLhcK+MKxYLGp5eVmSlM1mlUqlDlxr7wH6vc4Ha9l73tjm5uaRgzDbttvP5+bmjrTWqFUqFTUaDQWDQc3OznrdTl/Gtfdx7Vuid6/Quzfo3Rvj2vu49i3Ru1fo3Rv07o1x7X1c+5bo3SutwZjdGcGoEIT5lGmayufzSiaT7fcVCgXNz88rn8+rWCwqnU7LMAxls1nFYrFD19p7xtigE2GlUqmPz+BgwWCw/Yu6vb195PVGKRKJaGdnR7Ozs/Q+IuPat0TvXqF3b9C7N8a193HtW6J3r9C7N+jdG+Pa+7j2LdG7V44dO9bOBkbN0yDs2Wef9fLyHQKBgP7H//gfXrfRIZFIKJvNtqe+pGao1ZoUW11dPXQ75G69Dtt3qlarubKOJDUaDUUikSOvs7KyopWVFRc6AgAAAAAA3ayvr2t9ff3I6zQaDRe6GYynQVipVFIgEPBkFK6ldf1AIOBZD920tjvuDsNaisWiarVaz62OgwZYe9d98ODBQOscxo0zwur1ugudAAAAAACAXur1+tif9z31WyO9DOGcSqVS2traUi6X63i/ZVk6c+aMPv300yOf3eWEmxNhkjtnhIVCIRc6AQAAAAAAvYRCIVf+lvcyTPM8CLNtW4lEwvG5VW6r1Wr7Aia/WV5eVi6XUyqVUrFY7NjqWKvVFI1GtbGxceg5YeFw2JUQy+kh+06M4xlhAAAAAABMM7eOJ2qdEeYFz4OwXC6nP/uzP/O0h0QioR/96Eee9nCYeDyuYrHYcR5YKxjbW3dYGHbq1ClXgrBTp04deQ0AAAAAAACvjP54/j263e1wVM6dO+d1CweKRqMqFouKxWIdh+Jns1nl8/l99clk8sDAa9BJrr1ruTkRBgAAAAAAMGqeB2FMGR0smUzKsixJzeBrr0Qioa2trY5wqlarKZ1O76tdXFzseNvpdNjew/Hn5+cdvQ4AAAAAAMCPPN0a6dc7NXqtXC6rUChIkkzTPPT8NNM0tbW11RFQ5XI5ZTKZjoAsGo3uW9/J4fqlUqnjbTem9/7gD/5ADx8+1B/8wR8ceS04t7Kyonq9PnY3FxjXviV6R//G+etO794Y597H1Th/zendG+Pc+zgb56/7uPY+rn1L4937OPMyGwjYHt42MRgMqlaref4D9/DhQ83MzHh2UNte6XRaa2trkpp3jDxoImy3XC6n5eXl9tt7zwqzLKsjDMvn80okEj372H0WWTgcVrVa7evzOEgkEtHOzo7m5ubG7rD8ce4dmBb8ngL+x+8pMB74XQX8b5x/T73s3dOtkZlMxvMQTJJOnjzZcQaX13ZvXXSyHTGVSnW8vfuuklJzcmz3hNidO3cc9bG5udl+vnd7JQAAAAAAwLjxNAj76U9/OtLrff7554d+bNS9dLP7bK692xMP02ur44ULF9rPW2eP9bK77qCzxwAAAAAAAMaJ54flj9Lu7YN+tvtMsL3TXf2+vmX3514sFnuusbvGMAxf3N0TAAAAAADgKKYqCNu91c/PLl682H6+ubnp6C6PuwOzg0Ir0zQ73t86jP8w+Xy+/ZxpMAAAAAAAMAk8vWvkoD7//HNH4VBLuVxWNpvt6zVeaoVWxWJRtVpN169f73qGWatOUte6bDbbPnPs+vXrhx6YX6vV2ofkx2KxfWeQAQAAAAAAjKOxCMI+//xzZTIZFYvFgbYKSpJt2woEAi53Njz5fF5nzpxRrVbT2tqa5ufnDwykyuWyksmkpGZotbq6euiahmEon88rmUzKsiytra0dWL+0tCSpeafI3ZNhAAAAAAAA48z3WyOvXbum+fl55XI5lUol2bY90GPchMNh3b9/vx1+LS8vKx6PK5fLybIsFYtFpdNpzc/Pq1arKZPJaGNjo+e6iURCGxsbCofDSqfT7VCsVqupWCwqGo3KsiyZpqn79+933G0SAAAAAABgnPl6IuzXv/51e6tfIBBQIBAYy1BrUOFwWNlsVsvLy8pmsyoWi+1D7w3DkGmaWl1d1bVr1/oKrGKxmKrVqtbW1nTr1i0tLS2pVqspHA5rcXFR+Xz+0G2T02xlZUX1el2hUMjrVgAcgt9TwP/4PQXGA7+rgP/xezqYgO3jZGlxcVGWZbUDsFb407or4unTp3uu8cUXX6hWq+n27dt6+PChHj9+POy20UUkEtHOzo7m5ua0vb3tdTsAAAAAAGDEvMwGfD0R1grBJGljY6N9dtUgVldX9eyzz7rVGgAAAAAAAMaMr88Ia233W11dPVIIJjW3Ep45c8aFrgAAAAAAADCOfB2EtbZAnjt3zpX1stmsK+sAAAAAAABg/Pg6CGtNgT148MDV9QAAAAAAADB9fB2E/exnP5Nt27Isy5X1fvnLX7qyDgAAAAAAAMaPr4OwkydP6uc//7lu3bqlf/iHfzjyemyNBAAAAAAAmF4B27Ztr5voJR6P69ixY/rtb3878Br379/XD37wAz1
"text/plain": [
"<Figure size 1200x900 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"bins = np.linspace(-200, 200, 50)\n",
"sns.regplot(\n",
" x=ak.to_numpy(array[\"ideal_state_770_y\"]),\n",
" y=ak.to_numpy(array[\"yDiffOut\"]),\n",
" x_bins=bins,\n",
" fit_reg=None,\n",
" x_estimator=np.mean,\n",
")\n",
"plt.xlabel(\"y (VELO)\")\n",
"plt.ylabel(\"$y_{corr}$ [mm]\")\n",
"mplhep.lhcb.text(\"Simulation\")\n",
"plt.show()"
]
},
{
"cell_type": "code",
2024-02-29 15:54:19 +01:00
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
2024-02-29 15:54:19 +01:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAABMIAAAOWCAYAAAANzz7PAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB5o0lEQVR4nOz9T2xk2X0n+P6CVRZKkMyMzGo0QJDAOIOuB2k1zmBWA0+76Qy6Z+GC7S4yq2FBpUU7SbWf3iw4EKmUdg+QUsyGuen33CIlb8ooYLLIctuQFxoxpNk9oFHJKM3mQRiDUWUgiQAarWRUWIILtiriLdKkSCYZESRvxL3B+/kAAfHP4bk/BnmTFV/9zjmFTqfTCQAAAAC44sbSLgAAAAAAhkEQBgAAAEAuCMIAAAAAyAVBGAAAAAC5IAgDAAAAIBcEYQAAAADkgiAMAAAAgFwQhAEAAACQC4IwAAAAAHJBEAYAAABALgjCAAAAAMgFQRgAAAAAuSAIAwAAACAXBGEAAAAA5IIgDAAAAIBcEIQBAAAAkAuCMAAAAAByQRAGAAAAQC4IwgAAAADIBUEYAAAAALkgCAMAAAAgFwRhAAAAAOSCIAwAAACAXBCEAQAAAJALgjAAAAAAckEQBgAAAEAuCMIAAAAAyAVBGAAAAAC58GLaBZAvn/nMZ+Ljjz+OF154If7lv/yXaZcDAAAADNl/+2//LT755JN46aWX4pe//OVQr13odDqdoV6RXBsbGwu/cgAAAEChUIh2uz3Ua+oIY6iOhmCTk5MpVnJ+jUYj2u12jI2NxcTERNrlnMuo1j6qdUeoPS1qT4fa0zGqtY9q3RFqT4va06H2dIxq7aNad4Ta07K3txcRkUqjjCCMoRobGzu8UZ88eZJ2OecyNTUVe3t7MTExofYhGdW6I9SeFrWnQ+3pGNXaR7XuCLWnRe3pUHs6RrX2Ua07Qu1peeGFFw6zgWFLNQh75ZVX0rz8MYVCIf6v/+v/SrsMAAAAAAYk1SBsd3c3CoVCqntGHVy/UCikVgMAAAAAgzf8HrSMsXE7AAAAQD6kvkdYp9OJubm5KJVKqVy/2WzGxsZGKtcGAAAAYHhSD8I2Njbij//4j1OtYW5uLv7Nv/k3qdYAAAAAwGClvjSyUqmkXUK8+uqraZcAAAAAwICl3hF248aNtEsgBe12O6amprqOWVpaiqWlpSFVBAAAAFzW2tparK2tdR3TbreHVM3zUg3CnNSYP7/5m78ZH330UURE7O3tdR3barWGUVIuLC0tRavVivHx8bRLOZdRrTtC7ZzfKD/vak/HKNc+qkb5OVd7Oka59lE2ys/7qNY+qnVHjHbtWdVqtXq+3j/wm7/5mwOu5nmFTorHJo6NjUWz2Uz9F+6jjz6K69evp5pI5sXU1FTs7e3F2NhYTExMdB2btY6wg9onJyfjyZMnaZcDnMJ9CtnnPoXR4F6F7MvqfdpPR1ij0Yh2u51K7al2hK2urqYegkVEXLt2LVZXV9MuI1cmJiYydaMCAAAAl9dPU8tBiJeGVDfL/9rXvjbU63344Ydnfm7YtQAAAAAwXKmfGjlMi4uLaZcAAAAAQEpyFYQ9fvw47RIAAAAASEmqe4Rd1IcffhjNZrPv8fV6PdbX18/1NQAAAABcLSMRhH344Yexuroa1Wo16vX6hebodDpRKBQSrgwAAACAUZH5IOz+/fvx8OHDiHgWZgEAAADARWQ6CHv33XdjdXU1IiIKhUIUCgVhGKlZWlqKVqsV4+PjaZcCnMF9CtnnPoXR4F6F7HOfXkyhk+Fk6fbt21Gr1Q4DsFKpFOVyOUqlUkREvPzyyz3n+PnPfx7NZjPeeeed+Oijj+KTTz4ZdNl0MTU1FXt7ezE5ORlPnjxJuxwAAABgyNLMBjLdEXYQgkVEbG9vx507dy481/LycrzyyitJlQYAAADAiBlLu4BuisViRDwLsS4TgkVElEqluHnzZgJVAQAAADCKMh2EHSyBfPXVVxOZb319PZF5AAAAABg9mQ7CDrrAnj59muh8AAAAAORPpoOwb3zjG9HpdKJWqyUy3/e///1E5gEAAABg9GQ6CLt27Vp85zvfiUePHsXf//3fX3o+SyMBAAAA8ivTQVjEs43yy+VyzM/PX2qeDz74ILHOMgAAAABGz4tpF9CP7e3tuH37drzyyiuxsrISN27c6Ovrnj59Gs1mM3Z3d+Odd94ZcJUAAAAAZNlIBGE/+clPIiJid3c3FhcXLzRHp9OJQqGQZFkAAAAAjJDMB2H/4T/8h9jY2IiIiEKhEJ1O59xzCMC4iE4not3j122sEOHXCwAAAEZDpoOw733ve4cb3F80BIuIC38d+dbuRLz4ze5jfvWtiBcEYQAAADASMh2EnQzBKpVKzM7ORrFYPPc+Yd/97nfjww8/HGC1AAAAAGRZpoOwWq0WhUIhisViPH78OG7evHnhue7duxcvv/xygtUBAAAAMErG0i6gm2KxGBER9+/fv1QIdjDXrVu3EqgKAAAAgFGU6SBsZmYmIiJKpVIi833ve99LZB4AAAAARk+mg7DFxcXodDpRr9cTmU9HGEl7+/20K+A8ms1m2iUAAACQokwHYXNzc3Hr1q149OhRIvP95V/+ZSLzkA/9hFxf3hKGZV2tVouVlZW4fv163Lt3L+1yAAAASFGmN8uPiNjc3Izf/u3fjv/yX/5L/OEf/uGl5nrw4EH823/7bxOqjMtoNBoxNTXVdczS0lIsLS0NqaLjPvo44k/+ur+xX/mriNc+HzH+0kBLupJqtVo8evQoarVa1Ov1w1Nei8VilEqluHHjRpTL5ZidnY1KpfLc18/Pz0e9Xo+dnZ3DjzWbzahWq/Ho0aOoVqu6wAAAAIZobW0t1tbWuo5pNBpDquZ5mQ/CSqVSfPe7340//uM/vlQQ9tFHH0WtVkuwMi6j3W7H3t5e1zGtVmtI1TzvL2oRv/ynPgb+dCN+8cPFuPb/Ot/85XL5WHhzUYVC4dxfUywWY39//9TPzc/Px9bW1rnn3N7ePjWoOk2z2YwHDx7Ew4cPu445uF+r1erh2HK5HJVKJaanp2N7ezu2traiXC4f+9p79+5FvV53vwMAAKSg1Wr1fL2fpkwHYQdLGf/Fv/gXcf369XjllVdibm7u3PM0m8145513ki6PSxgbG4uJiYmuY8bHx4dUzfN+8LM+B/7OQsTn7sb/fbwek/+/B11DpEqlEisrK3H79u3DE1Ev62APvWq1GisrK2d2P5VKpVhZWYlKpdL18InNzc1oNpvx+PHjWF1djWq1eubY5eXleOONN6JUKvX9/aysrJwagM3NzR3OdTBfvV6Per0em5ub8c477xyGYycDrpPf8+bm5uHb09PTie0xCAAAQG/j4+MxOTnZdUyj0Yh2uz2kio4rdDqdTipX7sPt27fj/fd/vQFTp9O5UAfM0a/95JNPkiqPC5iamoq9vb2YnJyMJ0+epF3Omf7V/yfivXOU96+mIv7r/+PZAQ8bGxvPfX5ubu5YQDMI9Xo9pqenT/3c/v7+hcK3mZmZUzur1tfXY2Fh4Vy1zc7OPhdKLSwsxPr6el9zPHz4MFZWVp77eKlUit3d3VO/5uTPYxg/BwAAALpLMxvI9Gb5d+/ejU6nEwdZ3UVDMDiv65++2PjFxcVTP3/Wx5NUKpWeWyYY8Sz8uWgH2ll1nycEq1arz3VmFYvF2NnZ6TsEi3jWgba7u/vc9/L06dMzvyapzjsAAACuhkwHYQcvwguFwmEgdtEHnMdrnzvf+N/75/FnLTvsthwxSadd5zLXPu1rzxMu1Wq1mJ2dfe7rf/zjH58a2vVTzwcffHCsBpvhAwAA0K9M7xF27dq1KJfL8f7778fq6mqUy+W4cePGueep1+vx7W9/O376058mXyRX0pfKEV//YX8b5n/2UxFv/nOmc1ZIdJHf21FXr9fjzp07z318c3PzQiHYgYNusrOWgY6yg83/hxWcjppqtdr3oQwAAACnyXQQFhHxxhtvxPT0dHzta1+78By3bt2K119/PZdhBBdz7aWIP/v9iC/3cYDid/8gYvylgZfUl9N+x19++eV
"text/plain": [
"<Figure size 1200x900 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"bins = np.linspace(-200, 200, 50)\n",
"sns.regplot(\n",
" x=ak.to_numpy(array[\"ideal_state_9410_y\"]),\n",
" y=ak.to_numpy(array[\"yDiffOut\"]),\n",
" x_bins=bins,\n",
" fit_reg=None,\n",
" x_estimator=np.mean,\n",
")\n",
"plt.xlabel(\"y (T)\")\n",
"plt.ylabel(\"$y_{corr}$ [mm]\")\n",
"mplhep.lhcb.text(\"Simulation\")\n",
"plt.show()"
]
},
{
"cell_type": "code",
2024-02-29 15:54:19 +01:00
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
2024-02-29 15:54:19 +01:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAABMIAAAOWCAYAAAANzz7PAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB27ElEQVR4nOz9T2yb+Z0n+H8oO0EFyciUqzGAIB3aVNcijT1sm3LNYW47pnbmkGwyU5Td6EySw7SlJBjMAsKEaieHnT1UVaiZ6DJAuqTqPmxlc3BRSSdI9tARq+Y8sMTK3oKZiFUN2BAwQFksdYIUkirxd/BPbEnWH0oi+ZB8Xi+AgGh9+Dwfil/T4tvf7/fJNJvNZgAAAADAkBtJugEAAAAA6AVBGAAAAACpIAgDAAAAIBUEYQAAAACkgiAMAAAAgFQQhAEAAACQCoIwAAAAAFJBEAYAAABAKgjCAAAAAEgFQRgAAAAAqSAIAwAAACAVBGEAAAAApIIgDAAAAIBUEIQBAAAAkAqCMAAAAABSQRAGAAAAQCoIwgAAAABIBUEYAAAAAKkgCAMAAAAgFQRhAAAAAKSCIAwAAACAVBCEAQAAAJAKgjAAAAAAUkEQBgAAAEAqCMIAAAAASAVBGAAAAACpIAgDAAAAIBUEYQAAAACkwtWkGyBdPv3pT8eHH34YV65ciX/8j/9x0u0AAAAAPfY//sf/iI8//jiee+65+M1vftPTc2eazWazp2ck1a5cuRJ7e3tJtwEAAAAkbGRkJD7++OOentOMMHpqPwgbGRmJ8fHxrpxje3u76+dw7uTPm9Zze87peM5JnjuNzznJc3vO6XjOSZ47jc85yXN7zul4zkmeO43POclze87dO+/+ea5cudK1c5yoCT00MTHRjIjmxMTEQJ/DuZM/b1rP7Tn3VhrPncbnnOS5PefeSuO50/ickzy359xbaTx3Gp9zkuf2nAf/PMdJdEbYCy+8kOTpD8lkMvHf/tt/S7oNAAAAALok0SBsa2srMplMNBPcpmz//JlMJrEeAAAAAOi+kaQbSFqSIRwAAAAAvZP4ZvnNZjOKxWLkcrlEzt9oNGJ1dTWRcwMAAADQO4kHYaurq/Hnf/7nifZQLBbjn//zf55oDwAAAAB0V+JLIwuFQtItxIsvvph0CwAAAAB0WeJB2PXr15NuAQAAAIAUSHRppCs10g0LCwuxu7sbo6OjSbfSU0k97yR/3mk9d1KMsfScOylp/Hmn8TknKY0/7zQ+5ySl8eedxuecpLT+vNM4vpOShuecaSZ42cSRkZFoNBqJ/4A/+OCDGBsbi729vUT7SIPJycl4/PhxTExMxKNHj5Jup+OG/fmRPGOMbjPG6DZjjF4wzug2Y4xuG/YxluTzS3RpZLlcTjwEi4i4du1alMvlpNsAAAAAoIsSDcK++c1v9vR877333onf63UvAAAAAPRW4pvl99L8/HzSLQAAAACQkFQFYRsbG0m3AAAAAEBCEr1q5EW999570Wg02q6v1+uxsrJyrscAAAAAMFwGIgh77733olwuR7VajXq9fqFjNJvNyGQyHe4MAAAAgEHR90sj79+/H1NTU7G6uhpbW1vRbDYvdBt2U1NTkclkYm1t7dyPXVpaipmZmRgbG4tMJhNTU1MxOzsb1Wq1C50CAAAAJKOvZ4T98Ic/jHK5HBERmUwmMplMKkKt81pcXLzQTLlqtRqzs7PRaDSiUChEpVKJXC4XtVotFhcXY2ZmpvXn2Wy2840PoYWFhdjd3Y3R0dGkW2FIGWN0mzFGtxlj9IJxRrcZY3SbMdY9mWYfJ0u3bt2KWq3WCsByuVzk8/nI5XIREfH888+feYz3338/Go1GvPnmm/HBBx/Exx9/3O22e6parcbMzEzrfqVSiWKxeK7Hzc3NxcrKyjM109PTUavVIpfLxebmZkfCsMnJyXj8+HFMTEzEo0ePLn08AAAAYLAkmQ309Yyw/RAsImJ9fT1u37594WOVSqV44YUXOtVaX2g0GjE7O3upx+VyuWNDsIinodrU1FTU6/WYnZ2N9fX1S/ULAAAAkKS+3iNsfwZSqVS6VAgW8TTwuXHjRge66h/37t2L69evn3um1v5yyIinyypPksvlWrPLqtVqrK6uXrRVAAAAgMT1dRC2vwTyxRdf7MjxTpr5NIhWV1djbW0tKpXKuR5Xr9cPbYJ/586dU+vv3r3b+vq00AwAAACg3/V1ELY/C+zJkycdPd6gq9frMT8/H6VSKfL5/Lkeu3/xgYiIQqFw5myyg/uNNRqNC12VEgAAAKAf9HUQ9q1vfSuazWbUarWOHO+v/uqvOnKcpM3OzkY+nz8UarXr4PLGdkO0/Zl5EREPHjw49zkBAAAA+kFfB2HXrl2L73znO/HgwYP4+7//+0sfbxiWRi4uLkatVjv3ksiIeCZQbHfJ6cHAzIwwAAAAYFD1dRAWEa3lfxe5OuJB7777bsdmliWlVqvF0tJSrKysHJql1a6De4NFRNvHOFo36D9HAAAAIJ2uJt1AO9bX1+PWrVvxwgsvxOLiYly/fr2txz158iQajUZsbW3Fm2++2eUuu+/27dtRLBZjbm7uQo9/+PDhofvtXm3y+eefP3R/Y2Pj3HuTAQAAACRtIIKwt99+OyIitra2Yn5+/kLHaDabkclkOtlWT+3PiHv99dcvfIx6vX7o/kVnhG1tbV24BwAAAICk9H0Q9vWvf721wXsmk4lms3nuYwxyABbxdF+utbW1WF9fb3sW13GOBmEX1Wg0OnIcAAAAoD3NZsTeGZHISCZiwCOQruvrIOz1119vbXB/0RAsIi78uH7QaDRidnY25ubmolAoXPpYF3E0fHvy5Mml+gAAAADOZ68ZcfXbp9d89HLEFUHYqfo6CDsaghUKhZiZmYlsNnvufcJee+21eO+997rYbXfcvn07crlcX13xshMzwra3t2NycvLSx1lYWIiFhYVLHwcAAAA43fLyciwvL1/6ONvb2x3o5mL6Ogir1WqRyWQim83GxsZG3Lhx48LHunfv3jObvve7paWlqNVqsbm52ZHjZbPZjoRYl1meuW9vby8eP3586ePs7u5e+hgAAADA2XZ3dzvyWT5JfR2EZbPZ+OCDD+L+/fuXCsH2j3Xz5s0OddZ9tVotFhcXo1wud+wKjdevX+9IENbubLzTjIyMxPj4+KWPMzo6euljAAAAAGcbHR2NiYmJSx9ne3s79vb2OtDR+fV1EDY9PR1vv/1221c3PMtlrrjYa7Ozs5HP56NUKnXsmBedyXU0POvEjLDx8fF49OjRpY8DAAAA9EantieanJxMbGZZXwdh8/Pz8dZbb3XsaoeDMiNsaWkp6vV6FAqFmJ2dPbP+YFD16quvxoMHD1r37969G8ViMSIibt26FbVa7dDj2gm1jm6OPzU1deZjoB+1O+YBAAAYTn0dhBWLxbh582Y8ePAg/v2///eXPt6PfvSj+Ff/6l91oLPuev/99yMiolqtnvuxtVrtUNiVy+VaQdj09PSh2nq93tayy62trUP3L3v1SuilWq0WDx48iNXV1SgUClGpVJJuCQAAgIT0dRAWEVGpVOKP/uiP4m/+5m/iX/7Lf3mpY7366qsDEYR1y61btw7dbzcIOzjjLJvNdmypKuyHVLVaLer1eusqr/vj7Pr165HP52NmZubYAHZ2djbq9fqhC0o0Go2oVqvx4MGDqFarHdkXDwAAgOEwknQDZ8nlcvHaa6/Fn//5n1/qOB988MGhmVL9rFwuR7PZbPt2MJiqVCqHvlcul1vfy+fzh5aFPXz4sK1+NjY2Wl8fDdPSbnV1NTKZzLlvR2fnXdRFzj02Nnbi8WZnZy90zPPMXmw0GrG4uNj6OSwtLUW1Wo16vd4KrRqNRtRqtahWq7G0tBQzMzOt+sXFxVhdXY3Z2dlYW1t75vj37t2LV199NdbW1oRgAAAAHNLXM8J+9KMfRUTEH/zBH8TY2Fi88MILrWV+59FoNOLNN9/sdHsD6c6dO7G6uhoR0XYweLBucXGxK30Nqrm5ubhz507U6/VW+HKSQqEQi4uLcevWrY7tU9VsNqNer0e1Wo3FxcU
"text/plain": [
"<Figure size 1200x900 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"bins = np.linspace(-1, 1, 50)\n",
"sns.regplot(\n",
" x=ak.to_numpy(array[\"ideal_state_9410_ty\"]),\n",
" y=ak.to_numpy(array[\"yDiffOut\"]),\n",
" x_bins=bins,\n",
" fit_reg=None,\n",
" x_estimator=np.mean,\n",
")\n",
"plt.xlabel(\"ty (T)\")\n",
"plt.ylabel(\"$y_{corr}$ [mm]\")\n",
"mplhep.lhcb.text(\"Simulation\")\n",
"plt.show()"
]
},
{
"cell_type": "code",
2024-02-29 15:54:19 +01:00
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
2024-02-29 15:54:19 +01:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAABMEAAAOVCAYAAABgbPdiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAACTd0lEQVR4nOz9X2xj55nn+/7IqhjlHTeLUg0GIMiL1FJ74GBf9IhUzUX21bTI9ADb3nHaZLkRt+2LjignyO4GOBPKFV/bCtXTBM7pOYnJSg7Qdgx0mXTGgT1Axlr2XO0BBiWxPHfGJGQ5gAgCg05pFds5NhwX17ngiK0/FLkocpGLS98PQESUHr7r0Z9KN39532cFbNu2BQAAAAAAAPhYcNYNAAAAAAAAAG4jBAMAAAAAAIDvEYIBAAAAAADA9wjBAAAAAAAA4HuEYAAAAAAAAPA9QjAAAAAAAAD4HiEYAAAAAAAAfI8QDAAAAAAAAL5HCAYAAAAAAADfIwQDAAAAAACA7xGCAQAAAAAAwPcIwQAAAAAAAOB7hGAAAAAAAADwPUIwAAAAAAAA+B4hGAAAAAAAAHyPEAwAAAAAAAC+RwgGAAAAAAAA3yMEAwAAAAAAgO8RggEAAAAAAMD3CMEAAAAAAADge4RgAAAAAAAA8D1CMAAAAAAAAPgeIRgAAAAAAAB8jxAMAAAAAAAAvkcIBgAAAAAAAN8jBAMAAAAAAIDvEYIBAAAAAADA9wjBAAAAAAAA4HuEYAAAAAAAAPA9QrA5UavVtL6+rqWlJQUCAQUCAS0tLWljY0OWZZ153a2tLaVSKS0sLPTWzGQyMk3TU2sCAAAAAACMgxDM4yzLUiaTUSKRULlcVqPR6H2t0Whoa2tLCwsLKpfLI61rmqYWFha0sbEhSapUKqrX6yoUCqrVakqlUkqlUiMFbG6sCQAAAAAAMAkB27btWTeB/izLUiKROBJ8DZLNZlUqlYbWmaapVCo18DWJREK1Wk2GYWh3d1fhcHjqawIAAAAAAEwKIZiHpVIpmaapeDyuGzduKB6PS+oejbx9+7a2trZOvKZSqSidTp+6pmVZunr1qizLkmEYqtfrfesajYaWlpYkSclkUtvb21NdcxRf/vKX9dlnn+nChQv65//8n09kTQAAAAAAIP3P//k/9eDBA126dEm/+93vZt3OeGx4UqlUsiXZ+Xz+1Jp6vW7H43FbUu8RDocHrptMJnu1pVJpYG06nXZU68aaowgGg0d+Bjx48ODBgwcPHjx48ODBgwePyT6CweBE3sPPEjvBPGppaUmGYQzdLXV4d9WB7e1tJZPJobX7+/sDjyRWq1VlMhlJUjgc1v7+/lTWHNVDDz2k3//+9woGg4pEIpKkVqulTqdz5HPTNOvr04M3ru+FHmZ9fXrwxvXpwRvXpwdvXJ8evHF9evDG9enBG9f3Qg+zvj49eOP6p/Vw8LkvfelL+vzzz2fS18TMOoXDSbu7u7Yke39/31F9oVA4ks4WCoW+ddlstleTTCYdrX143UqlMpU1RxWNRm1JdjQaHfi5aZr19enBG9f3Qg+zvj49eOP69OCN69ODN65PD964Pj144/r04I3re6GHWV+fHrxx/dN68EJfk8LdIT3o1q1bymazjgfHH9/19dvf/rZv3eE7SB7MFxvGMIwjfU1jTQAAAAAAgEm7OMuLP/roo7O8/BGBQED/43/8j1m3IUl6+umnjwRFwxwPn44fj5S6w/QPu3btmuO1D+5OWa1WXV8TAAAAAADADTMNwer1ugKBgOwZjiU7uH4gEJhZD8c53VF1wLKsI8/7BWimaQ6t6ed4Xa1W6/XnxpoAAAAAAABuOPfHIWcZwE3Kwa6qA/2G4t++ffvIc6dHLa9cuXLk+c7OjqtrAgAAAAAAuGGmO8GkbgiVTqdHOv43SZZlHZlrNY8Oh0jZbLZvzfGg7Ky7tur1uqtrAgAAAAAAuGHmIVi5XNa3v/3tmfaQTqf1J3/yJzPtYRylUqn38cbGRt+a44HVWR0+eunGmuNotVqKxWK9j49/zqlcLqdcLjeRngAAAAAAmIVisahisTjSa/q9lz74nB/MPATrd3Rv2pwOdPeiRqPRG1BfKBRO3Y111qDp+BHHe/fuubrmODqdjprN5tDPDdNutyfSDwAAAAAAs9Jut0d+P3zgLO+l58HMQ7DFxcVZtzDXCoWCpO4Rw3w+7/r1JrVry401g8GgIpGIpG5S3el0jnzOqVAoNHYvuVxO7XZ7ImvRw/xe3wu88DOgh9lf3ytm/XOY9fXpwTu88DOYdQ+zvr5Xepg1L/wM6GH21/cCL/wM6GHy1w+FQopGoyO9pt976YPP+UHAnuFk+AsXLmh/f3/m/2Vz//59LS4u6sGDBzPtY1S1Wk2JRELhcFi7u7sDZ3ItLCwcCZuc/tpN01Qqleo9T6fTqlQqrq15FrFYTM1mU9FoVHt7e6d+DpgF/hbhFfwtwiv4W4SX8PcIr+BvEV7h9/fXM707pJfuzOilXpxaW1uTJL3//vtDh9JPasfd4XXcWBMAAAAAAMANMw3BCoXCzHeBSdLly5d7xwrnxfr6umq1miqViuLx+ND643O4nDp+VPHwOm6sCQAAAAAA4IaZhmDf//73p3q9jz/++NSvTbuXcZTLZZXLZZVKJaXTaUevWVlZOfLc6Ryu40Prl5aWXF0TAAAAAADADTMNwaZtfX191i2MzTRNra+vq1QqKZvNOn5dIpE48rzRaDh6Xb1eP/L88N083VgTAAAAAADADecqBNvZ2Zl1C2Op1WpKpVIqFAojBWDSyV1bTgOrw7u7wuHwkdljbqwJAAAAAADghouzbuAsPv74Y8dH76RuOFMqlUZ6jdc0Gg2trq4qn88rn8+P/Pp4PK5wONz7Gdy+fdvRUcrDweHx0MuNNQEAAAAAANwwFyHYxx9/rEKhINM0He82Os62bQUCgQl3Nh2NRkOJRELZbNbxAP9Go6FqtXokMLt+/brK5bKk7q4yJw7XbWxsnPi6G2sCAAAAAABMmudDsBs3bmhra0tSN8g6byzLUiqV0vXr10e6g2Umk9HNmzePfG59fb0XWJmmOXSNwzWGYfSd3eXGmpOQy+XUbrc9cfdRnG/8LcIr+FuEV/C3CC/h7xFewd8ivMLvf4sB28PJ0ltvvaVMJiNJvV1c47QbCAT04MGDifQ2DZZlKZFIyDAMlUolR69pNBq93VW7u7snvp5KpXpBVKVSGXh88XDANWgQvxtrjiIWi6nZbCoajWpvb2/s9QAAAAAAQJef3nN7OgRbWVlRrVZTIBCQbdsyDEPxeLw3SP3KlStD1/jtb38ry7L05ptv6v79+3MVgiUSCcdHDI87LWBqNBpaWlqS1J3p1S8ok7oB3MLCgqTu3Ru3t7dPvZYba47CT/8gAQAAAADwEj+95/Z0CBYMBns7wN577z2trq6eea1Go6FHH310bkKwcQIwafCOuWq12tthVygU+g7aP7h+OBzW3bt3FQ6HB17PjTWd8tM/SAAAAAAAvMRP77mDs25gkIOQJJ/PjxWASd35U1evXp1AV+7LZDJjBWDDjhim02ltb28rHA5rY2Ojdz3LsmSaZi+sisfjjsMqN9YEAAAAAACYFE8PxjcMQ3fu3NG1a9cmsp7TuVqzVqlUXL9GMpnU/v6+tra2dOvWLa2ursqyLIXDYa2srAyd7TWtNQEAAAAAACbB08chNzY29O///b9XqVTSt7/97Vm3A4/y09ZMAAAAAAC8xE/vuT19HPIHP/iBbNse62jgYT/5yU8msg4AAAAAAADmi6d3gknS1taWCoWCPv74Y/3BH/zBWGtdu3ZNt2/fnlBn8IqDVDoYDCoSiQyszeVyyuVyU+oMAAAAAADvKhaLKhaLA2tarZY6nY4vdoJ5eiaY1B2Kv729rUwmo1/+8pdnXufu3bsT21EGb+p0Omo2mwNr2u32lLoBAAAAAMDb2u320PfRfuL5EEyStre3tbKyokcffVQbGxtaXFx09Lp79+7JsizV63W9+eabLneJWXOyEywUCk2pm6NsW+oM2XMZDEiBwHT6AQAAAAAgFAopGo0OrDn
"text/plain": [
"<Figure size 1200x900 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
2024-02-29 15:54:19 +01:00
"bins = np.linspace(-0.5, 0.5, 50)\n",
"sns.regplot(\n",
" x=ak.to_numpy(array[\"ideal_state_770_ty\"]),\n",
" y=ak.to_numpy(array[\"yDiffOut\"]),\n",
" x_bins=bins,\n",
" fit_reg=None,\n",
" x_estimator=np.mean,\n",
")\n",
2024-02-29 15:54:19 +01:00
"plt.xlabel(\"$t_y$\")\n",
"plt.ylabel(\"$y_{corr}$ [mm]\")\n",
"mplhep.lhcb.text(\"Simulation\")\n",
2024-02-29 15:54:19 +01:00
"plt.show()\n",
"# plt.savefig(\n",
"# \"/work/cetin/LHCb/reco_tuner/parameterisations/plots/bend_y_ty_dist.pdf\",\n",
"# format=\"PDF\",\n",
"# )"
]
},
{
"cell_type": "code",
2024-02-29 15:54:19 +01:00
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
2024-02-29 15:54:19 +01:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAABMIAAAOWCAYAAAANzz7PAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB/EklEQVR4nOz9X2xk550feP+Kagut2FtdbAUBCBLY6eJooWBvpoutXDhXSZPvXKy99sRkazAeyReTJjWGdxZgYtKysZeWTAYhsMi+HpH2vECk8UWLbceGFWBjln0ZYNPN0uRO2CxLGoBEAUHUrK6xYcGWqvaiX9L8WyySVXWq6nw+ACH++dXz/Mg6XeL58jnPyTQajUYAAAAAwIAbSroBAAAAAOgGQRgAAAAAqSAIAwAAACAVBGEAAAAApIIgDAAAAIBUEIQBAAAAkAqCMAAAAABSQRAGAAAAQCoIwgAAAABIBUEYAAAAAKkgCAMAAAAgFQRhAAAAAKSCIAwAAACAVBCEAQAAAJAKgjAAAAAAUkEQBgAAAEAqCMIAAAAASAVBGAAAAACpIAgDAAAAIBUEYQAAAACkgiAMAAAAgFQQhAEAAACQCoIwAAAAAFJBEAYAAABAKgjCAAAAAEgFQRgAAAAAqSAIAwAAACAVBGEAAAAApMKVpBuAy/r0pz8dH330UTz11FPxD/7BP0i6HQAAAKCJ//pf/2t88skncfXq1fjVr37V1bkzjUaj0dUZoc2eeuqpqNfrSbcBAAAAnMPQ0FB88sknXZ3TijD63l4QNjQ0FCMjI20du1KpdGxsc/XnPObqn3nM1T/zmKt/5jFX/8xjrv6Zx1z9M8+gzjWI39OgzjVI39Pe+E899VTbxz5TA/rc6OhoIyIao6OjfTW2ufpzHnP1zzzm6p95zNU/85irf+YxV//MY67+mWdQ5xrE72lQ5xqk76mb38tRia4Ie+6555Kc/pBMJhP/9//9fyfdBgAAAAAdkmgQtrW1FZlMJhoJblO2N38mk0msBwAAAAA6byjpBpKWZAgHAAAAQPckvll+o9GI6enpyOfzicxfrVZjbW0tkbkBAAAA6J7Eg7C1tbX45//8nyfaw/T0dPzhH/5hoj0AAAAA0FmJXxo5OTmZdAvxwgsvJN0CAAAAAB2W+Iqw69evJ90CA6JSqcTY2FjTmvn5+Zifn+9SRwAAAJAuKysrsbKy0rSmUql0qZvjEg3C3KmRdqrX67Gzs9O0plarnWvM+fn5qNVqkc1mL9NaqufqlkH9+Q3qXN0yqD+/bs01iMdEhOeqn+bqJs9V/8zVLYP68xvUubplEH9+g/g9dZvn6rBarXbmuXmSMo0Eb5s4NDQU1Wo18X8Ejx8/juHh4ajX64n2wcWMjY3Fzs5ODA0NxcjISNNaK8K6b+/5GR0dje3t7aTboQc4JjiJ44KTOC44ieOCoxwTnMRxkZxWV4TV6/VEnp9EV4QtLS0lHoJFRFy7di2WlpaSboNLGhkZ8QIHAAAACWplAcpeUJmERDfL//rXv97V+T744INTv9btXgAAAADorsTvGtlNc3NzSbcAAAAAQEJSFYQ9fPgw6RYAAAAASEiie4Rd1AcffBDVarXl+nK5HKurq+d6DAAAAACDpS+CsA8++CCWlpaiWCxGuVy+0BiNRiMymUybOwMAAACgX/R8EPbqq6/G8vJyRDwJswAAAADgIno6CPvhD38YS0tLERGRyWQik8kIw6DPzM/PR61Wi2w2m3Qr9AjHBCdxXHASxwUncVxwlGOCkzguOE2m0cPJ0q1bt6JUKu0HYPl8PgqFQuTz+YiIePbZZ88c48MPP4xqtRpvv/12PH78OD755JNOt02XjY2Nxc7OToyOjsb29nbS7QAAAABNJHke39MrwvZCsIiIjY2NuH379oXHWlhYiOeee65drQEAAADQZ4aSbqCZXC4XEU9CrMuEYBER+Xw+bty40YauAAAAAOhHPR2E7V0C+cILL7RlvNXV1baMAwAAAED/6ekgbG8V2KNHj9o6HgAAAADp09NB2De/+c1oNBpRKpXaMt73v//9towDAAAAQP/p6SDs2rVr8Z3vfCfu3bsXf/d3f3fp8VwaCQAAAJBePR2ERTzZKL9QKMTMzMylxnn//ffbtrIMAAAAgP5zJekGWrGxsRG3bt2K5557LhYXF+P69estPe7Ro0dRrVZja2sr3n777Q53CQAAAEAv64sg7Be/+EVERGxtbcXc3NyFxmg0GpHJZNrZFgAAAEDXNBoR9UbzmqFMhPjjdD0fhP35n/95rK2tRUREJpOJRuOMZ/wEAjAAAACg39UbEVe+1bzm429HPCUGOVVP7xH2ve99L1ZXV/fDr4uEYJd5HAAAAACDo6dXhO3d5XFvJdjk5GRMTU1FLpc79z5hb7zxRnzwwQcd7BYAAACAXtbTQVipVIpMJhO5XC4ePnwYN27cuPBYd+/ejWeffbaN3QEAAADQT3r60shcLhcREa+++uqlQrC9sW7evNmGrgAAAADoRz0dhE1MTERERD6fb8t43/ve99oyDgAAAAD9p6cvjZybm4uf//znUS6X2zKeFWGDrVKpxNjYWNOa+fn5mJ+f71JHDLpqtbq/chUAAICIlZWVWFlZaVpTqVS61M1xPR2ETU9Px82bN+PevXvxL//lv7z0eD/60Y/in/2zf9aGzuhF9Xo9dnZ2mtbUarUudcOgKpVKce/evVhbW4vJyclYX19PuiUAAICeUavVzjw3T1JPB2EREevr6/H7v//78e/+3b+LP/qjP7rUWK+//rogbIANDQ3FyMhI05psNtulbkjCXkhVKpWiXC7v3zU2l8tFPp+P69evR6FQiKmpqZicnDz2+JmZmSiXy7G5ubn/uWq1GsViMe7duxfFYjGq1WoXvyMAAID+ks1mY3R0tGlNpVKJer3epY4OyzQajUYiM5/D2tpavPrqq/Hhhx9eeIzHjx/H9evX45NPPmljZ/SCsbGx2NnZidHR0dje3k66nUtZW1uLubm5cz+uUCgcCm8uKpPJnPsxuVwudnd3T/zazMxM3L9//9xjbmxsnBhUnaRarcbrr78ey8vL556nUCjE5ORkjI+Px8bGRty/f//Yz3IvHCuVSscePz09bUUYAADQNZ/UI658q3nNx9+OeKqnd4RP9jy+p1eE/ehHP4qIiL//9/9+DA8Px3PPPRfT09PnHqdarcbbb7/d7vag7WZnZ+POnTtRLpfj9ddfbxoiTU5OxuLiYty6datt+1Q1Go0ol8tRLBZjcXHx1NVP+Xw+FhcXY3JysunNLNbX16NarcbDhw9jaWkpisXiqbULCwvx4osvRj6fb/n7WVxcPDEAm56e3h9rb7xyuRzlcjnW19fj7bffjmq1GqVS6VjAdfR7Phh0jY+Pt23PQgAAALqvp1eE3bp1K9599939jxuNxoVWrBx8rBVhg2eQVoQdNTc3F2tra8c+342VSOVyOcbHx0/82u7u7oXCt4mJiRNXVq2ursbs7Oy5epuamjoWSs3Ozsbq6mpLYywvL8fi4uKxz+fz+dja2jrxMUefDyvCAACAbrIi7PJ6+kdz586daDQasZfVXTQEg3512mWSF7l88rzy+XwUCoVjn5+enr7wCrTT+j5PCFYsFo+tzMrlcrG5udlyCBbxZAXa1tbWse/l0aNHpz7GHSIBAAD6W08HYXsnzZlMZj8Qu+gb9KPTLjtsdjlip+e/zNwnPfY84VKpVIqpqaljj//5z39+YmjXSj/vv//+oR5shg8AAPSzH7x7dk2a9fQeYdeuXYtCoRDvvvtuLC0tRaFQiOvXr597nHK5HK+99lr8zd/8TfubhA46LSS6yL+Dflcul+P27dvHPr++vn6hEGzP3mqy0y4D7Vd7G/93KzTtN8ViseUbMgAAQK9oJeT6yv0nl0Z++Wbn++lHPR2ERUS8+OKLMT4+Hl//+tcvPMbNmzfjS1/6UirDA7iMk/7NPPvss20dr9V/lzMzM8dWa83OzrYlzMjn87GwsHChO092W6lUitXV1SgWi1EulyOXy8X169cjl8tFPp+PF198MSIi7t69G0tLS00vOy0
"text/plain": [
"<Figure size 1200x900 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
2024-02-29 15:54:19 +01:00
"bins = np.linspace(-0.5, 0.5, 100)\n",
"sns.regplot(\n",
" x=ak.to_numpy(array[\"dSlope_yEndT\"]),\n",
" y=ak.to_numpy(array[\"yDiffOut\"]),\n",
" x_bins=bins,\n",
" fit_reg=None,\n",
" x_estimator=np.mean,\n",
")\n",
2024-02-29 15:54:19 +01:00
"plt.xlabel(\"$\\Delta t_y$\")\n",
"plt.ylabel(\"$y_{corr}$ [mm]\")\n",
"mplhep.lhcb.text(\"Simulation\")\n",
2024-02-29 15:54:19 +01:00
"plt.show()\n",
"# plt.savefig(\n",
"# \"/work/cetin/LHCb/reco_tuner/parameterisations/plots/bend_y_deltaty_dist.pdf\",\n",
"# format=\"PDF\",\n",
"# )"
]
},
{
"cell_type": "code",
2024-02-29 15:54:19 +01:00
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
2024-02-29 15:54:19 +01:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAABMIAAAOWCAYAAAANzz7PAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB8sUlEQVR4nOz9X2yj6X0n+P6oqgTtdcKiqoMBtBKwLioNODgXs0VV70X2aqeonb2IT5JpqXsQr52LSUlOEMw50MKSy8ZeHXSXVYsI2LMHmZbavjgODGy35EyC5ACJRXuv9mJPS2zvnTGJ2B2gNAQGcInNSU6MpJs8FzWSpSqJoqRXfF/y/XwAokTp0fP+JPERi189fwrdbrcbAAAAADDixtIuAAAAAAAGQRAGAAAAQC4IwgAAAADIBUEYAAAAALkgCAMAAAAgFwRhAAAAAOSCIAwAAACAXBCEAQAAAJALgjAAAAAAckEQBgAAAEAuCMIAAAAAyAVBGAAAAAC5IAgDAAAAIBcEYQAAAADkgiAMAAAAgFwQhAEAAACQC4IwAAAAAHJBEAYAAABALgjCAAAAAMgFQRgAAAAAuSAIAwAAACAXBGEAAAAA5IIgDAAAAIBcEIQBAAAAkAuCMAAAAAByQRAGAAAAQC4IwgAAAADIBUEYAAAAALlwM+0CyJfPfvaz8dOf/jRu3LgR/+Sf/JO0ywEAAAAG7D/8h/8Qn376abz00kvxd3/3dwO9dqHb7XYHekVybWxsLDzkAAAAgEKhEJ1OZ6DXNCOMgToegk1OTqZYycU1m83odDoxNjYWExMTaZdzIcNa+7DWHaH2tKg9HWpPx7DWPqx1R6g9LWpPh9rTMay1D2vdEWpPy/7+fkREKhNlBGEM1NjY2NFAffLkSdrlXMjU1FTs7+/HxMSE2gdkWOuOUHta1J4OtadjWGsf1roj1J4WtadD7ekY1tqHte4Itaflxo0bR9nAoKUahL3yyitpXv6EQqEQ/+7f/bu0ywAAAADgmqQahO3t7UWhUEh1z6jD6xcKhdRqAAAAAOD6DX4OWsbYuB0AAAAgH1LfI6zb7cbc3FyUy+VUrt9qtWJjYyOVawMAAAAwOKkHYRsbG/E7v/M7qdYwNzcX//yf//NUawAAAADgeqW+NLJaraZdQrz66qtplwAAAADANUs9CLt9+3baJQAAAACQA6kujXRSY/784i/+Ynz88cfxi7/4i2mXkitLS0vRbrejWCymXcqFDGvdEWrn4ob5+672dAxz7cNqmL/nak/HMNc+zIb5+z6stQ9r3RHDXfswSzMbKHRTPDZxbGwsWq1W6g+4jz/+OMbHx6PT6aRaRx5MTU3F/v5+TE5OxpMnT9Iu50KGuXbIC+MUss84heFgrEL2DfM4TbP2VGeEra6uph6CRUTcunUrVldX0y4jV5rNZkxNTfVss7S0FEtLSwOqCAAAALiqtbW1WFtb69mm2WwOqJoXpRqEffWrXx3o9T766KP43Oc+l4la8q7T6cT+/n7PNu12e0DVAAAAAElot9vnvt5PU6pB2KAtLi7GX/7lX6ZdBvFsWezExETPNlmYLQgAAAD0r1gsxuTkZM82zWYzte2pchWE7ezspF0C/8nExMTQrWEGAAAAeutnm6PDPcLSMJRB2EcffRStVqvv9o1GI9bX1y/0OQAAAACMlqEIwj766KNYXV2NWq0WjUbjUn10u90oFAoJVwYAAADAsMh8EPbw4cN4/PhxRDwLswAAAADgMjIdhH3ve9+L1dXViIgoFApRKBSEYaRmaWkp2u22Tfwhw4xTyD7jFIaDsQrZZ5xeTqGb4WTp3r17Ua/XjwKwcrkclUolyuVyRES8/PLL5/bxk5/8JFqtVrz33nvx8ccfx6effnrdZdPD4YZ4k5OTNssHAACAHEozG8j0jLDDECwiYnt7O+7fv3/pvpaXl+OVV15JqjQAAAAAhsxY2gX0UiqVIuJZiHWVECwiolwux507dxKoCgAAAIBhlOkg7HAJ5KuvvppIf+vr64n0AwAAAMDwyXQQdjgL7OnTp4n2BwAAAED+ZDoI+/rXvx7dbjfq9Xoi/X3rW99KpB8AAAAAhk+mg7Bbt27FN7/5zXj33XfjP/7H/3jl/iyNBAAAAMivTAdhEc82yq9UKjE/P3+lfj788MPEZpYBAAAAMHxupl1AP7a3t+PevXvxyiuvxMrKSty+fbuvz3v69Gm0Wq3Y29uL995775qrBAAAACDLhiII++EPfxgREXt7e7G4uHipPrrdbhQKhSTLAgAAAGCIZD4I+93f/d3Y2NiIiIhCoRDdbvfCfQjAAAAAgEHrdiM658QYY4UIscXgZDoIe+edd442uL9sCBYRl/48AAAAgMvqdCNufqN3m0/ejLghCBuYTAdhz4dg1Wo1Zmdno1QqXXifsLfffjs++uija6wWAAAAgCzLdBBWr9ejUChEqVSKnZ2duHPnzqX7evDgQbz88ssJVgcAAADAMBlLu4BeSqVSREQ8fPjwSiHYYV93795NoCoAAAAAhlGmg7CZmZmIiCiXy4n098477yTSDwAAAADDJ9NB2OLiYnS73Wg0Gon0Z0YY5Fur1Uq7BAAAAFKU6SBsbm4u7t69G++++24i/f3xH/9xIv1wdc1mM6ampnre1tbW0i6TEVCv12NlZSXGx8fjwYMHaZcDAAAw0tbW1s59vd9sNlOrL9Ob5UdEbG5uxi//8i/Hv/23/zZ+8zd/80p9PXr0KP7Fv/gXCVXGVXQ6ndjf3+/Zpt1uD6ga0lKv1+Pdd9+Ner0ejUbj6JTXUqkU5XI5bt++HZVKJWZnZ6Narb7w+fPz89FoNGJ3d/fofa1WK2q1Wrz77rtRq9XMAgMAABigdrt97uv9NGU+CCuXy/H222/H7/zO71wpCPv444+jXq8nWBlXMTY2FhMTEz3bFIvFAVVzeRsbG7G4uHjhz6tUKifCm8sqFAoX/pxSqRQHBwenfmx+fj62trYu3Of29vapQdVpWq1WPHr0KB4/ftyzzeF4rdVqR20rlUpUq9WYnp6O7e3t2NraikqlcuJzHzx4EI1Gw3gHAABIQbFYjMnJyZ5tms1mdDqdAVV0UqHb7XZTuXIfji9lXF5ejkKhEHNzcxfup9VqxXvvvRetVis+/fTTJEvkgqampmJ/fz8mJyfjyZMnaZeTiFarFY1GIx49etQzRKpWq7GyshL37t07OhE1CY1GI2q1WqysrJw5+6lcLsfKykpUq9VzD59otVqxs7MTq6urUavVzmy3vLwcb7zxRpTL5b6/npWVlVMDsLm5uaO+DvtrNBrRaDRic3PzaPye9bXt7e2d+rHp6ekTewzOzc3F5uZmX7UCAABc1aediJvf6N3mkzcjbmR646rkpZkNZDoIu3fvXnzwwQdH97vd7qVmwBz/XEFYukYxCDtucXExNjY2Xnj/IAKYRqMR09PTp37s4ODgUuHbzMzMqTOr1tfXY2Fh4UK1zc7OvnDwxcLCQqyvr/fVx+PHj2NlZeWF9/cKwp7/eQjCAACAQRKEnS7NbCDT3+rXX389ut1uHGZ1lw3BYFDOWiZ5meWTF1Uul19YJhjxLPy57Ay0s+q+SAhWq9VemJlVKpVid3e37xAs4tkMtL29vRe+lqdPn575OUnOvAMAAGD4ZToIO3wRXigUjgKxy95gEM5adnjecsTrvP5Vrn3a514kXKrX6zE7O/vC5//gBz84NbTrp54PP/zwRA02wwcAAKBfmd4s/9atW1GpVOKDDz6I1dXVqFQqcfv27Qv302g04q233oof/ehHyRcJx5wVEl3mcTvsGo1G3L9//4X3b25uXioEO3Q4m+ysZaDD7HDz/0EFp8OmVqv1fSgDAADAaTIdhEVEvPHGGzE9PR1f/epXL93H3bt347XXXstlGEG+nPYYf/nllxPtr99xND8//8JsrYWFhUSCjHK5HMvLyz1PnsySer0e6+vrUavVotFoRKlUitu3b0epVIpyuRxvvPFGRDw78XJ1dbXn0tNarRbz8/Nx+/btWF9fH+lg6PAgiM3NzaODG7I0wzdPPwsAABgVmQ/CqtVqvP/++4n0de/
"text/plain": [
"<Figure size 1200x900 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
2024-02-29 15:54:19 +01:00
"bins = np.linspace(-1.0, 1.0, 50)\n",
"sns.regplot(\n",
" x=ak.to_numpy(array[\"dSlope_xEndT\"]),\n",
" y=ak.to_numpy(array[\"yDiffOut\"]),\n",
" x_bins=bins,\n",
" fit_reg=None,\n",
" x_estimator=np.mean,\n",
")\n",
2024-02-29 15:54:19 +01:00
"plt.xlabel(\"$\\Delta t_x$\")\n",
"plt.ylabel(\"$y_{corr}$ [mm]\")\n",
"mplhep.lhcb.text(\"Simulation\")\n",
2024-02-29 15:54:19 +01:00
"plt.show()\n",
"# plt.savefig(\n",
"# \"/work/cetin/LHCb/reco_tuner/parameterisations/plots/bend_y_deltatx_dist.pdf\",\n",
"# format=\"PDF\",\n",
"# )"
]
},
{
"cell_type": "code",
2024-02-29 15:54:19 +01:00
"execution_count": 55,
2023-12-19 13:00:59 +01:00
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"intercept= 0.16257747434812642\n",
"coef= {'dSlope_yEndT': 4379.8079670973475, 'dSlope_xEndT dSlope_yEndT': 37.76975769672208, 'ideal_state_770_ty dSlope_xEndT^2': 1224.0422108322969, 'ideal_state_770_ty dSlope_yEndT^2': -18042.13200642719}\n",
"r2 score= 0.8441034202616307\n",
"RMSE = 10.185808113150346\n"
2023-12-19 13:00:59 +01:00
]
}
],
"source": [
"features = [\n",
" \"ideal_state_770_ty\",\n",
" # \"ideal_state_9410_ty\",\n",
" # \"ideal_state_770_y\",\n",
" # \"ideal_state_9410_y\",\n",
" \"dSlope_xEndT\",\n",
" \"dSlope_yEndT\",\n",
"]\n",
"target_feat = \"yDiffOut\"\n",
"order = 3\n",
"\n",
"data = np.column_stack([ak.to_numpy(array[feat]) for feat in features])\n",
"target = ak.to_numpy(array[target_feat])\n",
"X_train, X_test, y_train, y_test = train_test_split(data,\n",
" target,\n",
" test_size=0.2,\n",
" random_state=42)\n",
"\n",
"poly = PolynomialFeatures(degree=order, include_bias=False)\n",
"X_train_model = poly.fit_transform(X_train)\n",
"X_test_model = poly.fit_transform(X_test)\n",
"\n",
"poly_features = poly.get_feature_names_out(input_features=features)\n",
"\n",
"# keep = [\n",
"# # \"ideal_state_770_ty\", #\n",
"# # \"dSlope_xEndT\", #\n",
"# \"dSlope_yEndT\", # keep\n",
"# # \"ideal_state_770_ty^2\", #\n",
"# # \"ideal_state_770_ty dSlope_xEndT\", #\n",
"# # \"ideal_state_770_ty dSlope_yEndT\", #\n",
"# # \"dSlope_xEndT^2\", # do not keep\n",
"# \"dSlope_xEndT dSlope_yEndT\", # keep\n",
"# # \"dSlope_yEndT^2\", #\n",
"# # \"ideal_state_770_ty^3\",\n",
"# # \"ideal_state_770_ty^2 dSlope_xEndT\", #\n",
"# # \"ideal_state_770_ty^2 dSlope_yEndT\",\n",
"# \"ideal_state_770_ty dSlope_xEndT^2\", # keep\n",
"# # \"ideal_state_770_ty dSlope_xEndT dSlope_yEndT\", #\n",
"# \"ideal_state_770_ty dSlope_yEndT^2\", # keep\n",
"# # \"dSlope_xEndT^3\", # do not keep\n",
"# # \"dSlope_xEndT^2 dSlope_yEndT\", #\n",
"# # \"dSlope_xEndT dSlope_yEndT^2\", #\n",
"# # \"dSlope_yEndT^3\",\n",
"# # \"ideal_state_770_ty^4\",\n",
"# # \"ideal_state_770_ty^3 dSlope_xEndT\",\n",
"# # \"ideal_state_770_ty^3 dSlope_yEndT\",\n",
"# # \"ideal_state_770_ty^2 dSlope_xEndT^2\", #\n",
"# # \"ideal_state_770_ty^2 dSlope_xEndT dSlope_yEndT\",\n",
"# # \"ideal_state_770_ty^2 dSlope_yEndT^2\",\n",
"# # \"ideal_state_770_ty dSlope_xEndT^3\", # do not keep\n",
"# # \"ideal_state_770_ty dSlope_xEndT^2 dSlope_yEndT\", #\n",
"# # \"ideal_state_770_ty dSlope_xEndT dSlope_yEndT^2\",\n",
"# # \"ideal_state_770_ty dSlope_yEndT^3\",\n",
"# # \"dSlope_xEndT^4\", #\n",
"# \"dSlope_xEndT^3 dSlope_yEndT\", # keep\n",
"# # \"dSlope_xEndT^2 dSlope_yEndT^2\", #\n",
"# # \"dSlope_xEndT dSlope_yEndT^3\",\n",
"# # \"dSlope_yEndT^4\",\n",
"# # \"ideal_state_770_ty^5\",\n",
"# # \"ideal_state_770_ty^4 dSlope_xEndT\",\n",
"# # \"ideal_state_770_ty^4 dSlope_yEndT\",\n",
"# # \"ideal_state_770_ty^3 dSlope_xEndT^2\",\n",
"# # \"ideal_state_770_ty^3 dSlope_xEndT dSlope_yEndT\",\n",
"# # \"ideal_state_770_ty^3 dSlope_yEndT^2\",\n",
"# # \"ideal_state_770_ty^2 dSlope_xEndT^3\", #\n",
"# # \"ideal_state_770_ty^2 dSlope_xEndT^2 dSlope_yEndT\",\n",
"# # \"ideal_state_770_ty^2 dSlope_xEndT dSlope_yEndT^2\",\n",
"# # \"ideal_state_770_ty^2 dSlope_yEndT^3\",\n",
"# \"ideal_state_770_ty dSlope_xEndT^4\", # keep\n",
"# # \"ideal_state_770_ty dSlope_xEndT^3 dSlope_yEndT\", #\n",
"# # \"ideal_state_770_ty dSlope_xEndT^2 dSlope_yEndT^2\",\n",
"# # \"ideal_state_770_ty dSlope_xEndT dSlope_yEndT^3\",\n",
"# # \"ideal_state_770_ty dSlope_yEndT^4\",\n",
"# # \"dSlope_xEndT^5\", #\n",
"# \"dSlope_xEndT^4 dSlope_yEndT\", # keep\n",
"# # \"dSlope_xEndT^3 dSlope_yEndT^2\", #\n",
"# # \"dSlope_xEndT^2 dSlope_yEndT^3\",\n",
"# # \"dSlope_xEndT dSlope_yEndT^4\",\n",
"# # \"dSlope_yEndT^5\",\n",
"# ]\n",
"\n",
"keep = [\n",
" \"dSlope_yEndT\", # keep\n",
" \"dSlope_xEndT dSlope_yEndT\", # keep\n",
" \"ideal_state_770_ty dSlope_xEndT^2\", # keep\n",
" \"ideal_state_770_ty dSlope_yEndT^2\", # keep\n",
" # \"dSlope_xEndT^3 dSlope_yEndT\", # do not keep\n",
" # \"ideal_state_770_ty dSlope_xEndT^4\", # keep\n",
" # \"dSlope_xEndT^4 dSlope_yEndT\", # keep\n",
"]\n",
"\n",
"# keep = [\n",
"# \"ideal_state_770_ty dSlope_yEndT^2\",\n",
"# \"ideal_state_770_ty dSlope_xEndT^2\",\n",
"# ]\n",
"remove = [i for i, f in enumerate(poly_features) if f not in keep]\n",
"\n",
"X_train_model = np.delete(X_train_model, remove, axis=1)\n",
"X_test_model = np.delete(X_test_model, remove, axis=1)\n",
"poly_features = np.delete(poly_features, remove)\n",
"# print(poly_features)\n",
"lin_reg = LinearRegression()\n",
"# lin_reg = Lasso(fit_intercept=False, alpha=0.000001)\n",
"# lin_reg = Lasso(alpha=0.1)\n",
"# lin_reg = LassoCV(max_iter=2000)\n",
"# lin_reg = ElasticNet(alpha=0.1)\n",
"# lin_reg = Ridge(alpha=0)\n",
"lin_reg.fit(X_train_model, y_train)\n",
"y_pred_test = lin_reg.predict(X_test_model)\n",
"print(\"intercept=\", lin_reg.intercept_)\n",
"print(\"coef=\", dict(zip(poly_features, lin_reg.coef_)))\n",
"print(\"r2 score=\", lin_reg.score(X_test_model, y_test))\n",
"print(\"RMSE =\", mean_squared_error(y_test, y_pred_test, squared=False))"
]
},
{
"cell_type": "code",
2024-02-29 15:54:19 +01:00
"execution_count": 45,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
2024-02-29 15:54:19 +01:00
"['dSlope_xEndT^2', 'dSlope_xEndT^3']"
]
},
2024-02-29 15:54:19 +01:00
"execution_count": 45,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"koeffs = dict(zip(poly_features, lin_reg.coef_))\n",
2023-12-19 13:00:59 +01:00
"\n",
"remove = []\n",
"for itr in koeffs.items():\n",
" if abs(itr[1]) < 1.0:\n",
" remove.append(itr[0])\n",
"remove"
]
},
{
"cell_type": "code",
2024-02-29 15:54:19 +01:00
"execution_count": 50,
"metadata": {},
"outputs": [],
"source": [
"ls_koeffs = []\n",
"for itr in koeffs.items():\n",
" ls_koeffs.append(itr[0])\n",
"# ls_koeffs"
]
},
{
"cell_type": "code",
"execution_count": 65,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"['1', 'dSlope_xEndT^2', 'dSlope_xEndT^3', 'dSlope_xEndT^6']"
]
},
"execution_count": 65,
"metadata": {},
"output_type": "execute_result"
}
],
"source": []
},
{
"cell_type": "code",
2024-02-29 15:54:19 +01:00
"execution_count": 59,
"metadata": {},
"outputs": [
{
"data": {
2024-02-29 15:54:19 +01:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAABMIAAAOWCAYAAAANzz7PAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAADlFklEQVR4nOz9X2zb6X3nfX8oT7sp2og/yT3YhedGrZ/s3ew+uDfWj3KAFs8ejEVm9sE2aCcm5UHzJD2YEWmnxR44MSlNsEc3xhLliQ520cakZk5SBKhE2tOivXEnJqWcPECxY5FynhtY5O6EPyl4bOwN7FqkmBYN2ox+z4FCjqg/FCn++FfvF0BUfy5e15em7Eaf+V7X5XEcxxEAAAAAAAAw4Ia6XQAAAAAAAADQCQRhAAAAAAAAOBcIwgAAAAAAAHAuEIQBAAAAAADgXCAIAwAAAAAAwLlAEAYAAAAAAIBzgSAMAAAAAAAA5wJBGAAAAAAAAM4FgjAAAAAAAACcCwRhAAAAAAAAOBcIwgAAAAAAAHAuEIQBAAAAAADgXCAIAwAAAAAAwLlAEAYAAAAAAIBzgSAMAAAAAAAA5wJBGAAAAAAAAM4FgjAAAAAAAACcCwRhAAAAAAAAOBcIwgAAAAAAAHAuEIQBAAAAAADgXCAIAwAAAAAAwLlAEAYAAAAAAIBzgSAMAAAAAAAA5wJBGAAAAAAAAM4FgjAAAAAAAACcCwRhAAAAAAAAOBcIwgAAAAAAAHAuEIQB6LpSqdTtEgAAAAAA5wBBGICuyOfzisViGhkZ0czMTLfLAQAAAACcA690uwAAvSefz2tlZUX5fF62bWtnZ0elUkmGYcg0TY2OjsqyLAUCAfn9/iPPD4VCsm1buVyu+rVSqaRsNquVlRVls1m6wAAAAAAAHedxHMfpdhGA25LJpCKRSNPPsyyrJrw5K4/H0/RzDMNQsVg89nuhUEjpdLrpOTOZzLFB1XFKpZLm5+e1uLjY9DqWZcnv92t8fFyZTEbpdPrIn2UlHMvn80eeHwwGlUqlml4XAAAAAIBmsDUSAykcDqtYLCqXyykYDNYd6/f7lclkquPd4DiOCoWCEomEDMM4cZxpmkokEioUCieGYJKUSqVULBYbCrai0ahyuZyKxWLDIVhli+LhEKwSUFXmq7yuTCajcDhcfW35fF6Li4uKRCLVwO5wx1dlHsdxZJpmQ3UBAAAAAOAmOsJwLkQiESWTySNf70Qnkm3bGh8fP/Z7xWKxblB2Ep/Pd2xnVSKRUDgcbqq2QCAg27Zrvh4Oh5VIJBqaY3FxUbFY7MjXTdNUoVA49jmH3w86wgAAAAAAnUBHGM6Fk7ZJnmX7ZLNM05RlWUe+HgwGzxSCSSfX3UwIls1mNT4+XhOCGYahXC7XcAgm7XegFQqFI69lZ2fnxOec9XUDAAAAANAKgjCcCydtxevUFr3j1mll7eOe20y4lM/nFQgEjjx/bW3t2NCukXq2trZqauAwfAAAAABAryEIw7lwUkg0Ojra2UJ6gG3bmpqaOvL1VCp1phCsotJNNmjS6fSRraP4VDab7XYJAAAAANAwgjCgA44L3C5evOjqfI2GeqFQ6Ei3Vjgcbvhg/XpM01Q0Gm15nk7I5/OKRCIaHx+Xx+PRyMiIxsfH5fP5qrd0ptNpzczMnBr2ZLPZ6vMHPRiybVvJZFKBQEAej+dIZ2G3naf3AgAAAEDzXul2AQA6J5vNHnvIfjwed22Nubm56u2TpVKp584DK5VKmpmZqd5uaZpm9WZR27aVz+eVz+er3688p55KuFgqlRSJRE68JKCf2batUCgk27Z7etvreXgvAAAAAJwdHWHAOXLcIfutHNp/HMMwmjq0v9Mq3V7S/i2bhUJBqVRKqVRKuVxOhULhSHfcaWHKedhia5qmcrmcisViT3f9nYf3AgAAAMDZEYQB50Q2mz32rKu5uTnX1wqFQpLq3xxZTzabVSgUqtm2GAgElEwmW6prcXGxul0uGo0eG9iZpqlMJlPTJXfaGWGJREKmaco0TaVSqZZq7Ae9th3yoPP2XgAAAABoDlsjgXPipFCglQPyT+L3+xWPx5vuzimVSpqamjqyfbNUKimbzSqbzSoej5/5YP/5+fnqx8d1xx0UjUb18uVLLS4unhro+f3+c7UFr1tdV5FIRIlEou6Y8/ZeAAAAAGgOHWHAOXHcweFuHJB/kmg02tSWS9u2NTY2duwZZofH+Xy+pg9CP8vZVvF4XJZlcWtkD0gmky13BAIAAAAAQRhwDti2fWyY045usLMolUrVQ86DwaASiYRyuZwymcyJ51EFAoGmAqrDIVijQdrc3FxPHw5/Hti2fWoHHwAAAAA0gq2RwDlwUmA0Pj7e4UqOVwmlMpnMkS41v9+vSCRybPAVCoWUy+UaWsM0zZrPY7GYpqenT+1aq9wo2Ys3YJ4HpVKpp88kAwAAANBf6AgDzoGTgrBeumHvuBCsonKA/WH5fL7hzi7DMGqCrFKpJJ/P11BXmeM4hGBdUDkzjq2pAAAAANxCEAZ0SSwWk8fjOdPD5/M1tdZJW/t6JdwxDOPU88pM06y5ybHiuK+d5PBY27Y1Pj7uytlT+XxekUjk1C67UqmkZDIpn89Xs25l+9/IyIg8Ho/Gx8e1uLh45Pm2bSsWi1Vv1BwfH1ckEjn2Pa6sc/jn57BAIFCd7+DDzS2hyWSyZp2RkRH5fL5jX6O0/+d53JlxB+s76c+60feiolQqaXFxUYFAoPrnX6kvFos1FMTZtq3FxUWNj4/X1HzwPRgfH69uAQYAAADQJQ5wTkg68igWix1ZOxwOH7u+mw/TNJteP5PJdOT1HxaNRmvq8Pv9DT/XMIwjr6MZfr//2D8Ly7KcXC7X1Fy5XM4Jh8M1NRmGcWRcsVh04vG4Y1lWzZqJRMJxnKN/Hgcf4XC4Ok+9caZpnvjzfPj9P0k8Hm/o70cul2v4zz+VSlX/XFKpVM0cpmmeWHuxWHRyudyRmnK5XM3j4HyNvBeHJRKJ6vh4PO7kcjmnUCg4qVSq5v0KBoNHnlssFp1oNFp9HQdrLBaLR97vg3UVCoVTawMAAADgPjrCgC4Jh8PK5XJneiQSiabWOqnzq1c6U5rpTDuuc+y0myYPymQy1XO/Ds/h8/ka7tixbVsrKysNdx1ZlnXkwPdcLiefz6d8Pq9UKqVCoaBisVhzQUAymVQ+n1cgEFA2m1Umk1GxWFSxWFQmk6n+2dm2rfn5+WPXDoVCDdV43J9LK7LZbHXtubm5mvkty1IqlZK0X/vMzEzNcw3DkGVZRy50qHzt4PeafS8qYrGYIpGIDMNQoVBQNBqVZVkyTVPBYFC5XK5aczqd1vj4eM3Pxs7OjgKBwJGfyY2NDY2NjVW39BYKBSUSiep7VSqVOPwfAAAA6JZuJ3FAp6jHOsLi8fiZ5zvckaNTOsIOd9VUHpWOpE473Nl0XLfNSSodRq2+jpP+TM4y58G5TutCOtwldFJX3uHOtYOdYQcd/PM4ae1MJtNQB1exWHS1IywYDJ76837aPI3WXtHoe3GwE+y09/pgx9dJ3YsHazQM49j39fDPbqf+/QEAAADwKTrCgHOg1zvCmnH49kfpbK8jGo2qWCye2AVVuamykbmPq+kkBy8oCIfDJ56NdvimxJO6AA/WXyqVevY9Peln8ODX3TgUv9H3IhaLVdcPh8MNjZX2u9yOu6Dh4OtYXl4+9n09/LO2sbHRUK0AAAAA3EMQBpwDJ4UDhUKhw5W07rhA5eXLl2eeK5VKKZfLHdmCJ+2HHlNTU66GS41uAz34njWzdbSXblhcXl5WPB5XKpU6NWySOhfMJpPJ6lqNBGfhcLjmPWjmgobDDq7XS+8VAAAAcF4QhAHnwEldR8d1tvS644KLixcvtjSnZVknnr2Wz+ePnF/VCWe90XNnZ8fdQlpgGIai0eiJZ7I1c7abmzKZTPXjRjvIDv4d6se/NwAAAAD2EYQB58RxHU+2bffsVrpmNLM1sZ5wOKxCoXAkhEqn0z0dfpw1NOs027aVTCaVzWaPPQi/k3U06/DP2CD8vQEAAADOI4Iw4Jy4devWsV/v5YD
"text/plain": [
"<Figure size 1200x900 with 2 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
2024-02-29 15:54:19 +01:00
"bins = np.linspace(-40, 40, 41)\n",
"ax = sns.regplot(\n",
" x=y_test,\n",
2024-02-29 15:54:19 +01:00
" y=y_test - y_pred_test,\n",
" x_bins=bins,\n",
" fit_reg=None,\n",
" x_estimator=np.mean,\n",
" label=\"bla\",\n",
")\n",
"ax2 = ax.twinx()\n",
2024-02-29 15:54:19 +01:00
"ax2.hist(\n",
" y_test,\n",
" bins=30,\n",
" range=[-40, 40],\n",
" color=\"#2A9D8F\",\n",
" alpha=0.5,\n",
" align=\"mid\",\n",
" density=True,\n",
")\n",
"ax.set_xlabel(r\"$y_{\\mathrm{corr}}$ [mm]\")\n",
"ax.set_ylabel(\n",
" r\"$\\left\\langle y_{\\mathrm{corr}}-y_{\\mathrm{corr}}^{\\mathrm{pred}}\\right\\rangle$ [mm]\"\n",
")\n",
"ax2.set_ylabel(\"Number of Tracks (normalised)\")\n",
"mplhep.lhcb.text(\"Simulation\", loc=0)\n",
2024-02-29 15:54:19 +01:00
"# plt.show()\n",
"plt.savefig(\n",
" \"/work/cetin/LHCb/reco_tuner/parameterisations/plots/bend_y_regression_plot.pdf\",\n",
" format=\"PDF\",\n",
")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
2024-02-29 15:54:19 +01:00
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
2023-12-19 13:00:59 +01:00
"features = [\n",
" \"ideal_state_770_ty\",\n",
" \"ideal_state_770_tx\",\n",
2023-12-19 13:00:59 +01:00
" \"dSlope_xEndT\",\n",
" \"dSlope_yEndT\",\n",
" \"dSlope_xEndT_abs\",\n",
" \"dSlope_yEndT_abs\",\n",
"]\n",
"target_feat = \"yDiffEndT\"\n",
"\n",
"data = np.column_stack([ak.to_numpy(array[feat]) for feat in features])\n",
"target = ak.to_numpy(array[target_feat])\n",
2024-02-29 15:54:19 +01:00
"X_train, X_test, y_train, y_test = train_test_split(data,\n",
" target,\n",
" test_size=0.2,\n",
" random_state=42)\n",
2023-12-19 13:00:59 +01:00
"\n",
"poly = PolynomialFeatures(degree=6, include_bias=False)\n",
"X_train_model = poly.fit_transform(X_train)\n",
"X_test_model = poly.fit_transform(X_test)\n",
2023-12-19 13:00:59 +01:00
"poly_features = poly.get_feature_names_out(input_features=features)\n",
"keep = [\n",
" #'dSlope_xEndT',\n",
" \"dSlope_yEndT\", # keep\n",
2023-12-19 13:00:59 +01:00
" #'dSlope_yEndT_abs',\n",
" #'ty dSlope_xEndT',\n",
" #'ty dSlope_yEndT',\n",
" \"ideal_state_770_ty dSlope_xEndT_abs\", # keep\n",
" \"ideal_state_770_ty dSlope_yEndT_abs\", # keep\n",
" \"ideal_state_770_ty dSlope_yEndT^2\", # keep\n",
" \"ideal_state_770_ty dSlope_xEndT^2\", # keep\n",
2023-12-19 13:00:59 +01:00
" #'tx dSlope_xEndT',\n",
" #'tx dSlope_xEndT_abs',\n",
" #'tx dSlope_yEndT',\n",
" \"ideal_state_770_ty ideal_state_770_tx dSlope_xEndT\", # keep\n",
" \"ideal_state_770_tx^2 dSlope_yEndT\", # keep\n",
2023-12-19 13:00:59 +01:00
" #'ty^2 dSlope_xEndT',\n",
" #'ty^2 dSlope_yEndT',\n",
2023-12-19 13:00:59 +01:00
" #'ty^2 dSlope_xEndT_abs',\n",
" #'ty^2 tx dSlope_xEndT',\n",
" #'ty tx^2 dSlope_yEndT',\n",
" \"ideal_state_770_ty ideal_state_770_tx^2 dSlope_xEndT_abs\", # keep\n",
" \"ideal_state_770_ty^3 ideal_state_770_tx dSlope_xEndT\", # keep\n",
2023-12-19 13:00:59 +01:00
" #'ty tx^3 dSlope_xEndT',\n",
" #'ty^3 dSlope_yEndT_abs',\n",
"]\n",
"do_not_keep = [\n",
" \"dSlope_xEndT\",\n",
" \"dSlope_yEndT_abs\",\n",
" \"ideal_state_770_ty dSlope_xEndT\",\n",
" \"ideal_state_770_tx dSlope_xEndT\",\n",
" \"ideal_state_770_tx dSlope_xEndT_abs\",\n",
" \"ideal_state_770_tx dSlope_yEndT\",\n",
" \"ideal_state_770_ty^2 dSlope_xEndT\",\n",
" \"ideal_state_770_ty^3 dSlope_yEndT_abs\",\n",
" \"ideal_state_770_ty ideal_state_770_tx dSlope_yEndT\",\n",
" \"ideal_state_770_ty ideal_state_770_tx^3 dSlope_xEndT\",\n",
" \"ideal_state_770_ty ideal_state_770_tx^2 dSlope_yEndT\",\n",
2023-12-19 13:00:59 +01:00
"]\n",
"reduce = True\n",
"if reduce:\n",
2024-02-29 15:54:19 +01:00
" remove = [\n",
" i for i, f in enumerate(poly_features) if (keep and f not in keep)\n",
" ]\n",
" X_train_model = np.delete(X_train_model, remove, axis=1)\n",
" X_test_model = np.delete(X_test_model, remove, axis=1)\n",
" poly_features = np.delete(poly_features, remove)\n",
2023-12-19 13:00:59 +01:00
" print(poly_features)\n",
"if not reduce:\n",
" remove = [\n",
2024-02-29 15:54:19 +01:00
" i for i, f in enumerate(poly_features)\n",
" if (\"dSlope_\" not in f) or (\"EndT^\" in f) or (\"abs^\" in f) or (\n",
" \"EndT dSlope\" in f) or (\"abs dSlope\" in f)\n",
" ]\n",
" X_train_model = np.delete(X_train_model, remove, axis=1)\n",
" X_test_model = np.delete(X_test_model, remove, axis=1)\n",
" poly_features = np.delete(poly_features, remove)\n",
" # print(poly_features)\n",
2023-12-19 13:00:59 +01:00
" lin_reg = Lasso(fit_intercept=False, alpha=0.000001)\n",
"else:\n",
" lin_reg = LinearRegression(fit_intercept=False)\n",
"lin_reg.fit(X_train_model, y_train)\n",
"y_pred_test = lin_reg.predict(X_test_model)\n",
2023-12-19 13:00:59 +01:00
"print(\"intercept=\", lin_reg.intercept_)\n",
"print(\n",
" \"coef=\",\n",
" {\n",
" k: v\n",
" for k, v in zip(poly_features, lin_reg.coef_)\n",
" if abs(v) > 1.0 and k not in keep and k not in do_not_keep\n",
" },\n",
")\n",
2023-12-19 13:00:59 +01:00
"print(\"r2 score=\", lin_reg.score(X_test_model, y_test))\n",
"print(\"RMSE =\", mean_squared_error(y_test, y_pred_test, squared=False))\n",
"print(\n",
" \"straight RMSE =\",\n",
" mean_squared_error(\n",
" array[\"ideal_state_9410_y\"],\n",
2024-02-29 15:54:19 +01:00
" array[\"ideal_state_770_y\"] + array[\"ideal_state_770_ty\"] *\n",
" (array[\"ideal_state_9410_z\"] - array[\"ideal_state_770_z\"]),\n",
" squared=False,\n",
" ),\n",
")\n",
2023-12-19 13:00:59 +01:00
"print(format_array(\"y_xEndT_diff\", lin_reg.coef_))"
]
}
],
"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.12"
2023-12-19 13:00:59 +01:00
},
"orig_nbformat": 4,
"vscode": {
"interpreter": {
"hash": "a2eff8b4da8b8eebf5ee2e5f811f31a557e0a202b4d2f04f849b065340a6eda6"
}
}
},
"nbformat": 4,
"nbformat_minor": 2
}