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

213 lines
21 KiB
Plaintext
Raw Permalink Normal View History

2023-12-19 13:00:59 +01:00
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import uproot\n",
"import awkward as ak\n",
"import matplotlib\n",
"import matplotlib.pyplot as plt\n",
"input_tree = uproot.open({\"/work/guenther/reco_tuner/data/param_data_selected_all_p.root\": \"Selected\"})\n",
"array = input_tree.arrays()\n",
"array[\"dSlope_fringe\"] = array[\"tx_ref\"] - array[\"tx\"]\n",
"array[\"poqmag_gev\"] = 1. / ( array[\"signed_rel_current\"] * array[\"qop\"] * 1000. )\n",
"array[\"B_integral\"] = array[\"poqmag_gev\"] * array[\"dSlope_fringe\"]"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(array([1.00000e+00, 0.00000e+00, 0.00000e+00, 1.00000e+00, 0.00000e+00,\n",
" 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,\n",
" 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 1.00000e+00,\n",
" 1.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,\n",
" 0.00000e+00, 0.00000e+00, 1.00000e+00, 0.00000e+00, 0.00000e+00,\n",
" 2.00000e+00, 2.00000e+00, 4.00000e+00, 0.00000e+00, 3.00000e+00,\n",
" 0.00000e+00, 4.00000e+00, 3.00000e+00, 1.00000e+00, 4.00000e+00,\n",
" 4.00000e+00, 1.00000e+00, 1.00000e+01, 8.00000e+00, 8.00000e+00,\n",
" 2.30000e+01, 3.20000e+01, 6.30000e+01, 8.40000e+01, 1.22000e+02,\n",
" 2.03000e+02, 3.10000e+02, 3.98000e+02, 5.90000e+02, 8.65000e+02,\n",
" 1.08700e+03, 1.55300e+03, 2.01000e+03, 2.67400e+03, 3.49500e+03,\n",
" 4.46300e+03, 5.85300e+03, 7.16800e+03, 8.92900e+03, 1.12130e+04,\n",
" 1.37320e+04, 1.69420e+04, 2.09870e+04, 2.59060e+04, 3.24750e+04,\n",
" 3.94560e+04, 4.81480e+04, 6.21400e+04, 8.28540e+04, 1.20767e+05,\n",
" 2.15388e+05, 5.13195e+05, 5.20000e+02, 3.29000e+02, 1.83000e+02,\n",
" 1.32000e+02, 1.02000e+02, 8.30000e+01, 6.40000e+01, 5.00000e+01,\n",
" 2.70000e+01, 2.80000e+01, 2.20000e+01, 1.80000e+01, 2.00000e+01,\n",
" 1.50000e+01, 1.00000e+01, 4.00000e+00, 7.00000e+00, 6.00000e+00,\n",
" 5.00000e+00, 2.00000e+00, 2.00000e+00, 3.00000e+00, 0.00000e+00,\n",
" 2.00000e+00, 3.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00]),\n",
" array([-2.45885773, -2.44145641, -2.4240551 , -2.40665378, -2.38925247,\n",
" -2.37185115, -2.35444984, -2.33704853, -2.31964721, -2.3022459 ,\n",
" -2.28484458, -2.26744327, -2.25004195, -2.23264064, -2.21523933,\n",
" -2.19783801, -2.1804367 , -2.16303538, -2.14563407, -2.12823275,\n",
" -2.11083144, -2.09343013, -2.07602881, -2.0586275 , -2.04122618,\n",
" -2.02382487, -2.00642355, -1.98902224, -1.97162093, -1.95421961,\n",
" -1.9368183 , -1.91941698, -1.90201567, -1.88461436, -1.86721304,\n",
" -1.84981173, -1.83241041, -1.8150091 , -1.79760778, -1.78020647,\n",
" -1.76280516, -1.74540384, -1.72800253, -1.71060121, -1.6931999 ,\n",
" -1.67579858, -1.65839727, -1.64099596, -1.62359464, -1.60619333,\n",
" -1.58879201, -1.5713907 , -1.55398938, -1.53658807, -1.51918676,\n",
" -1.50178544, -1.48438413, -1.46698281, -1.4495815 , -1.43218018,\n",
" -1.41477887, -1.39737756, -1.37997624, -1.36257493, -1.34517361,\n",
" -1.3277723 , -1.31037098, -1.29296967, -1.27556836, -1.25816704,\n",
" -1.24076573, -1.22336441, -1.2059631 , -1.18856178, -1.17116047,\n",
" -1.15375916, -1.13635784, -1.11895653, -1.10155521, -1.0841539 ,\n",
" -1.06675258, -1.04935127, -1.03194996, -1.01454864, -0.99714733,\n",
" -0.97974601, -0.9623447 , -0.94494338, -0.92754207, -0.91014076,\n",
" -0.89273944, -0.87533813, -0.85793681, -0.8405355 , -0.82313418,\n",
" -0.80573287, -0.78833156, -0.77093024, -0.75352893, -0.73612761,\n",
" -0.7187263 ]),\n",
" <BarContainer object of 100 artists>)"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAicAAAGdCAYAAADJ6dNTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAj3ElEQVR4nO3df3AU9eH/8deJ5ACbnMaUQCDBVIWKkSCBYqggYWo0yg9FLVYnxg4wpYZaJnWslLEUp04orUinHEjaqdgZZ8w4VTpTGZl0qkIN1hATi6ZW6QQTSEIKhRw/aoLH+/OHX+7rkQC5y13uvbfPx8zNcHu7e+99Z9m88t73+70eY4wRAACAJS5JdAEAAAC+jHACAACsQjgBAABWIZwAAACrEE4AAIBVCCcAAMAqhBMAAGAVwgkAALDKpYkuQKTOnDmjtrY2paamyuPxJLo4AACgH4wxOn78uLKysnTJJRduG3FcOGlra1N2dnaiiwEAAKLQ2tqqsWPHXnAdx4WT1NRUSV8cXFpaWoJLAwAA+iMQCCg7Ozv0e/xCHBNO/H6//H6/gsGgJCktLY1wAgCAw/SnS4bHaQ/+CwQC8vl86urqIpwAAOAQkfz+ZrQOAACwimPCid/v18SJEzVt2rREFwUAAMQRt3UAAEDccVsHAAA4FuEEAABYxTHhhD4nAAC4A31OAABA3NHnBAAAOJZjwgm3dQAAcAdu6wAAgLjjtg4AAHAswgkAALCKY55KDACIr6ueeC3s/f61dyaoJHA7x7Sc0CEWAAB3cEw4KS8vV1NTk+rq6hJdFAAAEEeOCScAAMAdCCcAAMAqhBMAAGAVwgkAALCKY8IJo3UAAHAHx4QTRusAAOAOjgknAADAHQgnAADAKoQTAABgFcIJAACwCuEEAABYxTHhhKHEAAC4g2PCCUOJAQBwB8eEEwAA4A6EEwAAYBXCCQAAsArhBAAAWIVwAgAArEI4AQAAViGcAAAAq1ya6AIAAAbfVU+8lugiAOdFywkAALCKY8IJ09cDAOAOjgknTF8PAIA7OCacAAAAdyCcAAAAqxBOAACAVQgnAADAKoQTAABgFcIJAACwCuEEAABYhXACAACsQjgBAABWSUg4ufTSSzV58mRNnjxZS5YsSUQRAACApRLyVOLLL79cjY2NifhqAHAlnkIMJ+G2DgAAsErE4WTnzp2aN2+esrKy5PF4tG3btl7rbNq0Sbm5uRo2bJgKCgq0a9eusM8DgYAKCgp0880366233oq68AAAIPlEHE5Onjyp/Px8bdy4sc/Pq6urtWLFCq1atUoNDQ2aOXOmSkpK1NLSElpn//79qq+v13PPPaeHHnpIgUAg+iMAAABJJeJwUlJSop///OdauHBhn5+vX79eixcv1pIlS3Tddddpw4YNys7O1ubNm0PrZGVlSZLy8vI0ceJEffzxx+f9vu7ubgUCgbAXAABIXjHtc9LT06P6+noVFxeHLS8uLlZtba0k6ejRo+ru7pYkHThwQE1NTfra17523n1WVlbK5/OFXtnZ2bEsMgAAsExMw8nhw4cVDAaVmZkZtjwzM1MdHR2SpH/+85+aOnWq8vPzNXfuXP36179Wenr6efe5cuVKdXV1hV6tra2xLDIAALBMXIYSezyesPfGmNCyGTNmaO/evf3el9frldfrjWn5AACAvWLacpKRkaEhQ4aEWknO6uzs7NWaEim/36+JEydq2rRpA9oPAACwW0zDSUpKigoKClRTUxO2vKamRjNmzBjQvsvLy9XU1KS6uroB7QcAANgt4ts6J06c0L59+0Lvm5ub1djYqPT0dOXk5KiiokKlpaWaOnWqCgsLVVVVpZaWFi1btiymBQcA9I3ZYOF0EYeTPXv2qKioKPS+oqJCklRWVqatW7dq0aJFOnLkiJ566im1t7crLy9P27dv17hx4wZUUL/fL7/fr2AwOKD9AAAAu3mMMSbRhYhEIBCQz+dTV1eX0tLSEl0cALBOrFpO9q+9Myb7AaTIfn/zbB0AAGAVx4QTRusAAOAOjgknjNYBAMAdHBNOAACAO8RlhlgAwOBh6DCSjWNaTuhzAgCAOzgmnNDnBAAAd3BMOAEAAO5AOAEAAFYhnAAAAKs4JpzQIRYAAHfg2ToA4CCDOWyYZ+sglni2DgAAcCzCCQAAsArhBAAAWIVwAgAArOKYZ+v4/X75/X4Fg8FEFwUABg3PzYEbOablhOnrAQBwB8eEEwAA4A6EEwAAYBXCCQAAsIpjOsQCQLKj8yvwBVpOAACAVRwTTnjwHwAA7uCYcMJQYgAA3MEx4QQAALgD4QQAAFiF0ToAkCCMzgH6RssJAACwCuEEAABYhXACAACsQjgBAABWoUMsAAwCOr8C/eeYlhNmiAUAwB0cE06YIRYAAHdwTDgBAADuQDgBAABWIZwAAACrMFoHAOKA0TlA9Gg5AQAAViGcAAAAqxBOAACAVQgnAADAKnSIBYABovMrEFu0nAAAAKskLJycOnVK48aN02OPPZaoIgAAAAslLJw8/fTTmj59eqK+HgAAWCoh4eSTTz7RRx99pDvuuCMRXw8AACwWcTjZuXOn5s2bp6ysLHk8Hm3btq3XOps2bVJubq6GDRumgoIC7dq1K+zzxx57TJWVlVEXGgAS6aonXgt7AYitiMPJyZMnlZ+fr40bN/b5eXV1tVasWKFVq1apoaFBM2fOVElJiVpaWiRJf/rTnzR+/HiNHz9+YCUHAABJKeKhxCUlJSopKTnv5+vXr9fixYu1ZMkSSdKGDRu0Y8cObd68WZWVlXrnnXf00ksv6eWXX9aJEyd0+vRppaWl6ac//Wmf++vu7lZ3d3fofSAQiLTIAADAQWLa56Snp0f19fUqLi4OW15cXKza2lpJUmVlpVpbW7V//3796le/0tKlS88bTM6u7/P5Qq/s7OxYFhkAAFgmpuHk8OHDCgaDyszMDFuemZmpjo6OqPa5cuVKdXV1hV6tra2xKCoAALBUXGaI9Xg8Ye+NMb2WSdLDDz980X15vV55vd5YFQ0AAFgupuEkIyNDQ4YM6dVK0tnZ2as1JVJ+v19+v1/BYHBA+wGASDAaBxh8Mb2tk5KSooKCAtXU1IQtr6mp0YwZMwa07/LycjU1Namurm5A+wEAAHaLuOXkxIkT2rdvX+h9c3OzGhsblZ6erpycHFVUVKi0tFRTp05VYWGhqqqq1NLSomXLlsW04AAAIDlFHE727NmjoqKi0PuKigpJUllZmbZu3apFixbpyJEjeuqpp9Te3q68vDxt375d48aNG1BBua0DAIA7eIwxJtGFiEQgEJDP51NXV5fS0tISXRwASc7NfU72r70z0UVAEonk93fCHvwHAADQl7gMJY4HbusAGAxubikBbOGYlhNG6wAA4A6OCScAAMAdCCcAAMAqjgknfr9fEydO1LRp0xJdFAAAEEcMJQbgWnR+vTCGEiOWGEoMAAAci3ACAACsQjgBAABWcUw4oUMsAADu4JhwwiRsAAC4g2PCCQAAcAfHPFsHAAaKocOAM9ByAgAArEI4AQAAVnFMOGG0DgAA7uCYcMJoHQAA3MEx4QQAALgDo3UAJCVG5gDORcsJAACwCuEEAABYhXACAACs4phwwlBiAADcwTHhhKHEAAC4g2PCCQAAcAeGEgNICgwdBpIHLScAAMAqhBMAAGAVwgkAALAK4QQAAFiFcAIAAKxCOAEAAFZxzFBiv98vv9+vYDCY6KIASDCGDQPJzTEtJ8wQCwCAOzgmnAAAAHcgnAAAAKsQTgAAgFUIJwAAwCqEEwAAYBXHDCUG4F4MHQbchZYTAABgFcIJAACwCuEEAABYhXACAACsMujh5Pjx45o2bZomT56sG264Qb/97W8HuwgAAMBigz5aZ8SIEXrrrbc0YsQInTp1Snl5eVq4cKGuvPLKwS4KAACw0KC3nAwZMkQjRoyQJH322WcKBoMyxgx2MQAAgKUiDic7d+7UvHnzlJWVJY/Ho23btvVaZ9OmTcrNzdWwYcNUUFCgXbt2hX1+7Ngx5efna+zYsXr88ceVkZER9QEAAIDkEvFtnZMnTyo/P1/f/e53dc8
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.hist(array[\"B_integral\"], bins=100, log=True)"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"['ty^2' 'tx^2' 'tx tx_ref' 'tx_ref^2' 'ty^4' 'ty^2 tx^2' 'ty^2 tx tx_ref'\n",
" 'ty^2 tx_ref^2' 'tx^4' 'tx^3 tx_ref' 'tx_ref^4']\n",
"intercept= -1.2094486121528516\n",
"coef= {'ty^2': -2.7897043324822492, 'tx^2': -0.35976930628193077, 'tx tx_ref': -0.47138558705675454, 'tx_ref^2': -0.5600847231491961, 'ty^4': 14.009315350693472, 'ty^2 tx^2': -16.162818973243674, 'ty^2 tx tx_ref': -8.807994419844437, 'ty^2 tx_ref^2': -0.8753190393972976, 'tx^4': 2.98254201374128, 'tx^3 tx_ref': 0.9625408279466898, 'tx_ref^4': 0.10200564097830103}\n",
"r2 score= 0.9916826041227943\n",
"RMSE = 0.006014471039836984\n",
"['ty^2', 'tx^2', 'tx tx_ref', 'tx_ref^2', 'ty^4', 'ty^2 tx^2', 'ty^2 tx tx_ref', 'ty^2 tx_ref^2', 'tx^4', 'tx^3 tx_ref', 'tx_ref^4']\n"
]
}
],
"source": [
"from sklearn.preprocessing import PolynomialFeatures\n",
"from sklearn.linear_model import LinearRegression, Lasso, Ridge\n",
"from sklearn.model_selection import train_test_split\n",
"from sklearn.pipeline import Pipeline\n",
"from sklearn.metrics import mean_squared_error\n",
"import numpy as np\n",
"\n",
"features = [\n",
" \"ty\", \n",
" \"tx\",\n",
" \"tx_ref\",\n",
"]\n",
"target_feat = \"B_integral\"\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, target, test_size=0.2, random_state=42)\n",
"\n",
"poly = PolynomialFeatures(degree=5, include_bias=False)\n",
"X_train_model = poly.fit_transform( X_train )\n",
"X_test_model = poly.fit_transform( X_test )\n",
"poly_features = poly.get_feature_names_out(input_features=features)\n",
"reduce = True\n",
"if reduce:\n",
" keep = [\n",
" #'ty',\n",
" #'tx',\n",
" #'tx_ref',\n",
" #'ty tx',\n",
" #'ty tx_ref',\n",
" 'ty^2',#keep\n",
" 'tx^2',#keep\n",
" 'tx tx_ref',#keep\n",
" 'tx_ref^2',#keep\n",
" 'ty^2 tx tx_ref',#keep\n",
" 'ty^2 tx^2',#keep\n",
" 'ty^2 tx_ref^2', #keep\n",
" 'tx^4',#keep\n",
" 'ty^4',#keep\n",
" 'tx_ref^4',#keep\n",
" #'tx_ref^5',\n",
" 'tx^3 tx_ref', #keep\n",
" #'tx tx_ref^3',\n",
" #'tx^2 tx_ref^2',\n",
" #'ty tx_ref^4',\n",
" #'tx tx_ref^4',\n",
" #'tx_ref^5',\n",
" ]\n",
" remove = [i for i, f in enumerate(poly_features) if (keep and f not in keep )]\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",
"if not reduce:\n",
" lin_reg = Lasso( alpha=0.0000001)#Lasso(fit_intercept=False, alpha=0.001)\n",
"else:\n",
" lin_reg = LinearRegression()\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))\n",
"print([key for key, val in dict(zip(poly_features, lin_reg.coef_)).items() if val != 0.0])"
]
},
{
"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"
},
"orig_nbformat": 4,
"vscode": {
"interpreter": {
"hash": "a2eff8b4da8b8eebf5ee2e5f811f31a557e0a202b4d2f04f849b065340a6eda6"
}
}
},
"nbformat": 4,
"nbformat_minor": 2
}