2023-09-25 16:36:13 +02:00
{
"cells": [
{
"cell_type": "code",
2023-10-10 16:17:10 +02:00
"execution_count": 1,
2023-09-25 16:36:13 +02:00
"metadata": {},
"outputs": [],
"source": [
"import uproot\n",
"import numpy as np\n",
"import sys\n",
"import os\n",
"import matplotlib\n",
"import matplotlib.pyplot as plt\n",
"from mpl_toolkits import mplot3d\n",
"import itertools\n",
"import awkward as ak\n",
"from scipy.optimize import curve_fit\n",
"%matplotlib inline"
]
},
{
"cell_type": "code",
2023-10-10 16:17:10 +02:00
"execution_count": 2,
2023-09-25 16:36:13 +02:00
"metadata": {},
2023-09-26 15:12:03 +02:00
"outputs": [
{
"data": {
"text/plain": [
2023-09-27 13:00:51 +02:00
"51"
2023-09-26 15:12:03 +02:00
]
},
2023-10-10 16:17:10 +02:00
"execution_count": 2,
2023-09-26 15:12:03 +02:00
"metadata": {},
"output_type": "execute_result"
}
],
2023-09-25 16:36:13 +02:00
"source": [
"file = uproot.open(\"tracking_losses_ntuple_Dst0ToD0EE.root:PrDebugTrackingLosses.PrDebugTrackingTool/Tuple;1\")\n",
"\n",
"\n",
"#selektiere nur elektronen von D*0->D0ee und nur solche mit einem momentum von unter 5 GeV \n",
"allcolumns = file.arrays()\n",
"found = allcolumns[(allcolumns.isElectron) & (~allcolumns.lost) & (allcolumns.fromSignal) & (allcolumns.p < 5e3)] #D: 2591\n",
"lost = allcolumns[(allcolumns.isElectron) & (allcolumns.lost) & (allcolumns.fromSignal) & (allcolumns.p < 5e3)] #D: 1908\n",
"\n",
2023-09-26 15:12:03 +02:00
"#ak.num(lost, axis=0)\n",
2023-09-27 10:38:06 +02:00
"ak.num(allcolumns[(allcolumns.fromPairProd) & (allcolumns.isElectron) & (~allcolumns.lost) & (allcolumns.fromSignal)],axis=0)\n",
2023-09-27 13:00:51 +02:00
"#found[0]"
2023-09-25 16:36:13 +02:00
]
},
{
"cell_type": "code",
2023-10-10 16:17:10 +02:00
"execution_count": 3,
2023-09-25 16:36:13 +02:00
"metadata": {},
"outputs": [
{
2023-10-10 11:00:25 +02:00
"name": "stdout",
"output_type": "stream",
"text": [
"eff all = 0.5759057568348522 +/- 0.007367987865301135\n"
]
2023-09-25 16:36:13 +02:00
}
],
"source": [
2023-10-10 11:00:25 +02:00
"def t_eff(found, lost, axis = 0):\n",
" sel = ak.num(found, axis=axis)\n",
" des = ak.num(lost, axis=axis)\n",
" return sel/(sel + des)\n",
"\n",
"def eff_err(found, lost):\n",
" n_f = ak.num(found, axis=0)\n",
" n_all = ak.num(found, axis=0) + ak.num(lost,axis=0)\n",
" return 1/n_all * np.sqrt(np.abs(n_f*(1-n_f/n_all)))\n",
"\n",
2023-09-25 16:36:13 +02:00
"\n",
2023-10-10 11:00:25 +02:00
"print(\"eff all = \", t_eff(found, lost), \"+/-\", eff_err(found, lost))"
2023-09-25 16:36:13 +02:00
]
},
{
"cell_type": "code",
2023-10-10 16:17:10 +02:00
"execution_count": 4,
2023-09-25 16:36:13 +02:00
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.7960893854748603"
]
},
2023-10-10 16:17:10 +02:00
"execution_count": 4,
2023-09-25 16:36:13 +02:00
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
2023-09-28 12:12:21 +02:00
"\n",
"#statistics\n",
"\n",
2023-09-25 16:36:13 +02:00
"nobrem_found = found[found[\"brem_photons_pe_length\"]==0]\n",
"nobrem_lost = lost[lost[\"brem_photons_pe_length\"]==0]\n",
"\n",
"\"\"\"\n",
"die effizienz mit der wir elektronen finden, die keine bremsstrahlung gemacht haben, ist nicht besonders gut, aber trotzdem besser als\n",
"für alle elektronen.\n",
"Auch hier handelt es sich um eine recht geringe sample size (~350)\n",
"\"\"\"\n",
"\n",
"t_eff(nobrem_found, nobrem_lost)\n",
"#ak.num(nobrem_lost, axis=0)"
]
},
{
"cell_type": "code",
2023-10-10 16:17:10 +02:00
"execution_count": 5,
2023-10-10 11:00:25 +02:00
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<pre>{energy: 2.04e+03,\n",
" photon_length: 4,\n",
" brem_photons_pe: [519, 484, 173, 466],\n",
" brem_vtx_z: [9.84e+03, 9.84e+03, 9.85e+03, 9.85e+03]}\n",
"------------------------------------------------------\n",
"type: {\n",
" energy: float64,\n",
" photon_length: int64,\n",
" brem_photons_pe: var * float64,\n",
" brem_vtx_z: var * float64\n",
"}</pre>"
],
"text/plain": [
"<Record {energy: 2.04e+03, ...} type='{energy: float64, photon_length: int6...'>"
]
},
2023-10-10 16:17:10 +02:00
"execution_count": 5,
2023-10-10 11:00:25 +02:00
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#try excluding all photons that originate from a vtx @ z>9500mm\n",
"#ignore all brem vertices @ z>9500mm \n",
"\n",
"#found\n",
"\n",
"brem_e_f = found[\"brem_photons_pe\"]\n",
"brem_z_f = found[\"brem_vtx_z\"]\n",
"e_f = found[\"energy\"]\n",
"length_f = found[\"brem_vtx_z_length\"]\n",
"\n",
"brem_f = ak.ArrayBuilder()\n",
"\n",
"for itr in range(ak.num(found,axis=0)):\n",
" brem_f.begin_record()\n",
" #[:,\"energy\"] energy\n",
" brem_f.field(\"energy\").append(e_f[itr])\n",
" #[:,\"photon_length\"] number of vertices\n",
" brem_f.field(\"photon_length\").integer(length_f[itr])\n",
" #[:,\"brem_photons_pe\",:] photon energy \n",
" brem_f.field(\"brem_photons_pe\").append(brem_e_f[itr])\n",
" #[:,\"brem_vtx_z\",:] brem vtx z\n",
" brem_f.field(\"brem_vtx_z\").append(brem_z_f[itr])\n",
" brem_f.end_record()\n",
"\n",
"brem_f = ak.Array(brem_f)\n",
"\n",
"#lost\n",
"\n",
"brem_e_l = lost[\"brem_photons_pe\"]\n",
"brem_z_l = lost[\"brem_vtx_z\"]\n",
"e_l = lost[\"energy\"]\n",
"length_l = lost[\"brem_vtx_z_length\"]\n",
"\n",
"brem_l = ak.ArrayBuilder()\n",
"\n",
"for itr in range(ak.num(lost,axis=0)):\n",
" brem_l.begin_record()\n",
" #[:,\"energy\"] energy\n",
" brem_l.field(\"energy\").append(e_l[itr])\n",
" #[:,\"photon_length\"] number of vertices\n",
" brem_l.field(\"photon_length\").integer(length_l[itr])\n",
" #[:,\"brem_photons_pe\",:] photon energy \n",
" brem_l.field(\"brem_photons_pe\").append(brem_e_l[itr])\n",
" #[:,\"brem_vtx_z\",:] brem vtx z\n",
" brem_l.field(\"brem_vtx_z\").append(brem_z_l[itr])\n",
" brem_l.end_record()\n",
"\n",
"brem_l = ak.Array(brem_l)\n",
"\n",
"\n",
"\n",
"\n",
"brem_f[0]"
]
},
{
"cell_type": "code",
2023-10-10 16:17:10 +02:00
"execution_count": 6,
2023-10-10 11:00:25 +02:00
"metadata": {},
"outputs": [],
"source": [
"acc_brem_found = ak.ArrayBuilder()\n",
"\n",
"for itr in range(ak.num(brem_f, axis=0)):\n",
" acc_brem_found.begin_record()\n",
" acc_brem_found.field(\"energy\").real(brem_f[itr,\"energy\"])\n",
" \n",
" acc_brem_found.field(\"brem_photons_pe\")\n",
" acc_brem_found.begin_list()\n",
" for jentry in range(brem_f[itr, \"photon_length\"]):\n",
" if brem_f[itr, \"brem_vtx_z\", jentry]>9500:\n",
" continue\n",
" else:\n",
" acc_brem_found.real(brem_f[itr,\"brem_photons_pe\", jentry])\n",
" \n",
" #acc_brem_found.field(\"brem_vtx_z\").real(brem_f[itr, \"brem_vtx_z\",jentry])\n",
" acc_brem_found.end_list()\n",
" \n",
" acc_brem_found.field(\"brem_vtx_z\")\n",
" acc_brem_found.begin_list()\n",
" for jentry in range(brem_f[itr, \"photon_length\"]):\n",
" if brem_f[itr, \"brem_vtx_z\", jentry]>9500:\n",
" continue\n",
" else:\n",
" acc_brem_found.real(brem_f[itr, \"brem_vtx_z\",jentry])\n",
" acc_brem_found.end_list()\n",
" \n",
"\n",
" \n",
" acc_brem_found.end_record()\n",
"\n",
"acc_brem_found = ak.Array(acc_brem_found)\n",
"\n",
"\n",
"\n",
"acc_brem_lost = ak.ArrayBuilder()\n",
"\n",
"for itr in range(ak.num(brem_l, axis=0)):\n",
" acc_brem_lost.begin_record()\n",
" acc_brem_lost.field(\"energy\").real(brem_l[itr,\"energy\"])\n",
" \n",
" acc_brem_lost.field(\"brem_photons_pe\")\n",
" acc_brem_lost.begin_list()\n",
" for jentry in range(brem_l[itr, \"photon_length\"]):\n",
" if brem_l[itr, \"brem_vtx_z\", jentry]>9500:\n",
" continue\n",
" else:\n",
" acc_brem_lost.real(brem_l[itr,\"brem_photons_pe\", jentry])\n",
" \n",
" #acc_brem_found.field(\"brem_vtx_z\").real(brem_f[itr, \"brem_vtx_z\",jentry])\n",
" acc_brem_lost.end_list()\n",
" \n",
" acc_brem_lost.field(\"brem_vtx_z\")\n",
" acc_brem_lost.begin_list()\n",
" for jentry in range(brem_l[itr, \"photon_length\"]):\n",
" if brem_l[itr, \"brem_vtx_z\", jentry]>9500:\n",
" continue\n",
" else:\n",
" acc_brem_lost.real(brem_l[itr, \"brem_vtx_z\",jentry])\n",
" acc_brem_lost.end_list()\n",
" \n",
" acc_brem_lost.end_record()\n",
"\n",
"acc_brem_lost = ak.Array(acc_brem_lost)\n"
]
},
{
"cell_type": "code",
2023-10-10 16:17:10 +02:00
"execution_count": 7,
2023-10-10 11:00:25 +02:00
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"cutoff = 50 MeV, sample size: 789\n",
"eff = 0.7427122940430925 +/- 0.0156\n",
"\n",
"cutoff = 100 MeV, sample size: 789\n",
"eff = 0.7427122940430925 +/- 0.0156\n",
"\n",
"cutoff = 150 MeV, sample size: 1025\n",
"eff = 0.7346341463414634 +/- 0.0138\n",
"\n",
"cutoff = 200 MeV, sample size: 1202\n",
"eff = 0.7346089850249584 +/- 0.0127\n",
"\n",
"cutoff = 250 MeV, sample size: 1337\n",
"eff = 0.7270007479431563 +/- 0.0122\n",
"\n",
"cutoff = 300 MeV, sample size: 1474\n",
"eff = 0.7164179104477612 +/- 0.0117\n",
"\n",
"cutoff = 350 MeV, sample size: 1599\n",
"eff = 0.7148217636022514 +/- 0.0113\n",
"\n",
"cutoff = 400 MeV, sample size: 1721\n",
"eff = 0.7094712376525276 +/- 0.0109\n",
"\n",
"cutoff = 450 MeV, sample size: 1854\n",
"eff = 0.7076591154261057 +/- 0.0106\n",
"\n",
"cutoff = 500 MeV, sample size: 1966\n",
"eff = 0.7004069175991862 +/- 0.0103\n",
"\n",
"cutoff = 550 MeV, sample size: 2065\n",
"eff = 0.6968523002421307 +/- 0.0101\n",
"\n",
"cutoff = 600 MeV, sample size: 2152\n",
"eff = 0.6909851301115242 +/- 0.01\n",
"\n",
"cutoff = 650 MeV, sample size: 2252\n",
"eff = 0.6873889875666075 +/- 0.0098\n",
"\n",
"cutoff = 700 MeV, sample size: 2346\n",
"eff = 0.6837169650468883 +/- 0.0096\n",
"\n",
"cutoff = 750 MeV, sample size: 2441\n",
"eff = 0.6780008193363376 +/- 0.0095\n",
"\n",
"cutoff = 800 MeV, sample size: 2534\n",
"eff = 0.6764009471191792 +/- 0.0093\n",
"\n",
"cutoff = 850 MeV, sample size: 2624\n",
"eff = 0.6714939024390244 +/- 0.0092\n",
"\n",
"cutoff = 900 MeV, sample size: 2707\n",
"eff = 0.6693756926486886 +/- 0.009\n",
"\n",
"cutoff = 950 MeV, sample size: 2785\n",
"eff = 0.6660682226211849 +/- 0.0089\n",
"\n",
"cutoff = 1000 MeV, sample size: 2839\n",
"eff = 0.6657273687918281 +/- 0.0089\n",
"\n",
"cutoff energy = 350MeV, sample size: 1599\n",
"eff = 0.7148 +/- 0.0113\n"
]
}
],
"source": [
"#finden wir die elektronen die keine bremsstrahlung gemacht haben mit hoher effizienz?\n",
"efficiencies_found = []\n",
"\n",
"\n",
"\n",
"for cutoff_energy in range(50,1050,50):\n",
"\tnobrem_f = acc_brem_found[ak.sum(acc_brem_found[\"brem_photons_pe\"],axis=-1,keepdims=False)<cutoff_energy]\n",
"\tnobrem_l = acc_brem_lost[ak.sum(acc_brem_lost[\"brem_photons_pe\"],axis=-1,keepdims=False)<cutoff_energy]\n",
"\n",
"\tif ak.num(nobrem_f,axis=0)+ak.num(nobrem_l,axis=0)==0:\n",
"\t\tcontinue\n",
"\teff = t_eff(nobrem_f,nobrem_l)\n",
"\tefficiencies_found.append(eff)\n",
"\n",
"\tprint(\"\\ncutoff = \",str(cutoff_energy) ,\"MeV, sample size: \",ak.num(nobrem_f,axis=0)+ak.num(nobrem_l,axis=0))\n",
"\tprint(\"eff = \",eff, \"+/-\", np.round(eff_err(nobrem_f, nobrem_l),4))\n",
"\n",
"\"\"\"\n",
"we see that a cutoff energy of xxxMeV is ideal because the efficiency drops significantly for higher values\n",
"\"\"\"\n",
"cutoff_energy = 350.0 #MeV\n",
"\n",
"\"\"\"\n",
"better statistics: cutoff=xxxMeV - sample size: xxx events and efficiency=xxxx\n",
"\"\"\"\n",
"nobrem_found = acc_brem_found[ak.sum(acc_brem_found[\"brem_photons_pe\"],axis=-1,keepdims=False)<cutoff_energy]\n",
"nobrem_lost = acc_brem_lost[ak.sum(acc_brem_lost[\"brem_photons_pe\"],axis=-1,keepdims=False)<cutoff_energy]\n",
"\n",
"print(\"\\ncutoff energy = 350MeV, sample size:\",ak.num(nobrem_found,axis=0)+ak.num(nobrem_lost,axis=0))\n",
"print(\"eff = \",np.round(t_eff(nobrem_found, nobrem_lost),4), \"+/-\", np.round(eff_err(nobrem_found, nobrem_lost),4))"
]
},
{
"cell_type": "code",
2023-10-10 16:17:10 +02:00
"execution_count": 8,
2023-09-25 16:36:13 +02:00
"metadata": {},
"outputs": [
{
"data": {
2023-10-10 11:00:25 +02:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAHICAYAAACyBMv/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABJbElEQVR4nO3deVxU9f4/8NcAw4CyuCCLirgvuIuaW5m7pmaZ5ZaSolcFLbU0/WZudUO7N691L5qaSmkulcttsZJMy61E1Ewx98J0EEUETYUB3r8//M1cx2FnZg6c83o+HjxqznzOmc/5zDi8+JzP53N0IiIgIiIiUgkXpStAREREZE8MN0RERKQqDDdERESkKgw3REREpCoMN0RERKQqDDdERESkKgw3REREpCoMN0RERKQqDDdERESkKgw3Cjt+/DiSkpKUrgYREZFqMNwobOXKlfj++++VrgYREZFqMNwo4Mcff8T48eORkpJi2ZaYmIjBgwfj6tWrCtaMiIio/GO4KaEdO3ZAp9NZftzc3FC1alU89thjePfdd5GZmZnvvmFhYQgICECLFi2wa9cu/Oc//0GvXr3Qt29fVKtWrVycg5bs2bPHqp0e/Pnpp59syh8/fhwRERGoV68ePD094enpiQYNGmDChAk4fPhwsV//6aefhqenJ27evJlvmZEjR0Kv1zs0HM+fPx86nQ7Xr1+32zEPHDiA+fPnF3hu5V1OTg78/f3xr3/9K98yjmjb8kDN7/+hQ4fQp08feHt7w8vLC926dcP+/futyhT3u+X27duYOnUqqlevDg8PD7Rq1QqbNm0qcTk1c1O6AuXVkSNHAABbt25FUFAQsrOzcfXqVezcuROvvPIKYmNjsXv3blSqVMlm34oVK2LhwoXIzc1FdHQ0XFxc8PHHH2PYsGHl5hy06K233kK3bt2stjVr1szq8YoVKzB58mQ0atQIL730Epo2bQqdTodTp05h48aNaNeuHc6dO4d69eoV+XUjIiKwfft2bNiwAZGRkTbPp6enY9u2bRgwYAACAgJKdnIKOXDgABYsWIAXXnhBtZ+zH3/8EdeuXcPgwYOVrkqZo9b3Pz4+Ho899hjat2+PdevWQUTw9ttvo0ePHti9ezc6duxoVb4o3y0AMHjwYMTHx2PRokVo2LAhNmzYgOHDhyM3NxcjRowodjlVEyqRp59+Wry8vCQ3N9fmudjYWAEgU6ZMyXPf48ePS8uWLaVPnz7yxBNPyLhx46RZs2byxBNPyLVr1xxddYvSnIOaXLlyRdLT0/N9fvfu3QJAPv300wKPs2/fPnFxcZGBAwdKZmZmnmU++eQTuXz5crHql52dLdWrV5ewsLA8n1++fLkAkC+++KJYxy2uefPmCQC7fkb/8Y9/CAC5ePGi3Y5Z1kRGRkrbtm0LLFOatv3rr79KWjXFldf3v7DvjD59+khAQIDVe5ORkSF+fn7SqVMny7aifreIiHz11VcCQDZs2GC1vVevXlK9enXJzs4uVjm1Y7gpoVq1akmXLl3yfb5mzZpSvXr1PJ+7dOmS5RdRVFSUrF27VrKzsyU2NlaysrJKXbclS5bItm3bCi1X0nPYu3ev9OrVS3x8fKRSpUryxBNPyJkzZ0pcrqgGDBggYWFhsnLlSmnRooV4eHhIzZo1Ze7cuZKTk1OsY924cUM++OAD6dGjh7i4uMjRo0fzLVvUL6AnnnhC9Hq9XLlypVh1OXPmjAwfPlyqVasm7u7u0rhxY/nPf/5jVWb27NkCQI4fP26zf/v27SUoKKhEX1rmX6pHjhyRp59+Wry9vcXHx0dGjhwpKSkpeZY9ceKEDBs2THx8fMTf31/GjBkjN2/etDn23r17pXv37uLl5SWenp7SsWNH+fLLL22O9/DP7t27i7R/ceuVkpIi48ePl5o1a4q7u7vlF01cXFy+7XPixAkBIJ988oll2+HDhwWAhIaGWpUdOHCgtGnTxmpbbm6uBAUFSXR0dL6v8eA5FPY+mMslJCTIM888I5UqVZLAwEDL80X5LJmP8csvv8iQIUPEx8dHKleuLNOmTROTySS//fab9OnTR7y8vCQkJEQWL15cYN0ftm3bNgEg3333nc1zy5Yts7x2fu//119/La1atZJ69epZvX9Go1ECAgKka9euRf6s53V8809xA1VxvjO8vLxk6NChNtsHDx4sACzfEcUJN+PGjRMvLy8xmUxW2zds2CAAZP/+/cUqJ1L6z4KjP0ulwXBTAtevXy+0V6N79+6i0+kK/aVrDjf2NGLECNHr9QUGnJKew7x588TFxUXGjh0rX331lXz22WfSvHlzCQ4Ollu3bhW7XHEEBQVJxYoVpUmTJrJu3TrZuXOnDBs2TADIqlWrCt3/r7/+kk2bNsmTTz4p7u7u4unpKc8884x8+umn+fa0iPzvC8jf319cXV3F29tbevfuLXv37rWUyc7OtvwCLo6TJ0+Kr6+vNG/eXD766CPZuXOnvPzyy+Li4iLz58+3lDt79qzodDqZOnWqzf4AZNasWcV6XTPzl1NISIjMmDFDvv32W1myZIlUrFhRWrdubRW2zWUbNWokc+fOlbi4OFmyZIkYDAYZM2aM1XH37Nkjer1ewsLCZPPmzbJ9+3bp3bu36HQ62bRpk4jcD/lTpkwRALJ161Y5ePCgHDx4UNLT04u0f3Hr1adPH6lWrZqsXLlS9uzZI9u3b5e5c+faHO9hQUFB8re//c3yeNGiReLp6SkALL1wJpNJfHx8ZObMmVb77tu3TwAUGuqL+j48WO7VV1+VuLg42b59u4gU/bP0YHu98cYbEhcXJzNnzhQAMnnyZGncuLG89957EhcXJ2PGjBEAsmXLlgLr/yCTyST+/v4ycuRIm+fat29vCYAFvf9nzpwRb29vGTx4sIiI5OTkSPfu3cXf379YfzyYj2n++f7776VGjRoSGBhYYM+LWUm/M9zd3WX06NE224cPHy4A5NtvvxWRon23mHXo0EHatWtns90cwFesWFGsciKl/yw4+rNUGgw3JbBz504BIGvWrMm3TOfOnaVixYpOrNX/ZGdnFxpwSnIOX3zxhQCQt99+26rcmTNnBICsX7++WOWK488//xQAUrduXau/5rKysiQwMFAGDBiQ535ZWVny5ZdfyogRI6RixYri7u4uAwYMkPXr1xc5ZB05ckReeukl2bZtm/z444+yZs0aadKkibi6uso333wjIiLJyckCQIYNG2azf3Z2tphMJsvPg5cB+/TpIzVr1rT5op08ebJ4eHjIjRs3LNu6du0qfn5+VoHj5ZdfLtIvz/yYv5ymTZtmtf3jjz+2ea/MZR9+XyMjI8XDw8PqvDp06CD+/v5WbZydnS3NmjWTmjVrWsrmd1miqPsXp15eXl424bAonn/+ealbt67lcc+ePWX8+PFSuXJl+fDDD0VEZP/+/QJAdu7cabXv1KlTpXnz5oW+RlHfB3O5uXPn2hyjqJ8l8zHeeecdq3KtWrWyBA0zk8kk1apVs4SMopo+fbp4enpa/VtNTEwUAPLvf//bsq2gy1KbN28WALJ06VKZO3euuLi42LRvcWRnZ8ugQYPEy8tLEhIS8i1nj++MVq1aScOGDa3+MDSZTFK3bl2rS0ZF+W4xa9CggfTp08fmta5cuSIA5K233ipWOZHSfxac8VkqKYabEoiOjrZ0IecnICDApovaHq5du1ZgV+vDP3q9XpKTk+1yDq1bt5Z69epJZmam1S9rk8kknp6esnDhwmKVKw5zV/dHH31k81yXLl2kQ4cONtt/+eUXqVKliri6ukrv3r1l9erVkpaWVuzXzktaWprUrFlTWrRoISIFh5uWLVtavSf/+Mc/RETk7t274ubmJlOmTLFppx07dggA2bFjh+U4H330kQCQzz77TETuf1kEBATIo48+WuLzMH85HT582Gq7yWQSNzc3iYiIsCn722+/WZV9//33BYDlc3b79m3R6XQSGRlp83qLFy8WAHLq1CkRyfuXW3H2L069unfvLpUqVZI33nhDDh48WORLwGvXrhUAcuHCBbl79654eHjIli1bZPDgwfL888+LiMiCBQvEYDD
2023-09-25 16:36:13 +02:00
"text/plain": [
2023-10-10 11:00:25 +02:00
"<Figure size 640x480 with 1 Axes>"
2023-09-25 16:36:13 +02:00
]
},
2023-10-10 11:00:25 +02:00
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"x_ = np.arange(50,1050,step=50)\n",
"plt.scatter(x_,efficiencies_found)\n",
"plt.xlabel(\"cutoff energy [MeV]\")\n",
"plt.ylabel(r\"$\\epsilon$\")\n",
"plt.title(r'$D^\\ast\\rightarrow Dee$, $p<5$GeV, photons w/ brem_vtx_z$<9500$mm')\n",
"plt.ylim([0.5,0.95])\n",
"plt.yticks(np.arange(0.5,0.96,step=0.01),minor=True)\n",
"plt.xticks(np.arange(0,1100,step=50),minor=True)\n",
"plt.grid()\n",
"plt.show()"
]
},
{
"cell_type": "code",
2023-10-10 16:17:10 +02:00
"execution_count": 9,
2023-10-10 11:00:25 +02:00
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"eff = 0.4993 +/- 0.0093\n"
]
},
{
"data": {
"text/html": [
"<pre>[{energy: 4.19e+03, brem_photons_pe: [1e+03, ..., 158], brem_vtx_z: [...]},\n",
" {energy: 4.12e+03, brem_photons_pe: [278, 468], brem_vtx_z: [...]},\n",
" {energy: 4.53e+03, brem_photons_pe: [445], brem_vtx_z: [9.33e+03]},\n",
" {energy: 4.22e+03, brem_photons_pe: [426, ..., 207], brem_vtx_z: [...]},\n",
" {energy: 2.54e+03, brem_photons_pe: [362, ...], brem_vtx_z: [...]},\n",
" {energy: 3.59e+03, brem_photons_pe: [715], brem_vtx_z: [2.59e+03]},\n",
" {energy: 3.78e+03, brem_photons_pe: [198, ..., 180], brem_vtx_z: [...]},\n",
" {energy: 3.07e+03, brem_photons_pe: [1.26e+03, ...], brem_vtx_z: [...]},\n",
" {energy: 3.89e+03, brem_photons_pe: [486], brem_vtx_z: [5.2e+03]},\n",
" {energy: 4.66e+03, brem_photons_pe: [906, ..., 538], brem_vtx_z: [...]},\n",
" ...,\n",
" {energy: 3.08e+03, brem_photons_pe: [575, 337], brem_vtx_z: [...]},\n",
" {energy: 2.39e+03, brem_photons_pe: [163, 971], brem_vtx_z: [...]},\n",
" {energy: 2.57e+03, brem_photons_pe: [184, 798], brem_vtx_z: [...]},\n",
" {energy: 2.19e+03, brem_photons_pe: [161, 337], brem_vtx_z: [...]},\n",
" {energy: 3.62e+03, brem_photons_pe: [1.19e+03], brem_vtx_z: [3.98e+03]},\n",
" {energy: 3.96e+03, brem_photons_pe: [135, ...], brem_vtx_z: [...]},\n",
" {energy: 3.48e+03, brem_photons_pe: [3.17e+03], brem_vtx_z: [8.58e+03]},\n",
" {energy: 3.19e+03, brem_photons_pe: [240, ...], brem_vtx_z: [...]},\n",
" {energy: 3.6e+03, brem_photons_pe: [142, ...], brem_vtx_z: [...]}]\n",
"---------------------------------------------------------------------------\n",
"type: 1452 * {\n",
" energy: float64,\n",
" brem_photons_pe: var * float64,\n",
" brem_vtx_z: var * float64\n",
"}</pre>"
],
"text/plain": [
"<Array [{energy: 4.19e+03, ...}, ..., {...}] type='1452 * {energy: float64,...'>"
]
},
2023-10-10 16:17:10 +02:00
"execution_count": 9,
2023-09-25 16:36:13 +02:00
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#wie viel energie relativ zur anfangsenergie verlieren die elektronen durch bremstrahlung und hat das einen einfluss darauf ob wir sie finden oder nicht?\n",
2023-10-10 11:00:25 +02:00
"\n",
"cutoff_energy=350\n",
"\n",
"brem_found = acc_brem_found[ak.sum(acc_brem_found[\"brem_photons_pe\"],axis=-1,keepdims=False)>=cutoff_energy]\n",
2023-09-25 16:36:13 +02:00
"energy_found = ak.to_numpy(brem_found[\"energy\"])\n",
"eph_found = ak.to_numpy(ak.sum(brem_found[\"brem_photons_pe\"], axis=-1, keepdims=False))\n",
2023-10-10 11:00:25 +02:00
"residual_found = energy_found - eph_found\n",
2023-09-25 16:36:13 +02:00
"energyloss_found = eph_found/energy_found\n",
"\n",
2023-10-10 11:00:25 +02:00
"brem_lost = acc_brem_lost[ak.sum(acc_brem_lost[\"brem_photons_pe\"],axis=-1,keepdims=False)>=cutoff_energy]\n",
2023-09-25 16:36:13 +02:00
"energy_lost = ak.to_numpy(brem_lost[\"energy\"])\n",
"eph_lost = ak.to_numpy(ak.sum(brem_lost[\"brem_photons_pe\"], axis=-1, keepdims=False))\n",
2023-10-10 11:00:25 +02:00
"residual_lost = energy_lost - eph_lost\n",
2023-09-25 16:36:13 +02:00
"energyloss_lost = eph_lost/energy_lost\n",
"\n",
2023-10-10 11:00:25 +02:00
"print(\"eff = \", np.round(t_eff(brem_found,brem_lost),4), \"+/-\", np.round(eff_err(brem_found, brem_lost),4))\n",
"brem_lost"
2023-09-25 16:36:13 +02:00
]
},
{
"cell_type": "code",
2023-10-10 11:00:25 +02:00
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
2023-10-10 16:17:10 +02:00
"execution_count": 10,
2023-09-25 16:36:13 +02:00
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
2023-10-10 11:00:25 +02:00
"mean energyloss relative to initial energy (found): 0.3135\n",
"mean energyloss relative to initial energy (lost): 0.4443\n"
2023-09-25 16:36:13 +02:00
]
}
],
"source": [
"mean_energyloss_found = ak.mean(energyloss_found)\n",
"mean_energyloss_lost = ak.mean(energyloss_lost)\n",
2023-10-10 11:00:25 +02:00
"print(\"mean energyloss relative to initial energy (found): \", np.round(mean_energyloss_found,4))\n",
"print(\"mean energyloss relative to initial energy (lost): \", np.round(mean_energyloss_lost,4))"
2023-09-25 16:36:13 +02:00
]
},
{
"cell_type": "code",
2023-10-10 16:17:10 +02:00
"execution_count": 11,
2023-09-25 16:36:13 +02:00
"metadata": {},
"outputs": [
{
"data": {
2023-10-10 11:00:25 +02:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAHMCAYAAAD7xYOxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABLmUlEQVR4nO3deVwU9f8H8Ncut9wqIIpHKKLgVV6h5pV3Hh1WHhl5lWkeaWpWJpZ91fpZaWmmqXR5lFdm5lGJt3lnhiaoJeaBigh4ICzv3x/G5rIL7A67zACv5+OxD52Zz3zmPfPZ483M5zOjExEBERERkQbp1Q6AiIiIKD9MVIiIiEizmKgQERGRZjFRISIiIs1iokJERESaxUSFiIiINIuJChEREWkWExUiIiLSLCYqREREpFlMVOzo6NGjOHv2rNphEBERlRpMVOxowYIF+OWXX9QOg4iIqNRgolJE27dvx9ChQ5GcnGycFx8fj8cffxyXLl1SMTIiIqKSj4kKgA0bNkCn0xlfzs7OqFChAlq3bo3Zs2cjMzMz33UbN26MoKAgNGjQAD///DM+/vhjdOzYEV26dEFAQECJ2IeyJC4uzuQ43fvau3evWfmjR49i8ODBqFmzJjw8PODh4YGwsDC88MILOHDggM3bf+yxx+Dh4YHU1NR8y/Tv3x8uLi4OTXRjYmKg0+lw5coVu9W5e/duxMTEFLhvJZ3BYEBgYCA++OCDfMs44tiWBKW5/fft24fOnTvD29sbXl5eaNeuHXbt2mVSxtbvloyMDIwZMwaVK1eGu7s7GjVqhOXLlysuV5o5qx2AFhw6dAgAsHr1agQHByM7OxuXLl3C5s2b8corryA2NhZbt26Fn5+f2bqenp546623kJOTg+nTp0Ov1+Prr79Gnz59Ssw+lEX/+9//0K5dO5N59erVM5n+9NNP8dJLLyE8PByjR49GZGQkdDodjh8/jmXLlqFp06ZITExEzZo1rd7u4MGDsXbtWixduhTDhw83W379+nWsWbMG3bt3R1BQkLKdU8nu3bsxdepUPPfcc6X2fbZ9+3ZcvnwZjz/+uNqhaE5pbf/9+/ejdevWaNasGb788kuICN599108/PDD2Lp1K6KiokzKW/PdAgCPP/449u/fjxkzZqB27dpYunQp+vbti5ycHPTr18/mcqWakDz22GPi5eUlOTk5ZstiY2MFgIwcOdLiukePHpWGDRtK586dpVu3bjJkyBCpV6+edOvWTS5fvuzo0I2Ksg+lyfnz5+X69ev5Lt+6dasAkG+//bbAenbu3Cl6vV569OghmZmZFst888038s8//9gUX3Z2tlSuXFkaN25scfknn3wiAOT777+3qV5bTZkyRQDY9T363nvvCQA5c+aM3erUmuHDh0uTJk0KLFOUY3vjxg2loamupLZ/Yd8ZnTt3lqCgIJO2SUtLk4oVK0qLFi2M86z9bhER+eGHHwSALF261GR+x44dpXLlypKdnW1TudKOiYqIVKtWTVq1apXv8pCQEKlcubLFZUlJScYflREjRsiSJUskOztbYmNj5c6dO0WO7f3335c1a9YUWk7pPuzYsUM6duwoPj4+4ufnJ926dZOTJ08qLmet7t27S+PGjWXBggXSoEEDcXd3l5CQEHnzzTfFYDDYVFdKSop89tln8vDDD4ter5fDhw/nW9baL5Nu3bqJi4uLnD9/3qZYTp48KX379pWAgABxdXWVOnXqyMcff2xSZtKkSQJAjh49arZ+s2bNJDg4WNEXUO4P5KFDh+Sxxx4Tb29v8fHxkf79+0tycrLFsseOHZM+ffqIj4+PBAYGysCBAyU1NdWs7h07dkj79u3Fy8tLPDw8JCoqStavX29WX97X1q1brVrf1riSk5Nl6NChEhISIq6ursYfjS1btuR7fI4dOyYA5JtvvjHOO3DggACQiIgIk7I9evSQBx54wGReTk6OBAcHy/Tp0/Pdxr37UFg75JY7ePCgPPHEE+Ln5yeVKlUyLrfmvZRbx2+//Sa9e/cWHx8f8ff3l5dfflmysrLkxIkT0rlzZ/Hy8pLq1avLzJkzC4w9rzVr1ggA+emnn8yWzZs3z7jt/Nr/xx9/lEaNGknNmjVN2u/ChQsSFBQkbdq0sfq9bqn+3JetyZEt3xleXl7y9NNPm81//PHHBYDxO8KWRGXIkCHi5eUlWVlZJvOXLl0qAGTXrl02lRMp+nvB0e+loijzicqVK1cKPdvQvn170el0hf6A5iYq9tSvXz9xcXEpMFlRug9TpkwRvV4vgwYNkh9++EFWrlwp9evXl6pVq0p6errN5WwRHBwsnp6eUrduXfnyyy9l8+bN0qdPHwEgCxcuLHT9GzduyPLly6Vnz57i6uoqHh4e8sQTT8i3336b7xkQkf++TAIDA8XJyUm8vb2lU6dOsmPHDmOZ7Oxs44+pLf744w/x9fWV+vXryxdffCGbN2+WcePGiV6vl5iYGGO5hIQE0el0MmbMGLP1Acirr75q03Zz5X7RVK9eXcaPHy+bNm2S999/Xzw9PeX+++83SZxzy4aHh8ubb74pW7Zskffff1/c3Nxk4MCBJvXGxcWJi4uLNG7cWFasWCFr166VTp06iU6nk+XLl4vI3YR95MiRAkBWr14te/bskT179sj169etWt/WuDp37iwBAQGyYMECiYuLk7Vr18qbb75pVl9ewcHB8vzzzxunZ8yYIR4eHgLAeHYsKytLfHx8ZMKECSbr7ty5UwAUmqBb2w73lps4caJs2bJF1q5dKyLWv5fuPV5vv/22bNmyRSZMmCAA5KWXXpI6derInDlzZMuWLTJw4EABIKtWrSow/ntlZWVJYGCg9O/f32xZs2bNjMlcQe1/8uRJ8fb2lscff1xERAwGg7Rv314CAwNt+kMgt87c1y+//CJVqlSRSpUqFXhGJJfS7wxXV1d59tlnzeb37dtXAMimTZtExLrvllwPPvigNG3a1Gx+bjL96aef2lROpOjvBUe/l4qizCcqmzdvFgCyePHifMu0bNlSPD09izGq/2RnZxearCjZh++//14AyLvvvmtS7uTJkwJAvvrqK5vK2eLcuXMCQEJDQ03+yrpz545UqlRJunfvbnG9O3fuyPr166Vfv37i6ekprq6u0r17d/nqq6+sTpgOHToko0ePljVr1sj27dtl8eLFUrduXXFycpKNGzeKiMjFixcFgPTp08ds/ezsbMnKyjK+7r3U1rlzZwkJCTH70nzppZfE3d1dUlJSjPPatGkjFStWNEkexo0bZ9UPYX5yv2hefvllk/lff/21WVvlls3brsOHDxd3d3eT/XrwwQclMDDQ5BhnZ2dLvXr1JCQkxFg2v1P/1q5vS1xeXl5miZ41nnnmGQkNDTVOd+jQQYYOHSr+/v7y+eefi4jIrl27BIBs3rzZZN0xY8ZI/fr1C92Gte2QW+7NN980q8Pa91JuHbNmzTIp16hRI2PSkCsrK0sCAgKMCYO1xo4dKx4eHiaf1fj4eAEgH330kXFeQZd+VqxYIQDkww8/lDfffFP0er3Z8bVFdna29OrVS7y8vOTgwYP5lrPHd0ajRo2kdu3aJn/kZWVlSWhoqMllGWu+W3KFhYVJ586dzbZ1/vx5ASD/+9//bConUvT3QnG8l5Qq84nK9OnTjadp8xMUFGR2GtgeLl++XODpzLwvFxcXuXjxol324f7775eaNWtKZmamyQ9vVlaWeHh4yFtvvWVTOVvknk7+4osvzJa1atVKHnzwQbP5v/32m5QvX16cnJykU6dOsmjRIrl27ZrN27bk2rVrEhISIg0aNBCRghOVhg0bmrTJe++9JyIit27dEmdnZxk5cqTZcdqwYYMAkA0bNhjr+eKLLwSArFy5UkTufvCDgoLkoYceUrwfuV80Bw4cMJmflZUlzs7OMnjwYLOyJ06cMCk7f/58AWB8n2VkZIhOp5Phw4ebbW/mzJkCQI4fPy4iln+obFnflrjat28vfn5+8vbbb8uePXusvsy6ZMkSASCnT5+WW7duibu7u6xatUoef/xxeeaZZ0REZOrUqeLm5iY3b940WbdatWomZzPyY20
2023-09-25 16:36:13 +02:00
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
2023-10-10 11:00:25 +02:00
"plt.hist(energyloss_lost, bins=100, density=True, alpha=0.5, histtype='bar', color=\"darkorange\", label=\"lost\")\n",
"plt.hist(energyloss_found, bins=100, density=True, alpha=0.5, histtype='bar', color=\"blue\", label=\"found\")\n",
2023-09-25 16:36:13 +02:00
"plt.xticks(np.arange(0,1.1,0.1), minor=True,)\n",
2023-10-10 11:00:25 +02:00
"plt.yticks(np.arange(0,5.5,0.5), minor=True)\n",
2023-09-25 16:36:13 +02:00
"plt.xlabel(r\"$E_\\gamma/E_0$\")\n",
"plt.ylabel(\"counts (normed)\")\n",
2023-10-10 11:00:25 +02:00
"plt.title(r'$D^\\ast\\rightarrow Dee$, $p<5$GeV, photons w/ brem_vtx_z$<9500$mm')\n",
2023-09-25 16:36:13 +02:00
"plt.legend()\n",
2023-10-10 11:00:25 +02:00
"plt.grid(alpha=0.5)\n",
2023-09-25 16:36:13 +02:00
"\n",
"\"\"\"\n",
"found: elektronen verlieren durchschnittlich 0.44 ihrer anfangsenergie durch bremsstrahlung\n",
"lost: elektronen verlieren durchschnittlich 0.59 ihrer anfangsenergie durch bremsstrahlung\n",
"\n",
"-> lost electrons lose slightly more energy than found electrons. This is however nowhere near as extreme as for the B decay\n",
"\"\"\"\n",
"\n",
2023-09-28 12:12:21 +02:00
"#check if pt near 1 have small eta; perhaps exclude pt with eta outside range; or 2d plot\n",
"\n",
2023-09-25 16:36:13 +02:00
"plt.show()"
]
},
{
"cell_type": "code",
2023-10-10 16:17:10 +02:00
"execution_count": 12,
2023-10-10 11:00:25 +02:00
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABiUAAAJOCAYAAADYhwTwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAACkVElEQVR4nOzdd3hUVeL/8c8lCWkkIaEk9CYgUgRBQ3ElobNEVBSUIIKVFRBRUOwJFhBchRVErOACAXRV1kYElGBBiiAWcLGhgBKKCaFIS3J+f/jN/BySkLk3mcmQvF/PM8+T3Dn9ljlnzpx7LWOMEQAAAAAAAAAAgJdVKe8CAAAAAAAAAACAyoFJCQAAAAAAAAAA4BNMSgAAAAAAAAAAAJ9gUgIAAAAAAAAAAPgEkxIAAAAAAAAAAMAnmJQAAAAAAAAAAAA+waQEAAAAAAAAAADwCSYlAAAAAAAAAACATzApAQAAAAAAAAAAfIJJCQAAAD/x1VdfaefOneVdDAAAAAAAvIZJCQAAAD/x/PPP68MPPyzvYgAAAAAA4DVMSgAAAJSjjz76SDfffLP27dvn2rZt2zYNGjRIe/fuLceSAQAAAABQ9piUAAAAKKX33ntPlmW5XoGBgapRo4YuueQS/etf/9KJEyeKjduxY0fFxsaqXbt2+uCDDzR79mz17t1b/fr1U61atc6KOlQmGRkZbu3019e6desKhf/qq6904403qlmzZgoNDVVoaKiaN2+uUaNG6fPPP7ed/xVXXKHQ0FAdPHiw2DDDhg1TUFCQVye1UlNTZVmWDhw4UGZprl27VqmpqWes29kuLy9PtWvX1owZM4oN4422PRtU5P2/YcMG9e3bVxEREapWrZoSExP16aefuoWxe205cuSIxo8fr7p16yokJETt27fXkiVLHIcDAADwpcDyLgAAAMDZbvPmzZKkN954Q3Xq1FFubq727t2rFStWaOLEiZo/f75Wr16t6tWrF4obHh6uhx9+WPn5+Zo6daqqVKmiRYsW6Zprrjlr6lAZTZkyRYmJiW7b2rRp4/b/c889p7Fjx6ply5a6/fbb1bp1a1mWpW+//VaLFy/WhRdeqB9++EHNmjXzON8bb7xRy5YtU1pamkaPHl3o/ZycHL355ptKSkpSbGyss8qVk7Vr12ry5MkaOXJkhT3OPvroI+3fv1+DBg0q76L4nYq6/zdu3KhLLrlEF110kRYsWCBjjKZPn66ePXtq9erV6tKli1t4T64tkjRo0CBt3LhRjz/+uFq0aKG0tDQNHTpU+fn5Sk5Oth0OAADAl5iUAAAAKKXNmzerWrVquvzyy2VZlmv7lVdeqa5du2rkyJF66KGH9PTTTxeK+/XXX2v48OGKi4vT3//+d9WtW1ePPfaYFixYoFdeeUU1a9b0+zpUJHv27FF4eLgiIyPPGK558+bq3Llzse9/+umnGj16tAYMGKD//Oc/qlq1quu9Hj16aMyYMXrttdcUGhpqq3z9+/dX3bp19fLLLxc5KbF48WIdO3ZMN954o6104Rv/+c9/1KlTJzVq1Mgr6f/xxx8KCwvzStooWknXjAcffFDVq1dXenq6a9/06tVLTZs21cSJEwutmCjp2iL9ubJt5cqVrgkGSUpMTNQvv/yiu+66S1dffbUCAgI8DgcAAOBr3L4JAACglDZt2qT27du7fZlfYMSIEapfv75ef/31IuNGR0fr0UcfVXp6upo0aaJu3bppy5YtGjJkiKKiokpdthkzZmjZsmUlhnNah08++UR9+vRRVFSUoqOjNWDAAH3//feOw3nq0ksvVadOnfTCCy/o/PPPV2hoqBo0aKCUlBTl5+fbSis7O1svvfSSevXqpfr16+unn35yXK4CU6ZMUUBAgJ577jm3CYm/Gjx4sOrWreu27fvvv1dycrJq166t4OBgtWrVSs8884zr/YCAAI0YMUKbNm3S119/XSjNefPmqU6dOurfv7/tMhfcNuiLL77QoEGDFBkZqaioKF177bXav39/kXH27t2roUOHKioqSrGxsbrhhhuUk5NTKNwnn3yinj17KiIiQmFhYerataveffddt7zvuusuSVKTJk1ct63JyMjwKP7pddi6desZy7V//37dcsstatCggYKDg1WrVi1169ZNq1atKrZ9tm7dKsuy9Nprr7m2bdq0SZZlqXXr1m5hBw4cqI4dO7ptM8bozTff1JVXXllsHn+1a9euM+6Hgrpu3rxZV111laKjo91W3ZR0LP01ja+++kqDBw9WVFSUYmJidOeddyo3N1fbt29Xv379FBERocaNG2v69Okelb3AsmXLZFmWPvjgg0LvPfvss668i9v/6enp6tChg8455xy3/ZeZmam4uDglJCQoLy/Po7IUd2sky7L0888/26qXnWvGp59+qoSEBLfJooiICF1yySVau3at9uzZYytvSXrzzTdVrVo1DR482G379ddfr99++03r16+3FU4q/bHg7WMJAABULExKAAAAlMLvv/+unTt3qkOHDsWGadGihfbs2VPkl+X169dXUlKS27aCL56DgoJKXb7PP/9cQ4YMOePEhNM6pKamqnv37mrQoIEWL16sF198Ubt27VLPnj115MgR2+Hs2LRpk/73v/9pxowZuuuuu/TWW2/p4osv1sMPP6yXX365xPh//PGHli5dqssuu0xxcXG67bbbVL16dS1dulTnnXdeifHHjBmjwMBARUZGqm/fvvrkk09c7+Xl5Wn16tXq1KmT6tSp43Gdtm3bpgsvvFDffPONnnzySb3zzjsaMGCAxo0bp8mTJ7vC3XDDDbIsq1A9t23bpg0bNmjEiBGl+vXzFVdcoXPOOUf/+c9/lJqaqmXLlqlv3746depUobBXXnmlWrRooddff1333HOP0tLSdMcdd7iFWbNmjXr06KGcnBy99NJLWrx4sSIiInTppZdq6dKlkqSbbrpJt912m6Q/byH22Wef6bPPPtMFF1zgUXy75Ro+fLiWLVumhx56SCtWrNCLL76oXr166ffffy+2XVq3bq06deq4TVysWrVKoaGh2rZtm3777TdJUm5urtasWaNevXq5xS/4AtrTSQlP98OgQYN0zjnn6LXXXtPcuXMleX4sFRgyZIjOP/98vf7667r55ps1Y8YM3XHHHbr88ss1YMAAvfnmm+rRo4cmTZqkN954w6PyS1JSUpJq166tefPmFXpv/vz5uuCCC9SuXbti93/Xrl316quvat++fbrhhhskSfn5+Ro2bJiMMVq8eLHHx3pBmgWvDz/8UPXq1VNcXJxiYmJKjO/0mnHy5EkFBwcX2l6w7fTJxTNdWwp88803atWqlQID3W980K5dO9f7dsL9VWmPBW8dSwAAoIIxAAAAcGzFihVGknn55ZeLDdOtWzcTHh7uw1L9f7m5uSY5OdkEBQWZN998s8gwTurw9ttvG0lm+vTpbuG+++47I8ksXLjQVjg7du/ebSSZpk2bmoMHD7q2nzx50sTFxZmkpKQi4508edK88847Jjk52YSHh5uqVauapKQks3DhQnP48GGP8t68ebO5/fbbzZtvvmk++ugj8/LLL5tWrVqZgIAAk56ebowxJjMz00gy11xzTaH4ubm55tSpU65Xfn6+672+ffua+vXrm5ycHLc4Y8eONSEhISYrK8u1rXv37qZmzZrm5MmTrm0TJkwwksx3333nUV1Ol5KSYiSZO+64w237okWLCu2rgrCn79fRo0ebkJAQt3p17tzZ1K5d262Nc3NzTZs2bUz9+vVdYZ944gkjyezYscMtTU/j2ylXtWrVzPjx4+00jzHGmGuvvdY0bdrU9X+vXr3MzTffbKKjo80rr7xijDHm008/NZLMihUr3OKOHz/etG3btsQ8PN0PBeEeeuihQml4eiwVpPHkk0+6hWvfvr2RZN544w3XtlOnTplatWqZQYMGlViHv7rzzjtNaGio27m6bds2I8nMmjXLta24/W+MMUuXLjWSzMyZM81DDz1kqlSpUqh97cjNzTWXXXaZqVatmtm0aVOx4crimtG+fXvTokULk5eX59p26tQp07RpUyPJpKWlGWM8u7YUaN6
"text/plain": [
"<Figure size 2000x600 with 3 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"#energyloss in abh von der energie der elektronen\n",
"\n",
"fig, ((ax0, ax1)) = plt.subplots(nrows=1, ncols=2, figsize=(20,6))\n",
"\n",
"a0=ax0.hist2d(energyloss_found, energy_found, bins=(np.linspace(0,1,50), np.linspace(0,6e3,50)), cmap=plt.cm.jet, cmin=1, vmax=8)\n",
"ax0.set_ylim(0,6e3)\n",
"ax0.set_xlim(0,1)\n",
"ax0.set_xlabel(r\"energyloss $E_\\gamma/E_0$\")\n",
"ax0.set_ylabel(r\"$E_0$\")\n",
"ax0.set_title(\"found energyloss wrt electron energy\")\n",
"\n",
"a1=ax1.hist2d(energyloss_lost, energy_lost, bins=(np.linspace(0,1,50), np.linspace(0,6e3,50)), cmap=plt.cm.jet, cmin=1, vmax=8)\n",
"ax1.set_ylim(0,6e3)\n",
"ax1.set_xlim(0,1)\n",
"ax1.set_xlabel(r\"energyloss $E_\\gamma/E_0$\")\n",
"ax1.set_ylabel(r\"$E_0$\")\n",
"ax1.set_title(\"lost energyloss wrt electron energy\")\n",
"\n",
"fig.colorbar(a1[3],ax=ax1)\n",
"fig.suptitle(r\"$D^\\ast\\rightarrow D ee$, $p<5$GeV, photons w/ brem_vtx_z$<9500$mm\")\n",
"\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
2023-10-10 16:17:10 +02:00
"execution_count": 13,
2023-09-25 16:36:13 +02:00
"metadata": {},
"outputs": [],
"source": [
"#ist die shape der teilspur im scifi anders? (koenntest du zum beispiel durch vergleich der verteilungen der fit parameter studieren,\n",
"#in meiner thesis findest du das fitmodell -- ist einfach ein polynom dritten grades)\n",
"z_ref=8520 #mm\n",
"\n",
"def scifi_track(z, a, b, c, d):\n",
" return a + b*(z-z_ref) + c*(z-z_ref)**2 + d*(z-z_ref)**3\n",
"\n",
"def z_mag(xv, zv, tx, a, b):\n",
" \"\"\" optical centre of the magnet is defined as the intersection between the trajectory tangents before and after the magnet\n",
"\n",
" Args:\n",
" xv (double): velo x track\n",
" zv (double): velo z track\n",
" tx (double): velo x slope\n",
" a (double): ax parameter of track fit\n",
" b (double): bx parameter of track fit\n",
"\n",
" Returns:\n",
" double: z_mag\n",
" \"\"\"\n",
" return (xv-tx*zv-a+b*z_ref)/(b-tx)"
]
},
{
"cell_type": "code",
2023-10-10 16:17:10 +02:00
"execution_count": 14,
2023-09-25 16:36:13 +02:00
"metadata": {},
"outputs": [],
"source": [
"scifi_found = found[found[\"scifi_hit_pos_x_length\"]>3]\n",
"scifi_lost = lost[lost[\"scifi_hit_pos_x_length\"]>3]\n",
"\n",
"scifi_x_found = scifi_found[\"scifi_hit_pos_x\"]\n",
"scifi_z_found = scifi_found[\"scifi_hit_pos_z\"]\n",
"\n",
"tx_found = scifi_found[\"velo_track_tx\"]\n",
"\n",
"scifi_x_lost = scifi_lost[\"scifi_hit_pos_x\"]\n",
"scifi_z_lost = scifi_lost[\"scifi_hit_pos_z\"]\n",
"\n",
"tx_lost = scifi_lost[\"velo_track_tx\"]\n",
"\n",
"xv_found = scifi_found[\"velo_track_x\"]\n",
"zv_found = scifi_found[\"velo_track_z\"]\n",
"\n",
"xv_lost = scifi_lost[\"velo_track_x\"]\n",
"zv_lost = scifi_lost[\"velo_track_z\"]\n",
"\n",
"\n",
"\n",
"#ak.num(scifi_found[\"energy\"], axis=0)\n",
"#scifi_found.snapshot()"
]
},
{
"cell_type": "code",
2023-10-10 16:17:10 +02:00
"execution_count": 15,
2023-09-25 16:36:13 +02:00
"metadata": {},
"outputs": [],
"source": [
"scifi_fitpars_found = ak.ArrayBuilder()\n",
"\n",
"for i in range(0,ak.num(scifi_found[\"energy\"], axis=0)):\n",
" popt, pcov = curve_fit(scifi_track,ak.to_numpy(scifi_z_found[i,:]),ak.to_numpy(scifi_x_found[i,:]))\n",
" scifi_fitpars_found.begin_list()\n",
" scifi_fitpars_found.real(popt[0])\n",
" scifi_fitpars_found.real(popt[1])\n",
" scifi_fitpars_found.real(popt[2])\n",
" scifi_fitpars_found.real(popt[3])\n",
" scifi_fitpars_found.end_list()\n",
"\n",
"scifi_fitpars_lost = ak.ArrayBuilder()\n",
"\n",
"for i in range(0,ak.num(scifi_lost[\"energy\"], axis=0)):\n",
" popt, pcov = curve_fit(scifi_track,ak.to_numpy(scifi_z_lost[i,:]),ak.to_numpy(scifi_x_lost[i,:]))\n",
" scifi_fitpars_lost.begin_list()\n",
" scifi_fitpars_lost.real(popt[0])\n",
" scifi_fitpars_lost.real(popt[1])\n",
" scifi_fitpars_lost.real(popt[2])\n",
" scifi_fitpars_lost.real(popt[3])\n",
" scifi_fitpars_lost.end_list()\n",
"\n",
"\n",
"scifi_fitpars_lost = scifi_fitpars_lost.to_numpy()\n",
"scifi_fitpars_found = scifi_fitpars_found.to_numpy()\n",
"\n",
"\n",
"\n",
"dtx_found = scifi_fitpars_found[:,1] - tx_found\n",
"dtx_lost = scifi_fitpars_lost[:,1] - tx_lost\n"
]
},
{
"cell_type": "code",
2023-10-10 16:17:10 +02:00
"execution_count": 16,
2023-09-25 16:36:13 +02:00
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABOwAAANVCAYAAADLJ2G2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAADqDElEQVR4nOzdeVzU1f7H8ffIjguKC4gLkmWuWUkq5IKZGKlZ6dX0XtxNMzOlm7kmWmmrl8zUFotsUW+praZSiUuQpaJlaouplIGK5q6s398f/ZjrMAMM6wzwej4e86g5c77n+/l+vzPj4TPne47JMAxDAAAAAAAAAJxCNUcHAAAAAAAAAOB/SNgBAAAAAAAAToSEHQAAAAAAAOBESNgBAAAAAAAAToSEHQAAAAAAAOBESNgBAAAAAAAAToSEHQAAAAAAAOBESNgBAAAAAAAAToSEHQAAAAAAAOBESNgBFcDq1avVpk0beXl5yWQyac+ePYqOjpbJZLKot2TJEsXGxjomSCf33nvvKSYmxtFhlFh8fLxMJpPi4+MdHYrdYmNjZTKZdOTIEYvyWbNmqWnTpnJ1dVXt2rUlSWFhYQoLCyu0TWd4r5tMJk2cONGhMQAAKj76eSXniH5eWFiY2rZtW+rtNmvWTCNGjCj1dsuSyWRSdHS0RdmXX36p4OBgVa9eXSaTSR9++GG+fcK8EhISFB0drTNnzpRZzIUZMWKEatSo4bD9AxIJO8DpnTx5UpGRkWrevLk2bNigxMREtWjRQmPGjFFiYqJFXTpy+assCbuKqE+fPkpMTFTDhg3NZR999JGeeuopDRs2TFu2bNEXX3wh6e/38JIlSwptk/c6AKAyoJ9XOujnOVZiYqLGjBljfm4YhgYNGiQ3Nzd9/PHHSkxMVPfu3W32CW1JSEjQ3LlzHZqwA5yBq6MDAFCwn3/+WZmZmfrXv/6l7t27m8u9vb3VuHFjB0Zm7fLly/L09LT6Rbgyu3z5sry8vBwdRrFdunRJ3t7eZbqP+vXrq379+hZl+/btkyRNmjRJDRo0MJe3bt261PefmZkpk8kkV1f+yQMAOBf6ec6tovfzsrOzlZWVJQ8PjzLdT+fOnS2e//nnnzp9+rTuuece9ezZ0+K1vH3C0lDRrxOQH0bYAU5sxIgR6tKliyRp8ODBMplM5tsF894q0axZM/3444/asmWLTCaTTCaTmjVrJul/t1G+8847ioqKkr+/v7y8vNS9e3clJSVZ7HPnzp2677771KxZM3l5ealZs2YaMmSIjh49alEvd0j7pk2bNGrUKNWvX1/e3t5KT0/Xr7/+qpEjR+q6666Tt7e3GjVqpH79+umHH36waCM3rvfee0+PPfaYGjZsqBo1aqhfv346fvy4zp8/r/vvv1/16tVTvXr1NHLkSF24cMGiDcMwtGTJEt14443y8vJSnTp1NHDgQP3222/mOmFhYfrss8909OhR87m5+txlZGToySefVMuWLeXh4aH69etr5MiROnnypMW+mjVrpr59+2rt2rW66aab5Onpqblz5+Z7/eLi4tS/f381btxYnp6euvbaazVu3DilpaXlu83VDh48qDvuuEPe3t6qV6+exo8fr/Pnz9us+8UXX6hnz56qVauWvL29deutt+rLL7+0qJP7ntm9e7cGDhyoOnXqqHnz5vnu/9KlS/r3v/+toKAgeXp6ytfXV8HBwVq5cqVFvR07dqhfv36qW7euPD091bx5c02ePNn8et7bH5o1a6ZZs2ZJkvz8/Cxuo7Dnllh73utvv/22HnnkETVq1EgeHh769ddfdfLkSU2YMEGtW7dWjRo11KBBA912223atm2b1T7S09M1b948tWrVSp6enqpbt6569OihhISEfOMyDEMzZsyQm5ubXnvttQKPAQAA+nkVu5+Xa9u2bercubO8vLzUqFEjzZ49W9nZ2YVul5mZqalTp8rf31/e3t7q0qWLvv32W5t1U1NTNW7cODVu3Fju7u4KCgrS3LlzlZWVZa5z5MgRmUwmPfvss3ryyScVFBQkDw8Pbd68Od8Y3n//fXXq1Ek+Pj7y9vbWNddco1GjRlnUOXPmjB555BFdc8018vDwUIMGDXTnnXfq4MGD5jpX9+Wio6PNyebHHnvM4r1qzy2x0dHRevTRRyVJQUFB5uuZOx1MQdfp5ZdfVrdu3dSgQQNVr15d7dq107PPPqvMzEyr/WzYsEE9e/Y0H3urVq20YMGCfOOSpK+//lr16tVT3759dfHixQLrAqWB4QaAE5s9e7Y6duyoBx98UPPnz1ePHj1Uq1Ytm3XXrVungQMHysfHx3xLYd5f02bMmKGbb75Zr7/+us6ePavo6GiFhYUpKSlJ11xzjaS//7G//vrrdd9998nX11cpKSlaunSpbrnlFu3fv1/16tWzaHPUqFHq06eP3n77bV28eFFubm76888/VbduXT399NOqX7++Tp8+rbfeekudOnVSUlKSrr/+equ4evToodjYWB05ckT//ve/NWTIELm6uqp9+/ZauXKlkpKSNGPGDNWsWVOLFi0ybztu3DjFxsZq0qRJeuaZZ3T69GnNmzdPoaGh2rt3r/z8/LRkyRLdf//9OnTokNatW2ex75ycHPXv31/btm3T1KlTFRoaqqNHj2rOnDkKCwvTzp07LX6x2717tw4cOKBZs2YpKChI1atXz/f6HTp0SCEhIRozZox8fHx05MgRLVy4UF26dNEPP/wgNze3fLc9fvy4unfvLjc3Ny1ZskR+fn569913bc6Z9s4772jYsGHq37+/3nrrLbm5uemVV15R7969tXHjRqtfNu+9917dd999Gj9+fIGdjaioKL399tt68sknddNNN+nixYvat2+fTp06Za6zceNG9evXT61atdLChQvVtGlTHTlyRJs2bcq33XXr1unll1/W8uXLtWHDBvn4+BRpFIE97/Xp06crJCREy5YtU7Vq1dSgQQNzx3zOnDny9/fXhQsXtG7dOoWFhenLL780/5GUlZWliIgIbdu2TZMnT9Ztt92mrKwsffPNN0pOTlZoaKhVTOnp6RoxYoQ+++wzffLJJ7rjjjvsPh4AQNVEP69i9/OkvxNp9913n6ZNm6Z58+bps88+05NPPqm//vpLixcvLnDbsWPHasWKFfr3v/+tXr16ad++fbr33nutfpxNTU1Vx44dVa1aNT3++ONq3ry5EhMT9eSTT+rIkSN68803LeovWrRILVq00PPPP69atWrpuuuus7n/xMREDR48WIMHD1Z0dLQ8PT119OhRffXVV+Y658+fV5cuXXTkyBE99thj6tSpky5cuKCtW7cqJSVFLVu2tGp3zJgxat++ve6991499NBDGjp0aJFG+I0ZM0anT5/WSy+9pLVr15pvn736Loz8rtOhQ4c0dOhQBQUFyd3dXXv37tVTTz2lgwcP6o033jBvv3z5co0dO1bdu3fXsmXL1KBBA/3888/mO0Bs+e9//6thw4Zp1KhReumll+Ti4mL3MQHFZgBwaps3bzYkGe+//75F+Zw5c4y8H+E2bdoY3bt3z7eNm2++2cjJyTGXHzlyxHBzczPGjBmT7/6zsrKMCxcuGNWrVzdefPFFc/mbb75pSDKGDRtW6DFkZWUZGRkZxnXXXWdMmTLFKq5+/fpZ1J88ebIhyZg0aZJF+d133234+vqanycmJhqSjBdeeMGi3u+//254eXkZU6dONZf16dPHCAwMtIpt5cqVhiRjzZo1FuXfffedIclYsmSJuSwwMNBwcXExfvrpp0KPOa+cnBwjMzPTOHr0qCHJ+Oijjwqs/9hjjxkmk8nYs2ePRXmvXr0MScbmzZsNwzCMixcvGr6+vlbnMDs722jfvr3RsWNHc1nue+bxxx+3K+a2bdsad999d4F1mjdvbjRv3ty4fPlyvnVy3yuHDx+2iuXkyZMWdbt3727zPZxXYe/1bt26FdpGVlaWkZmZafTs2dO45557zOUrVqwwJBmvvfZagdtLMh588EHj1KlTRpcuXYxGjRpZXS8AAApCP+9/Klo/r3v37jb7dGPHjjWqVatmHD16NN9tDxw4YEiyOF+GYRj
"text/plain": [
"<Figure size 1500x1000 with 4 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, ((ax0, ax1), (ax2, ax3)) = plt.subplots(nrows=2, ncols=2, figsize=(15,10))\n",
"\n",
"ax0.hist(scifi_fitpars_found[:,0], bins=100, density=True, alpha=0.5, histtype='bar', color=\"blue\", label=r\"$a_x$ found\")\n",
"ax0.hist(scifi_fitpars_lost[:,0], bins=100, density=True, alpha=0.5, histtype='bar', color=\"darkorange\", label=r\"$a_x$ lost\")\n",
"ax0.set_xlabel(\"a\")\n",
"ax0.set_ylabel(\"normed\")\n",
"ax0.set_title(\"fitparameter a der scifi track\")\n",
"ax0.legend()\n",
"\n",
"ax1.hist(scifi_fitpars_found[:,1], bins=100, density=True, alpha=0.5, histtype='bar', color=\"blue\", label=r\"$b_x$ found\")\n",
"ax1.hist(scifi_fitpars_lost[:,1], bins=100, density=True, alpha=0.5, histtype='bar', color=\"darkorange\", label=r\"$b_x$ lost\")\n",
"ax1.set_xlabel(\"b\")\n",
"ax1.set_ylabel(\"normed\")\n",
"ax1.set_title(\"fitparameter b der scifi track\")\n",
"ax1.legend()\n",
"\n",
"ax2.hist(scifi_fitpars_found[:,2], bins=200, density=True, alpha=0.5, histtype='bar', color=\"blue\", label=r\"$c_x$ found\")\n",
"ax2.hist(scifi_fitpars_lost[:,2], bins=1000, density=True, alpha=0.5, histtype='bar', color=\"darkorange\", label=r\"$c_x$ lost\")\n",
"ax2.set_xlim([-4e-5,4e-5])\n",
"ax2.set_xticks(np.arange(-4e-5,4.5e-5,1e-5),minor=False)\n",
"ax2.set_xlabel(\"c\")\n",
"ax2.set_ylabel(\"normed\")\n",
"ax2.set_title(\"fitparameter c der scifi track\")\n",
"ax2.legend()\n",
"\n",
"ax3.hist(scifi_fitpars_found[:,3], bins=200, density=True, alpha=0.5, histtype='bar', color=\"blue\", label=r\"$d_x$ found\")\n",
"ax3.hist(scifi_fitpars_lost[:,3], bins=1000, density=True, alpha=0.5, histtype='bar', color=\"darkorange\", label=r\"$d_x$ lost\")\n",
"ax3.set(xlim=(-5e-8,5e-8))\n",
"#ax3.text(-4e-8,3e8,\"d negligible <1e-7\")\n",
"ax3.set_xlabel(\"d\")\n",
"ax3.set_ylabel(\"normed\")\n",
"ax3.set_title(\"fitparameter d der scifi track\")\n",
"ax3.legend()\n",
"\n",
"\"\"\"\n",
"a_x: virtual hit on the reference plane\n",
"\"\"\"\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "code",
2023-10-10 16:17:10 +02:00
"execution_count": 17,
2023-09-25 16:36:13 +02:00
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"found\n",
"a = 18.04180550457298\n",
"b = 0.0057651944487250376\n",
"c = 9.480277789440454e-08\n",
"d = -4.452015498411874e-11\n",
"lost\n",
"a = -35.48371609662243\n",
"b = -0.010382326385451382\n",
"c = -6.208301288913938e-07\n",
"d = 9.580933791481267e-11\n"
]
}
],
"source": [
"print(\"found\")\n",
"print(\"a = \", str(np.mean(scifi_fitpars_found[:,0])))\n",
"print(\"b = \", str(np.mean(scifi_fitpars_found[:,1])))\n",
"print(\"c = \", str(np.mean(scifi_fitpars_found[:,2])))\n",
"print(\"d = \", str(np.mean(scifi_fitpars_found[:,3])))\n",
"\n",
"print(\"lost\")\n",
"print(\"a = \", str(np.mean(scifi_fitpars_lost[:,0])))\n",
"print(\"b = \", str(np.mean(scifi_fitpars_lost[:,1])))\n",
"print(\"c = \", str(np.mean(scifi_fitpars_lost[:,2])))\n",
"print(\"d = \", str(np.mean(scifi_fitpars_lost[:,3])))"
]
},
{
"cell_type": "code",
2023-10-10 16:17:10 +02:00
"execution_count": 18,
2023-09-25 16:36:13 +02:00
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-1.5438992626615335e-08"
]
},
2023-10-10 16:17:10 +02:00
"execution_count": 18,
2023-09-25 16:36:13 +02:00
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.min(scifi_fitpars_found[:,3])"
]
},
{
"cell_type": "code",
2023-10-10 16:17:10 +02:00
"execution_count": 19,
2023-09-25 16:36:13 +02:00
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABQEAAAIhCAYAAAD+cXzUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdeZhcVYH38e+9VXVr36v3vTvd6U5n6ewLhCSELSQsAirosLg76MygjqOjr7sDijOj46voOCqoiOACKDskgSSQrbMvnXR3Or2v1bXvdbf3Dx/yGhMgOiIE7+d58kedOnXuuXWS9K/PveceQdd1HYPBYDAYDAaDwWAwGAwGg8HwliW+0R0wGAwGg8FgMBgMBoPBYDAYDK8vYxLQYDAYDAaDwWAwGAwGg8FgeIszJgENBoPBYDAYDAaDwWAwGAyGtzhjEtBgMBgMBoPBYDAYDAaDwWB4izMmAQ0Gg8FgMBgMBoPBYDAYDIa3OGMS0GAwGAwGg8FgMBgMBoPBYHiLMyYBDQaDwWAwGAwGg8FgMBgMhrc4YxLQYDAYDAaDwWAwGAwGg8FgeIszJgENBoPBYDAYDAaDwWAwGAyGtzhjEtBg+Bvz0EMP0d7ejt1uRxAEDhw48EZ36ay++MUvIgjCa9Z74IEH+Na3vvX6d+hV1NfXs2HDhje0DwDRaJQbb7yR0tJSBEHg2muvfd2PWV9fz2233XZa2f79+1m1ahVerxdBEPjWt77FCy+8gCAIvPDCC6/a3n333YcgCOzZs+c1j7169WpWr1596nU2m+WLX/ziax7DYDAYDIa/RS//jB0YGHhd2r/nnnu47777zrn+nXfeyaOPPvq69OVcDAwMIAgC//7v//6G9eFlZ8tOr6eXz/2Px+tsvyecaya/7bbbcLlc53R8QRD44he/eOp1V1cXX/ziF1+3v5sGg+H/M7/RHTAYDH894XCYm2++mSuuuIJ77rkHq9VKS0vLG92t/5UHHniAI0eOcMcdd7zRXXnDfeUrX+GRRx7hxz/+MU1NTQQCgdf9mI888ggej+e0sve+971kMhkefPBB/H4/9fX1OBwOduzYwaxZs/5ix77nnntOe53NZvnSl74EcNrkoMFgMBgMhtffPffcQygUOuPi4Cu58847ueGGG/4qFy3f7M6WnV5PFRUV7Nixg6amplNlr/R7wvvf/36uuOKKv+jxd+zYQXV19anXXV1dfOlLX2L16tWv+7kbDH/rjElAg+FvSE9PD7Is83d/93esWrXqje7OX52qqiiKgtVqfaO78ro4cuQITU1NvPvd7/6rHXP+/Pln7ccHPvAB1q1bd1r5smXL/qLH/ktOKBoMBoPBYHjzyuVy2Gy2c7oj7Xz0Stnp9WK1Ws/IZa/0e4LD4Thtwu4v4S+dCQ0Gw7kzlgMbDH8jbrvtNi688EIA3vnOdyIIwml3S/3ud79j+fLlOBwO3G43l156KTt27DijjbNdnTvbMgFBEPjoRz/Kz372M9ra2nA4HMybN4/HH3/8jM8/8cQTdHR0YLVaaWhoOOdlGatXr+aJJ55gcHAQQRBO/YH/v8zh7rvv5qtf/SoNDQ1YrVaef/558vk8n/jEJ+jo6MDr9RIIBFi+fDm//e1vzziGpmn83//7f+no6MBut+Pz+Vi2bBm/+93vXrVv99xzD2azmS984Qunyr73ve8xb948XC4Xbreb1tZWPvOZz7zmeUajUW6//XaqqqqQJInGxkY++9nPUigUTjvXjRs3cuzYsVPfw6sti928eTOrV68mGAxit9upra3l+uuvJ5vNnqpTKBT48pe/TFtbGzabjWAwyJo1a9i+ffupOn+4HPjlZUaKovC9733vtPE41+XAL0ulUvz93/89oVCIYDDIddddx9jY2Gl1/nA58MDAACUlJQB86UtfOnXsl/sWDof54Ac/SE1NDVarlZKSEi644AI2btx4Tv0xGAwGg+Gt6sc//jHz5s3DZrMRCAR429vexrFjx06rc/LkSW688UYqKyuxWq2UlZWxdu3aU4+Vqa+v5+jRo2zZsuXUz+BXu6NLEAQymQw/+clPTtV/+Wf6y3ni2Wef5b3vfS8lJSU4HA4KhQInTpzgPe95D83NzTgcDqqqqrjqqqs4fPjwGceIx+N84hOfoLGxEavVSmlpKVdeeSXHjx9/xX7Jssytt96Ky+U6lVmz2Sz//M//TENDw6nvaNGiRfziF794ze/2yJEjXHPNNfj9fmw2Gx0dHfzkJz859f6rZadXci55cnR09FTukSSJyspKbrjhBiYnJ4EzlwO/2u8J57oc+GUnTpzgyiuvxOVyUVNTwyc+8YlTmfVlf7gc+L777uPtb387AGvWrDn1Hbzct/3797NhwwZKS0uxWq1UVlayfv16RkZGzrlPBoPh/zPuBDQY/kZ87nOfY8mSJXzkIx/hzjvvZM2aNaeWcT7wwAO8+93v5rLLLuMXv/gFhUKBu+++m9WrV7Np06ZToeBP9cQTT9DZ2cmXv/xlXC4Xd999N29729vo7u6msbERgE2bNnHNNdewfPlyHnzwQVRV5e677z4VUl7NPffcwwc/+EH6+vp45JFHzlrn29/+Ni0tLfz7v/87Ho+H5uZmCoUC0WiUf/7nf6aqqopiscjGjRu57rrruPfee7nllltOff62227j/vvv533vex9f/vKXkSSJffv2veIzS3Rd55Of/CTf/va3+eEPf3hqEurBBx/k9ttv5x/+4R/493//d0RR5MSJE3R1db3qOebzedasWUNfXx9f+tKXmDt3Ltu2beOuu+7iwIEDPPHEE6eWdNx+++0kEgl+/vOfA698p9zAwADr169n5cqV/PjHP8bn8zE6OsrTTz9NsVjE4XCgKArr1q1j27Zt3HHHHVx88cUoisLOnTsZGhpixYoVZ7S7fv16duzYwfLly7nhhhv4xCc+8arn9mre//73s379eh544AGGh4f55Cc/yd/93d+xefPms9avqKjg6aef5oorruB973sf73//+wFOTQzefPPN7Nu3j3/7t3+jpaWFeDzOvn37iEQif3YfDQaDwWA4391111185jOf4aabbuKuu+4iEonwxS9+keXLl9PZ2UlzczMAV1555amMVltby/T0NNu3bycejwO/fzzIDTfcgNfrPfW4jldbebFjxw4uvvhi1qxZw+c+9zmAsz5eZP369fzsZz8jk8lgsVgYGxsjGAzyta99jZKSEqLRKD/5yU9YunQp+/fvZ+bMmcDvLyZeeOGFDAwM8KlPfYqlS5eSTqfZunUr4+PjtLa2ntGneDzOddddx7Fjx9iyZQsLFy4E4OMf/zg/+9nP+OpXv8r8+fPJZDIcOXLkNTNEd3c3K1asoLS0lG9/+9sEg0Huv/9+brvtNiYnJ/mXf/mXPzk7nUueHB0dZfHixciyzGc+8xnmzp1LJBLhmWeeIRaLUVZWdka7r/Z7wp9ClmWuvvpq3ve+9/GJT3yCrVu38pWvfAWv18vnP//5s35m/fr13HnnnXzmM5/hu9/9LgsWLACgqamJTCbDpZdeSkNDA9/97ncpKytjYmKC559/nlQq9Sf3z2AwALrBYPib8fzzz+uA/qtf/epUmaqqemVlpT5nzhxdVdVT5alUSi8tLdVXrFhxquzWW2/V6+rqzmj3C1/4gv7H/50AellZmZ5MJk+VTUxM6KIo6nfdddepsqVLl+qVlZV6Lpc7VZZMJvVAIHBGm2ezfv36s/apv79fB/Smpia9WCy+ahuKouiyLOvve9/79Pnz558q37p1qw7on/3sZ1/183V1dfr69ev1bDarX3/99brX69U3btx4Wp2PfvSjus/ne83z+WPf//73dUD/5S9/eVr517/+dR3Qn3322VNlq1at0tvb21+zzV//+tc6oB84cOAV6/z0pz/VAf1//ud/XrWturo6/dZbbz2tDNA/8pGPnFb28t+9559//lXbu/fee3VAv/32208rv/vuu3VAHx8fP1W2atUqfdWqVadeh8NhHdC/8IUvnNGuy+XS77jjjlc9tsFgMBgMb2Uv/4zt7+/XdV3XY7GYbrfb9SuvvPK0ekNDQ7r
"text/plain": [
"<Figure size 1500x600 with 2 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, ((ax0, ax1)) = plt.subplots(nrows=1, ncols=2, figsize=(15,6))\n",
"\n",
"for i in range(0,ak.num(scifi_found[\"energy\"], axis=0)):\n",
" z_coord = np.linspace(scifi_z_found[i,0],12000,300)\n",
" fit = scifi_track(z_coord, *scifi_fitpars_found[i])\n",
" ax0.plot(z_coord, fit, \"-\", lw=0.5)\n",
" ax0.errorbar(ak.to_numpy(scifi_z_found[i,:]),ak.to_numpy(scifi_x_found[i,:]),fmt=\".\",ms=2)\n",
"\n",
"#ax0.legend()\n",
"ax0.set_xlabel(\"z [mm]\")\n",
"ax0.set_ylabel(\"x [mm]\")\n",
"ax0.set_title(\"found tracks of scifi hits\")\n",
"ax0.set(xlim=(7e3,12000), ylim=(-4000,4000))\n",
"ax0.grid()\n",
"\n",
"for i in range(0,ak.num(scifi_lost[\"energy\"], axis=0)):\n",
" z_coord = np.linspace(scifi_z_lost[i,0],12000,300)\n",
" fit = scifi_track(z_coord, *scifi_fitpars_lost[i])\n",
" ax1.plot(z_coord, fit, \"-\", lw=0.5)\n",
" ax1.errorbar(ak.to_numpy(scifi_z_lost[i,:]),ak.to_numpy(scifi_x_lost[i,:]),fmt=\".\",ms=2)\n",
"\n",
"#ax1.legend()\n",
"ax1.set_xlabel(\"z [mm]\")\n",
"ax1.set_ylabel(\"x [mm]\")\n",
"ax1.set_title(\"lost tracks of scifi hits\")\n",
"ax1.set(xlim=(7e3,12000), ylim=(-4000,4000))\n",
"ax1.grid()\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "code",
2023-10-10 16:17:10 +02:00
"execution_count": 20,
2023-09-25 16:36:13 +02:00
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"found \n",
"zmag = 5318.452650242005\n",
"lost \n",
"zmag = 5425.137423441005\n"
]
}
],
"source": [
"#vergleich der zmag werte\n",
"zmag_found = z_mag(xv_found, zv_found, tx_found, scifi_fitpars_found[:,0], scifi_fitpars_found[:,1])\n",
"zmag_lost = z_mag(xv_lost, zv_lost, tx_lost, scifi_fitpars_lost[:,0], scifi_fitpars_lost[:,1])\n",
"zmag_lost = zmag_lost[~np.isnan(zmag_lost)]\n",
"zmag_found = zmag_found[~np.isnan(zmag_found)]\n",
"\n",
"print(\"found \\nzmag = \", str(np.mean(zmag_found)))\n",
"print(\"lost \\nzmag = \", str(np.mean(zmag_lost)))"
]
},
{
"cell_type": "code",
2023-10-10 16:17:10 +02:00
"execution_count": 21,
2023-09-25 16:36:13 +02:00
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAloAAAHNCAYAAADG0NWMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABOZ0lEQVR4nO3deXwURf7/8feQkwQI4QqJQAigknCIJhgIpwrhUAFBiStGUMTNgiDgsoLKgniAigrItbgIHqvw8wsoKiJBAUEiYrgFlVUUxMRwmXAIuer3B9/Ml0kmdzqT4/V8POZhprq6urqY7vlYVVNtM8YYAQAAoMzVcHUFAAAAqioCLQAAAIsQaAEAAFiEQAsAAMAiBFoAAAAWIdACAACwCIEWAACARQi0AAAALEKgBQAAYBECLQAAAIsQaAEAAFiEQAsVyvbt2zV9+nT98ccfRco/ffp02Ww2nTx5ssB8y5cvl81m088//1ys+hS1/IrC2Xnm16YlbRMUn9VtXZryi3vNFUfO9VOeVq5cqTZt2qhmzZqy2Wzas2dPsa4L0DZljUALFcr27dv11FNPlfkFfuuttyohIUGBgYFlWm5F4+w882vT6tImKJhV15wrnDhxQrGxsWrZsqXWr1+vhIQEXXPNNcW6LkDblDV3V1cAKA8NGzZUw4YNXV0NyxXnPKtLm6D6+OGHH5SRkaF7771XPXr0sKf7+PhUuM/6hQsX5OPj4+pqlKvqeM4SPVqVSk43/L59+3TXXXfJz89P9erV08SJE5WZmanvv/9effv2Ve3atdW8eXO98MILDvv/97//1f3336+rr75aPj4+uuqqq3T77bdr//79eY71wQcfqH379vLy8lKLFi00d+5cp8MAOWnffvut/vKXv8jPz08BAQF64IEHlJqa6pD38OHDuueee9SoUSN5eXkpNDRUCxYscChr0qRJkqSQkBDZbDbZbDZt3ry5WO303XffqUWLFoqMjFRKSook50Mr3333nf7yl78oICBAXl5eatasme677z5dunSp2OU7k9M2u3fv1uDBg1WnTh35+fnp3nvv1YkTJxzybtu2Tbfccotq164tHx8fRUVF6eOPP3bIc+LECT300ENq2rSpvLy81LBhQ3Xp0kUbN26058l9ngW1aX7DTUWpS3H+3Z05dOiQvS65X35+fjLGFFqG1Qr7fBTneipJ+SNGjFDz5s3z7FfU4bii1K+wa66wa/ZKH3/8sTp06CAvLy+FhIRo9uzZRWqHb7/9VjabTe+99549LTExUTabTW3atHHIO2DAAIWHhzstZ8SIEerataskKSYmRjabTT179pRUvOuiONdtUT8DOWXu2rVLd955p/z9/dWyZcsSlVHS+3+O0t6Hi/qZyO+ci3Ifq2ro0aqEhg4dqnvvvVd//etfFR8frxdeeEEZGRnauHGjRo8erb///e9655139Nhjj6lVq1YaPHiwJOm3335T/fr1NWvWLDVs2FCnT5/WG2+8ocjISO3evVvXXnutJGn9+vUaPHiwunfvrpUrVyozM1OzZ8/W77//nm+dhgwZopiYGI0cOVL79+/XlClTJEmvv/66JOngwYOKiopSs2bN9NJLL6lx48b69NNPNW7cOJ08eVLTpk3Tgw8+qNOnT+vVV1/V6tWr7d38YWFhRW6bLVu26I477lD37t31zjvv5Pt/T3v37lXXrl3VoEEDzZgxQ1dffbWSkpK0du1apaeny8vLq1TlX+mOO+7Q0KFDFRcXp2+//VZTp07VwYMHtWPHDnl4eGjLli3q3bu32rdvr6VLl8rLy0sLFy7U7bffrnfffVcxMTGSpNjYWO3atUvPPvusrrnmGv3xxx/atWuXTp06le+xC2pTZ/N5ilqXHIX9u+cnKChICQkJDmkbN27U1KlT9cADD5T7vJ7civL5KOr1VNLyS6so9Svo81GUazbHZ599poEDB6pz585asWKFsrKy9MILLxR4z8jRpk0bBQYGauPGjbrrrrskXf4s1KxZUwcPHtRvv/2moKAgZWZmasuWLYqLi3NaztSpU3XjjTdqzJgxeu6553TTTTepTp06TvMWdN45AUVh121R2/hKgwcP1t133624uDidP3++RGWU9P4vlf4+XJzPRH7nXJL7WKVnUGlMmzbNSDIvvfSSQ3qHDh2MJLN69Wp7WkZGhmnYsKEZPHhwvuVlZmaa9PR0c/XVV5sJEybY0zt27GiaNm1qLl26ZE87e/asqV+/vsn9kcmp0wsvvOCQPnr0aOPt7W2ys7ONMcb06dPHNGnSxKSmpjrke/jhh423t7c5ffq0McaYF1980UgyR44cKUKL/N/xT5w4Yd566y3j6elpxo0bZ7KyshzyLVu2zKHcm2++2dStW9ekpKSUSfkF7Xtl2xpjzH/+8x8jybz99tvGGGM6depkGjVqZM6ePWvPk5mZadq2bWuaNGlib8NatWqZ8ePHF3jM3OdpTP5t6ixvUetS1H/3olqzZo3x8vIyf//734u1n1WK+vm4Un7XkzEl+/wNHz7cBAcH50nPafuCyi9O/fL7fBT1mjXGmMjISBMUFGT+/PNPe1paWpqpV69enro6c++995oWLVrY3/fq1cuMGjXK+Pv7mzfeeMMYY8yXX35pJJkNGzbkW86mTZuMJPPee+85pBfnuijqdetMfm2cU+Y///nPfPctahmluf+X9j5cnM9EfudclPtYVcPQYSV02223ObwPDQ2VzWZTv3797Gnu7u5q1aqVfvnlF3taZmamnnvuOYWFhcnT01Pu7u7y9PTU4cOHdejQIUnS+fPn9c0332jQoEHy9PS071urVi3dfvvt+dZpwIABDu/bt2+vixcvKiUlRRcvXtRnn32mO+64Qz4+PsrMzLS/+vfvr4sXL+qrr74qVZs8++yzGjFihGbNmqW5c+eqRo38P9oXLlzQli1bNHTo0CLP2yhO+bkNGzbM4f3QoUPl7u6uTZs26fz589qxY4fuvPNO1apVy57Hzc1NsbGx+vXXX/X9999Lkm688UYtX75czzzzjL766itlZGQUuQ5FUZy65Cjo372o3nrrLd1111168skn9eKLL5buJMpAUT8fRbmeSlN+aZW0fpKKdc2eP39eO3fu1ODBg+Xt7W0vo3bt2gXeM650yy236KefftKRI0d08eJFbdu2TX379tVNN92k+Ph4SZd7uby8vOzDg1Yr6LrNUdw2HjJkSJ604pZR0vt/ae/DJd0/9zlbfR+riAi0KqF69eo5vPf09JSPj4/DTS4n/eLFi/b3EydO1NSpUzVo0CB9+OGH2rFjh3bu3KnrrrtOf/75pyTpzJkzMsYoICAgz3GdpeWoX7++w/ucoY8///xTp06dUmZmpl599VV5eHg4vPr37y9JpV4+4e2339ZVV12lu+++u9C8Z86cUVZWlpo0aWJJ+bk1btzY4b27u7vq16+vU6dO2dvb2S//goKCJMnepb5y5UoNHz5c//73v9W5c2fVq1dP9913n5KTk4tdJ2eKU5ccBf27F8X8+fN1//3365VXXtGTTz7psG3u3Lm66667dM8996hOnTqKjIxUcnKyxo0bp3r16qlt27b2L5Jjx46pX79+atiwoerWrauHHnpI2dnZki5/kf3jH/9Q48aN1bJlSy1evNg+Pya/dijK56Mo11Npyi+tktZPUrGu2TNnzig7OzvP51zK+9nPT69evSRdDqa2bdumjIwM3XzzzerVq5c+++wz+7YuXbqoZs2aRW6D0ijous1R3DZ2dm0Vt4yS3v9Lex8u6f65z9nq+1hFxBytauTtt9/Wfffdp+eee84h/eTJk6pbt64kyd/fXzabzencipJeCP7+/vZekTFjxjjNExISUqKyc6xfv14xMTHq1q2bPvvsMwUHB+ebt169enJzc9Ovv/5qSfm5JScn66qrrrK/z8zM1KlTp1S/fn35+/urRo0aSkpKyrPfb7/9Jklq0KCB/b9z5szRnDlzdPToUa1du1aTJ09WSkqK1q9fX+T65Kc4dSkLzz77rKZPn67XX39d9913X57t+/bt09dff63Vq1fr9ddfV1RUlG655RbNnTt
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.hist(zmag_found, bins=100, density=True, alpha=0.5, histtype='bar', color=\"blue\", label=\"found\")\n",
"plt.hist(zmag_lost, bins=3000, density=True, alpha=0.5, histtype=\"bar\",color=\"darkorange\", label=\"lost\")\n",
"plt.xlabel(\"$z_{mag}$ [mm]\")\n",
"plt.ylabel(\"normed\")\n",
"plt.title(\"magnet kick position $z_{mag}$ calculated w fitparameters\")\n",
"plt.legend()\n",
"plt.xticks(np.arange(5100,5805,5), minor=True)\n",
"plt.yticks(np.arange(0,0.011,0.001), minor=True)\n",
"plt.xlim(5100,5800)\n",
"\n",
"\"\"\"\n",
"wir können definitiv einen unterschied zwischen den z_mag werten, bzw deren verteilungen für lost and found, erkennen\n",
"\"\"\"\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "code",
2023-10-10 16:17:10 +02:00
"execution_count": 22,
2023-09-25 16:36:13 +02:00
"metadata": {},
"outputs": [
{
2023-09-26 12:18:03 +02:00
"name": "stdout",
"output_type": "stream",
"text": [
"2591\n",
"1908\n",
2023-09-27 14:05:40 +02:00
"4499\n",
"5767\n"
2023-09-26 12:18:03 +02:00
]
2023-09-25 16:36:13 +02:00
}
],
"source": [
2023-09-26 13:35:31 +02:00
"both = ak.concatenate([found,lost],axis=0)\n",
2023-09-26 12:18:03 +02:00
"print(ak.num(found,axis=0))\n",
"print(ak.num(lost,axis=0))\n",
2023-09-27 14:05:40 +02:00
"print(ak.num(both,axis=0))\n",
"\n",
"tracks = allcolumns[(allcolumns.fromSignal) & (allcolumns.isElectron) & (~allcolumns.lost)]\n",
"alltracks = allcolumns[(allcolumns.fromSignal) & (allcolumns.isElectron)]\n",
"print(ak.num(tracks,axis=0))"
2023-09-25 16:36:13 +02:00
]
},
{
"cell_type": "code",
2023-10-10 16:17:10 +02:00
"execution_count": 23,
2023-09-25 16:36:13 +02:00
"metadata": {},
"outputs": [
{
2023-09-26 15:12:03 +02:00
"name": "stdout",
"output_type": "stream",
"text": [
2023-09-27 14:05:40 +02:00
"#events w/ shared track electrons from found and lost: 35\n",
"event_count: [359, 359]\n",
"velo idx: [25, 25]\n",
"mcp_index: [2926, 2936]\n",
2023-09-26 15:12:03 +02:00
"\n",
2023-09-27 14:05:40 +02:00
"velo x: [5.89, 5.89]\n",
"velo y: [9.81, 9.81]\n",
2023-09-26 15:12:03 +02:00
"\n",
2023-09-27 14:05:40 +02:00
"velo tx: [0.00824, 0.00824]\n",
"velo ty: [0.0135, 0.0135]\n",
"percentage of e with shared tracks: 0.7488\n"
2023-09-26 15:12:03 +02:00
]
2023-09-25 16:36:13 +02:00
}
],
"source": [
2023-09-26 15:12:03 +02:00
"#versuche teilchen von denselben events mit shared tracks zu finden. \n",
"#idee: alle teilchen eines events sind durch event_count auffindbar. \n",
2023-09-27 14:05:40 +02:00
"a_f_itr = tracks[\"event_count\"].to_numpy()\n",
2023-09-26 15:12:03 +02:00
"f_itr = np.unique(a_f_itr)\n",
"\n",
"shared = ak.ArrayBuilder()\n",
2023-09-26 12:18:03 +02:00
"count = 0\n",
"\n",
"for itr in f_itr:\n",
2023-09-27 14:05:40 +02:00
" temp = alltracks[alltracks[\"event_count\"]==itr]\n",
2023-09-26 12:18:03 +02:00
" if ak.num(temp,axis=0)>1:\n",
" #iterate over cols in temp and append all with duplicate velo_track_idx, such that possibles is array with possible shared tracks particles\n",
2023-09-26 13:35:31 +02:00
" #idea: look at one event at a time and at one velo track at a time. if there are at least two e with the same velo_track_idx in the same event,\n",
" #concatenate to the array of other shared track particles\n",
2023-09-26 12:18:03 +02:00
" _jitr = temp[\"velo_track_idx\"].to_numpy()\n",
" jitr = np.unique(_jitr)\n",
" for jentry in jitr:\n",
" jtem = temp[temp[\"velo_track_idx\"]==jentry]\n",
" if ak.num(jtem,axis=0)>1:\n",
2023-09-26 15:12:03 +02:00
" shared.append(jtem)\n",
2023-09-26 12:18:03 +02:00
" else:\n",
" continue\n",
" else:\n",
" continue\n",
2023-09-26 15:12:03 +02:00
"shared = ak.Array(shared)\n",
"\n",
"idx=0\n",
"print(\"#events w/ shared track electrons from found and lost: \",ak.num(shared, axis=0))\n",
"\n",
"print(\"event_count: \", shared[idx,:,\"event_count\"])\n",
"print(\"velo idx: \" ,shared[idx,:,\"velo_track_idx\"])\n",
"print(\"mcp_index: \", shared[idx,:,\"mcp_idx\"])\n",
"\n",
"print(\"\\nvelo x: \" ,shared[idx,:,\"velo_track_x\"])\n",
"print(\"velo y: \" ,shared[idx,:,\"velo_track_y\"])\n",
"\n",
"print(\"\\nvelo tx: \" ,shared[idx,:,\"velo_track_tx\"])\n",
2023-09-27 13:00:51 +02:00
"print(\"velo ty: \" ,shared[idx,:,\"velo_track_ty\"])\n",
"\n",
"\n",
"\n",
2023-09-28 12:12:21 +02:00
"print(\"percentage of e with shared tracks: \", np.round((ak.num(shared,axis=0)*2)/(ak.num(alltracks,axis=0))*100,4)) #error for percentage"
2023-09-26 12:18:03 +02:00
]
},
{
"cell_type": "code",
2023-10-10 16:17:10 +02:00
"execution_count": 24,
2023-09-26 12:18:03 +02:00
"metadata": {},
2023-09-26 13:35:31 +02:00
"outputs": [
{
2023-09-26 15:12:03 +02:00
"name": "stdout",
"output_type": "stream",
"text": [
"#velo_track_idx in all events: 217\n",
"velo idx: [0, 0, 0]\n",
"mcp_index: [1066, 1251, 666]\n",
"event_count: [5735, 7049, 7378]\n",
"\n",
"velo x: [-2.21, 9.45, -18.5]\n",
"velo y: [-21.6, -33.4, 17.4]\n",
"\n",
"velo tx: [-0.0022, 0.0153, -0.0224]\n",
"velo ty: [-0.0263, -0.0469, 0.0212]\n"
]
2023-09-26 13:35:31 +02:00
}
],
2023-09-26 12:18:03 +02:00
"source": [
2023-09-26 15:12:03 +02:00
"#electrons with same velo_track_idx from all events\n",
2023-09-26 13:35:31 +02:00
"temp_ = found[\"velo_track_idx\"].to_numpy()\n",
"temp = np.unique(temp_)\n",
"count=0\n",
"psb=ak.ArrayBuilder()\n",
"for jentry in temp:\n",
" jtem = found[found[\"velo_track_idx\"]==jentry]\n",
" if ak.num(jtem,axis=0)>1:\n",
" psb.append(jtem)\n",
" else:\n",
" continue\n",
"\n",
"psb = ak.Array(psb)\n",
"\n",
2023-09-26 15:12:03 +02:00
"print(\"#velo_track_idx in all events: \",ak.num(psb,axis=0))\n",
"idx = 0\n",
"print(\"velo idx: \" ,psb[idx,:,\"velo_track_idx\"])\n",
"print(\"mcp_index: \", psb[idx,:,\"mcp_idx\"])\n",
"print(\"event_count: \", psb[idx,:,\"event_count\"])\n",
"print(\"\\nvelo x: \" ,psb[idx,:,\"velo_track_x\"])\n",
"print(\"velo y: \" ,psb[idx,:,\"velo_track_y\"])\n",
"\n",
"print(\"\\nvelo tx: \" ,psb[idx,:,\"velo_track_tx\"])\n",
2023-09-27 13:00:51 +02:00
"print(\"velo ty: \" ,psb[idx,:,\"velo_track_ty\"])\n",
"\n"
2023-09-26 13:35:31 +02:00
]
},
{
"cell_type": "code",
2023-09-26 15:12:03 +02:00
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
2023-10-10 16:17:10 +02:00
"execution_count": 25,
2023-09-26 13:35:31 +02:00
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
2023-09-26 15:12:03 +02:00
"#events w/ shared track electrons from Photon Conversions: 951\n",
"shared_idx: 0\n",
"event_count: [11, 11]\n",
"velo idx: [27, 27]\n",
"mcp_index: [1211, 1215]\n",
2023-09-26 13:35:31 +02:00
"\n",
2023-09-26 15:12:03 +02:00
"velo x: [-20.4, -20.4]\n",
"velo y: [-24.4, -24.4]\n",
2023-09-26 13:35:31 +02:00
"\n",
2023-09-26 15:12:03 +02:00
"velo tx: [-0.0234, -0.0234]\n",
2023-09-27 13:00:51 +02:00
"velo ty: [-0.028, -0.028]\n",
"percentage of e with shared tracks: 5.59\n"
2023-09-26 13:35:31 +02:00
]
}
],
"source": [
2023-09-26 15:12:03 +02:00
"#generell: wie viele elektronen von photon conversions teilen sich einen track?\n",
"conv = allcolumns[(allcolumns.fromPairProd) & (allcolumns.isElectron) & (~allcolumns.lost)]\n",
2023-09-26 13:35:31 +02:00
"\n",
2023-09-26 15:12:03 +02:00
"conv_it = conv[\"event_count\"].to_numpy()\n",
"conv_itr = np.unique(conv_it)\n",
"\n",
"cshared = ak.ArrayBuilder()\n",
"count = 0\n",
"\n",
"for itr in conv_itr:\n",
" temp = conv[conv[\"event_count\"]==itr]\n",
" if ak.num(temp,axis=0)>1:\n",
" #iterate over cols in temp and append all with duplicate velo_track_idx, such that possibles is array with possible shared tracks particles\n",
" #idea: look at one event at a time and at one velo track at a time. if there are at least two e with the same velo_track_idx in the same event,\n",
" #concatenate to the array of other shared track particles\n",
" _jitr = temp[\"velo_track_idx\"].to_numpy()\n",
" jitr = np.unique(_jitr)\n",
" for jentry in jitr:\n",
" jtem = temp[temp[\"velo_track_idx\"]==jentry]\n",
" if ak.num(jtem,axis=0)>1:\n",
" cshared.append(jtem)\n",
" else:\n",
" continue\n",
" else:\n",
" continue\n",
"cshared = ak.Array(cshared)\n",
"\n",
"\n",
"print(\"#events w/ shared track electrons from Photon Conversions: \",ak.num(cshared, axis=0))\n",
"idx = 0\n",
"print(\"shared_idx: \", idx)\n",
"print(\"event_count: \", cshared[idx,:,\"event_count\"])\n",
"print(\"velo idx: \" ,cshared[idx,:,\"velo_track_idx\"])\n",
"print(\"mcp_index: \", cshared[idx,:,\"mcp_idx\"])\n",
"\n",
"print(\"\\nvelo x: \" ,cshared[idx,:,\"velo_track_x\"])\n",
"print(\"velo y: \" ,cshared[idx,:,\"velo_track_y\"])\n",
"\n",
"print(\"\\nvelo tx: \" ,cshared[idx,:,\"velo_track_tx\"])\n",
2023-09-27 13:00:51 +02:00
"print(\"velo ty: \" ,cshared[idx,:,\"velo_track_ty\"])\n",
"\n",
2023-09-28 12:12:21 +02:00
"print(\"percentage of e with shared tracks: \", np.round((ak.num(cshared,axis=0)*2)/(ak.num(conv,axis=0))*100,4))\n",
"\n",
"#constrain z pos of brem vtx; only interested in conversions in velo z<770\n"
2023-09-26 12:18:03 +02:00
]
},
{
"cell_type": "code",
2023-10-10 16:17:10 +02:00
"execution_count": 26,
2023-09-26 12:18:03 +02:00
"metadata": {},
2023-09-27 13:00:51 +02:00
"outputs": [
{
"data": {
"text/plain": [
"34025"
]
},
2023-10-10 16:17:10 +02:00
"execution_count": 26,
2023-09-27 13:00:51 +02:00
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ak.num(conv,axis=0)"
]
2023-09-25 16:36:13 +02:00
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
2023-09-26 12:18:03 +02:00
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
2023-09-25 16:36:13 +02:00
"metadata": {},
2023-09-26 15:12:03 +02:00
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
2023-09-25 16:36:13 +02:00
},
{
"cell_type": "code",
2023-10-10 16:17:10 +02:00
"execution_count": 27,
2023-09-25 16:36:13 +02:00
"metadata": {},
2023-10-10 16:17:10 +02:00
"outputs": [
{
"data": {
"text/plain": [
"'\\nwas meinst du? velo track sharen? wie überprüfe ich das? versuche ich einfach e teilchen mit identischen velo tracks und slopes zu finden?\\nbenutze ich mother key um elektronen von gleichen events zu finden? meinst du wie viele elektronen vom selben event teilen sich einen track,\\noder wie viele teilchen von unterschiedlichen events haben dieselben tracks?\\n'"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
2023-09-25 16:36:13 +02:00
"source": [
"\"\"\"\n",
"was meinst du? velo track sharen? wie überprüfe ich das? versuche ich einfach e teilchen mit identischen velo tracks und slopes zu finden?\n",
"benutze ich mother key um elektronen von gleichen events zu finden? meinst du wie viele elektronen vom selben event teilen sich einen track,\n",
"oder wie viele teilchen von unterschiedlichen events haben dieselben tracks?\n",
"\"\"\"\n"
]
},
{
"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": "env1",
"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.11.5"
},
"orig_nbformat": 4
},
"nbformat": 4,
"nbformat_minor": 2
}
