You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

893 lines
391 KiB

1 year ago
1 year ago
1 year ago
1 year ago
1 year ago
1 year ago
1 year ago
1 year ago
1 year ago
1 year ago
1 year ago
1 year ago
1 year ago
1 year ago
1 year ago
1 year ago
1 year ago
1 year ago
1 year ago
1 year ago
1 year ago
1 year ago
1 year ago
1 year ago
1 year ago
  1. {
  2. "cells": [
  3. {
  4. "cell_type": "code",
  5. "execution_count": 1,
  6. "metadata": {},
  7. "outputs": [],
  8. "source": [
  9. "import uproot\t\n",
  10. "import numpy as np\n",
  11. "import sys\n",
  12. "import os\n",
  13. "import matplotlib\n",
  14. "import matplotlib.pyplot as plt\n",
  15. "from mpl_toolkits import mplot3d\n",
  16. "import itertools\n",
  17. "import awkward as ak\n",
  18. "from scipy.optimize import curve_fit\n",
  19. "from mpl_toolkits.axes_grid1 import ImageGrid\n",
  20. "%matplotlib inline"
  21. ]
  22. },
  23. {
  24. "cell_type": "code",
  25. "execution_count": 2,
  26. "metadata": {},
  27. "outputs": [
  28. {
  29. "data": {
  30. "text/plain": [
  31. "10522"
  32. ]
  33. },
  34. "execution_count": 2,
  35. "metadata": {},
  36. "output_type": "execute_result"
  37. }
  38. ],
  39. "source": [
  40. "file = uproot.open(\"tracking_losses_ntuple_Bd2KstEE.root:PrDebugTrackingLosses.PrDebugTrackingTool/Tuple;1\")\n",
  41. "\n",
  42. "#selektiere nur elektronen von B->K*ee und nur solche mit einem momentum von ueber 5 GeV \n",
  43. "allcolumns = file.arrays()\n",
  44. "found = allcolumns[(allcolumns.isElectron) & (~allcolumns.lost) & (allcolumns.fromSignal) & (allcolumns.p > 5e3)] #B: 9056\n",
  45. "lost = allcolumns[(allcolumns.isElectron) & (allcolumns.lost) & (allcolumns.fromSignal) & (allcolumns.p > 5e3)] #B: 1466\n",
  46. "\n",
  47. "ak.num(found, axis=0) + ak.num(lost, axis=0)\n",
  48. "#ak.count(found, axis=None)"
  49. ]
  50. },
  51. {
  52. "cell_type": "code",
  53. "execution_count": 3,
  54. "metadata": {},
  55. "outputs": [
  56. {
  57. "name": "stdout",
  58. "output_type": "stream",
  59. "text": [
  60. "eff all = 0.8606728758791105 +/- 0.003375885792719708\n"
  61. ]
  62. }
  63. ],
  64. "source": [
  65. "def t_eff(found, lost, axis = 0):\n",
  66. " sel = ak.num(found, axis=axis)\n",
  67. " des = ak.num(lost, axis=axis)\n",
  68. " return sel/(sel + des)\n",
  69. "\n",
  70. "def eff_err(found, lost):\n",
  71. " n_f = ak.num(found, axis=0)\n",
  72. " n_all = ak.num(found, axis=0) + ak.num(lost,axis=0)\n",
  73. " return 1/n_all * np.sqrt(np.abs(n_f*(1-n_f/n_all)))\n",
  74. "\n",
  75. "\n",
  76. "print(\"eff all = \", t_eff(found, lost), \"+/-\", eff_err(found, lost))"
  77. ]
  78. },
  79. {
  80. "cell_type": "code",
  81. "execution_count": 4,
  82. "metadata": {},
  83. "outputs": [],
  84. "source": [
  85. "#try excluding all photons that originate from a vtx @ z>9500mm\n",
  86. "#ignore all brem vertices @ z>9500mm \n",
  87. "\n",
  88. "#found\n",
  89. "\n",
  90. "brem_e_f = found[\"brem_photons_pe\"]\n",
  91. "brem_z_f = found[\"brem_vtx_z\"]\n",
  92. "e_f = found[\"energy\"]\n",
  93. "length_f = found[\"brem_vtx_z_length\"]\n",
  94. "\n",
  95. "brem_f = ak.ArrayBuilder()\n",
  96. "\n",
  97. "for itr in range(ak.num(found,axis=0)):\n",
  98. " brem_f.begin_record()\n",
  99. " #[:,\"energy\"] energy\n",
  100. " brem_f.field(\"energy\").append(e_f[itr])\n",
  101. " #[:,\"photon_length\"] number of vertices\n",
  102. " brem_f.field(\"photon_length\").integer(length_f[itr])\n",
  103. " #[:,\"brem_photons_pe\",:] photon energy \n",
  104. " brem_f.field(\"brem_photons_pe\").append(brem_e_f[itr])\n",
  105. " #[:,\"brem_vtx_z\",:] brem vtx z\n",
  106. " brem_f.field(\"brem_vtx_z\").append(brem_z_f[itr])\n",
  107. " brem_f.end_record()\n",
  108. "\n",
  109. "brem_f = ak.Array(brem_f)\n",
  110. "\n",
  111. "#lost\n",
  112. "\n",
  113. "brem_e_l = lost[\"brem_photons_pe\"]\n",
  114. "brem_z_l = lost[\"brem_vtx_z\"]\n",
  115. "e_l = lost[\"energy\"]\n",
  116. "length_l = lost[\"brem_vtx_z_length\"]\n",
  117. "\n",
  118. "brem_l = ak.ArrayBuilder()\n",
  119. "\n",
  120. "for itr in range(ak.num(lost,axis=0)):\n",
  121. " brem_l.begin_record()\n",
  122. " #[:,\"energy\"] energy\n",
  123. " brem_l.field(\"energy\").append(e_l[itr])\n",
  124. " #[:,\"photon_length\"] number of vertices\n",
  125. " brem_l.field(\"photon_length\").integer(length_l[itr])\n",
  126. " #[:,\"brem_photons_pe\",:] photon energy \n",
  127. " brem_l.field(\"brem_photons_pe\").append(brem_e_l[itr])\n",
  128. " #[:,\"brem_vtx_z\",:] brem vtx z\n",
  129. " brem_l.field(\"brem_vtx_z\").append(brem_z_l[itr])\n",
  130. " brem_l.end_record()\n",
  131. "\n",
  132. "brem_l = ak.Array(brem_l)"
  133. ]
  134. },
  135. {
  136. "cell_type": "code",
  137. "execution_count": 5,
  138. "metadata": {},
  139. "outputs": [],
  140. "source": [
  141. "cut_brem_found = ak.ArrayBuilder()\n",
  142. "\n",
  143. "for itr in range(ak.num(brem_f, axis=0)):\n",
  144. " cut_brem_found.begin_record()\n",
  145. " cut_brem_found.field(\"energy\").real(brem_f[itr,\"energy\"])\n",
  146. " \n",
  147. " cut_brem_found.field(\"brem_photons_pe\")\n",
  148. " cut_brem_found.begin_list()\n",
  149. " for jentry in range(brem_f[itr, \"photon_length\"]):\n",
  150. " if brem_f[itr, \"brem_vtx_z\", jentry]>9500:\n",
  151. " continue\n",
  152. " else:\n",
  153. " cut_brem_found.real(brem_f[itr,\"brem_photons_pe\", jentry])\n",
  154. " \n",
  155. " #cut_brem_found.field(\"brem_vtx_z\").real(brem_f[itr, \"brem_vtx_z\",jentry])\n",
  156. " cut_brem_found.end_list()\n",
  157. " \n",
  158. " cut_brem_found.field(\"brem_vtx_z\")\n",
  159. " cut_brem_found.begin_list()\n",
  160. " for jentry in range(brem_f[itr, \"photon_length\"]):\n",
  161. " if brem_f[itr, \"brem_vtx_z\", jentry]>9500:\n",
  162. " continue\n",
  163. " else:\n",
  164. " cut_brem_found.real(brem_f[itr, \"brem_vtx_z\",jentry])\n",
  165. " cut_brem_found.end_list()\n",
  166. " \n",
  167. "\n",
  168. " \n",
  169. " cut_brem_found.end_record()\n",
  170. "\n",
  171. "cut_brem_found = ak.Array(cut_brem_found)\n",
  172. "\n",
  173. "\n",
  174. "\n",
  175. "cut_brem_lost = ak.ArrayBuilder()\n",
  176. "\n",
  177. "for itr in range(ak.num(brem_l, axis=0)):\n",
  178. " cut_brem_lost.begin_record()\n",
  179. " cut_brem_lost.field(\"energy\").real(brem_l[itr,\"energy\"])\n",
  180. " \n",
  181. " \n",
  182. " cut_brem_lost.field(\"brem_photons_pe\")\n",
  183. " cut_brem_lost.begin_list()\n",
  184. " for jentry in range(brem_l[itr, \"photon_length\"]):\n",
  185. " if brem_l[itr, \"brem_vtx_z\", jentry]>9500:\n",
  186. " continue\n",
  187. " else:\n",
  188. " cut_brem_lost.real(brem_l[itr,\"brem_photons_pe\", jentry])\n",
  189. " \n",
  190. " #cut_brem_found.field(\"brem_vtx_z\").real(brem_f[itr, \"brem_vtx_z\",jentry])\n",
  191. " cut_brem_lost.end_list()\n",
  192. " \n",
  193. " cut_brem_lost.field(\"brem_vtx_z\")\n",
  194. " cut_brem_lost.begin_list()\n",
  195. " for jentry in range(brem_l[itr, \"photon_length\"]):\n",
  196. " if brem_l[itr, \"brem_vtx_z\", jentry]>9500:\n",
  197. " continue\n",
  198. " else:\n",
  199. " cut_brem_lost.real(brem_l[itr, \"brem_vtx_z\",jentry])\n",
  200. " cut_brem_lost.end_list()\n",
  201. " \n",
  202. " cut_brem_lost.end_record()\n",
  203. "\n",
  204. "cut_brem_lost = ak.Array(cut_brem_lost)\n"
  205. ]
  206. },
  207. {
  208. "cell_type": "code",
  209. "execution_count": 6,
  210. "metadata": {},
  211. "outputs": [
  212. {
  213. "data": {
  214. "text/html": [
  215. "<pre>{energy: 9.36e+03,\n",
  216. " brem_photons_pe: [2.47e+03, 170, 224, 388, 3.23e+03, 809, 172, 224],\n",
  217. " brem_vtx_z: [400, 501, 638, 667, 677, 709, 8.58e+03, 9.28e+03]}\n",
  218. "---------------------------------------------------------------------\n",
  219. "type: {\n",
  220. " energy: float64,\n",
  221. " brem_photons_pe: var * float64,\n",
  222. " brem_vtx_z: var * float64\n",
  223. "}</pre>"
  224. ],
  225. "text/plain": [
  226. "<Record {energy: 9.36e+03, ...} type='{energy: float64, brem_photons_pe: va...'>"
  227. ]
  228. },
  229. "execution_count": 6,
  230. "metadata": {},
  231. "output_type": "execute_result"
  232. }
  233. ],
  234. "source": [
  235. "#data in cut_brem_found and cut_brem_lost\n",
  236. "\n",
  237. "cut_length_found = ak.num(cut_brem_found[\"brem_photons_pe\"],axis=-1)\n",
  238. "cut_length_lost = ak.num(cut_brem_lost[\"brem_photons_pe\"], axis=-1)\n",
  239. "\n",
  240. "cut_brem_found[1]\n"
  241. ]
  242. },
  243. {
  244. "cell_type": "code",
  245. "execution_count": 7,
  246. "metadata": {},
  247. "outputs": [
  248. {
  249. "data": {
  250. "text/plain": [
  251. "8"
  252. ]
  253. },
  254. "execution_count": 7,
  255. "metadata": {},
  256. "output_type": "execute_result"
  257. }
  258. ],
  259. "source": [
  260. "cut_length_found[1]"
  261. ]
  262. },
  263. {
  264. "cell_type": "markdown",
  265. "metadata": {},
  266. "source": [
  267. "### Split in Upstream and Downstream Events and analyse separately"
  268. ]
  269. },
  270. {
  271. "cell_type": "code",
  272. "execution_count": 8,
  273. "metadata": {},
  274. "outputs": [],
  275. "source": [
  276. "#try to find a split between energy lost before and after the magnet (z~5000mm)\n",
  277. "\n",
  278. "upstream_found = ak.ArrayBuilder()\n",
  279. "downstream_found = ak.ArrayBuilder()\n",
  280. "\n",
  281. "for itr in range(ak.num(cut_brem_found, axis=0)):\n",
  282. " upstream_found.begin_record()\n",
  283. " upstream_found.field(\"energy\").real(cut_brem_found[itr,\"energy\"])\n",
  284. " \n",
  285. " downstream_found.begin_record()\n",
  286. " downstream_found.field(\"energy\").real(cut_brem_found[itr,\"energy\"])\n",
  287. " \n",
  288. " upstream_found.field(\"brem_photons_pe\")\n",
  289. " downstream_found.field(\"brem_photons_pe\")\n",
  290. " upstream_found.begin_list()\n",
  291. " downstream_found.begin_list()\n",
  292. " for jentry in range(cut_length_found[itr]):\n",
  293. " if (cut_brem_found[itr, \"brem_vtx_z\", jentry]>5000):\n",
  294. " if cut_brem_found[itr, \"brem_vtx_z\", jentry]<=9500:\n",
  295. " downstream_found.real(cut_brem_found[itr,\"brem_photons_pe\",jentry])\n",
  296. " else:\n",
  297. " continue\n",
  298. " else:\n",
  299. " upstream_found.real(cut_brem_found[itr,\"brem_photons_pe\", jentry]) \n",
  300. " upstream_found.end_list()\n",
  301. " downstream_found.end_list()\n",
  302. " \n",
  303. " upstream_found.field(\"brem_vtx_z\")\n",
  304. " downstream_found.field(\"brem_vtx_z\")\n",
  305. " upstream_found.begin_list()\n",
  306. " downstream_found.begin_list()\n",
  307. " for jentry in range(cut_length_found[itr]):\n",
  308. " if cut_brem_found[itr, \"brem_vtx_z\", jentry]>5000:\n",
  309. " if cut_brem_found[itr,\"brem_vtx_z\",jentry]<=9500:\n",
  310. " downstream_found.real(cut_brem_found[itr,\"brem_vtx_z\",jentry])\n",
  311. " else:\n",
  312. " continue\n",
  313. " else:\n",
  314. " upstream_found.real(cut_brem_found[itr, \"brem_vtx_z\",jentry])\n",
  315. " upstream_found.end_list()\n",
  316. " downstream_found.end_list()\n",
  317. " upstream_found.end_record()\n",
  318. " downstream_found.end_record()\n",
  319. " \n",
  320. "\n",
  321. "upstream_found = ak.Array(upstream_found)\n",
  322. "downstream_found = ak.Array(downstream_found)\n",
  323. "\n",
  324. "\n",
  325. "upstream_lost = ak.ArrayBuilder()\n",
  326. "downstream_lost = ak.ArrayBuilder()\n",
  327. "\n",
  328. "for itr in range(ak.num(cut_brem_lost, axis=0)):\n",
  329. " upstream_lost.begin_record()\n",
  330. " upstream_lost.field(\"energy\").real(cut_brem_lost[itr,\"energy\"])\n",
  331. " \n",
  332. " downstream_lost.begin_record()\n",
  333. " downstream_lost.field(\"energy\").real(cut_brem_lost[itr,\"energy\"])\n",
  334. " \n",
  335. " upstream_lost.field(\"brem_photons_pe\")\n",
  336. " downstream_lost.field(\"brem_photons_pe\")\n",
  337. " upstream_lost.begin_list()\n",
  338. " downstream_lost.begin_list()\n",
  339. " for jentry in range(cut_length_lost[itr]):\n",
  340. " if (cut_brem_lost[itr, \"brem_vtx_z\", jentry]>5000):\n",
  341. " if cut_brem_lost[itr, \"brem_vtx_z\", jentry]<=9500:\n",
  342. " downstream_lost.real(cut_brem_lost[itr,\"brem_photons_pe\",jentry])\n",
  343. " else:\n",
  344. " continue\n",
  345. " else:\n",
  346. " upstream_lost.real(cut_brem_lost[itr,\"brem_photons_pe\", jentry]) \n",
  347. " upstream_lost.end_list()\n",
  348. " downstream_lost.end_list()\n",
  349. " \n",
  350. " upstream_lost.field(\"brem_vtx_z\")\n",
  351. " downstream_lost.field(\"brem_vtx_z\")\n",
  352. " upstream_lost.begin_list()\n",
  353. " downstream_lost.begin_list()\n",
  354. " for jentry in range(cut_length_lost[itr]):\n",
  355. " if cut_brem_lost[itr, \"brem_vtx_z\", jentry]>5000:\n",
  356. " if cut_brem_lost[itr,\"brem_vtx_z\",jentry]<=9500:\n",
  357. " downstream_lost.real(cut_brem_lost[itr,\"brem_vtx_z\",jentry])\n",
  358. " else:\n",
  359. " continue\n",
  360. " else:\n",
  361. " upstream_lost.real(cut_brem_lost[itr, \"brem_vtx_z\",jentry])\n",
  362. " upstream_lost.end_list()\n",
  363. " downstream_lost.end_list()\n",
  364. " upstream_lost.end_record()\n",
  365. " downstream_lost.end_record()\n",
  366. " \n",
  367. "\n",
  368. "upstream_lost = ak.Array(upstream_lost)\n",
  369. "downstream_lost = ak.Array(downstream_lost)\n"
  370. ]
  371. },
  372. {
  373. "cell_type": "code",
  374. "execution_count": 9,
  375. "metadata": {},
  376. "outputs": [
  377. {
  378. "data": {
  379. "text/html": [
  380. "<pre>{energy: 4.62e+04,\n",
  381. " brem_photons_pe: [3.26e+03, 4.45e+03, 178, 1.45e+04, 1.1e+03, 3.79e+03],\n",
  382. " brem_vtx_z: [162, 187, 387, 487, 1.34e+03, 2.32e+03]}\n",
  383. "-------------------------------------------------------------------------\n",
  384. "type: {\n",
  385. " energy: float64,\n",
  386. " brem_photons_pe: var * float64,\n",
  387. " brem_vtx_z: var * float64\n",
  388. "}</pre>"
  389. ],
  390. "text/plain": [
  391. "<Record {energy: 4.62e+04, ...} type='{energy: float64, brem_photons_pe: va...'>"
  392. ]
  393. },
  394. "execution_count": 9,
  395. "metadata": {},
  396. "output_type": "execute_result"
  397. }
  398. ],
  399. "source": [
  400. "upstream_found[0]"
  401. ]
  402. },
  403. {
  404. "cell_type": "code",
  405. "execution_count": null,
  406. "metadata": {},
  407. "outputs": [],
  408. "source": []
  409. },
  410. {
  411. "cell_type": "code",
  412. "execution_count": 13,
  413. "metadata": {},
  414. "outputs": [
  415. {
  416. "name": "stdout",
  417. "output_type": "stream",
  418. "text": [
  419. "\n",
  420. "upstream: cutoff energy = 350MeV, sample size: 1562\n",
  421. "eff = 0.9181 +/- 0.007\n"
  422. ]
  423. }
  424. ],
  425. "source": [
  426. "#plot efficiency against cutoff energy \n",
  427. "up_efficiencies = []\n",
  428. "up_deff = []\n",
  429. "\n",
  430. "\n",
  431. "for cutoff_energy in range(0,10050,200):\n",
  432. "\tup_nobrem_f = upstream_found[ak.sum(upstream_found[\"brem_photons_pe\"],axis=-1,keepdims=False)<cutoff_energy]\n",
  433. "\tup_nobrem_l = upstream_lost[ak.sum(upstream_lost[\"brem_photons_pe\"],axis=-1,keepdims=False)<cutoff_energy]\n",
  434. "\t\n",
  435. "\tif ak.num(up_nobrem_f,axis=0)+ak.num(up_nobrem_l,axis=0)==0:\n",
  436. "\t\tup_efficiencies.append(0)\n",
  437. "\t\tup_deff.append(0)\n",
  438. "\t\tcontinue\n",
  439. "\n",
  440. "\teff = t_eff(up_nobrem_f,up_nobrem_l)\n",
  441. "\tdeff = eff_err(up_nobrem_f,up_nobrem_l)\n",
  442. "\tup_efficiencies.append(eff)\n",
  443. "\tup_deff.append(deff)\n",
  444. "\n",
  445. "\t#print(\"\\ncutoff = \",str(cutoff_energy),\"MeV, sample size: \",ak.num(up_nobrem_f,axis=0)+ak.num(up_nobrem_l,axis=0))\n",
  446. "\t#print(\"eff = \",np.round(eff,4), \"+/-\", np.round(eff_err(up_nobrem_f, up_nobrem_l),4))\n",
  447. "\n",
  448. "\"\"\"\n",
  449. "we see that a cutoff energy of xxxMeV is ideal because the efficiency drops significantly for higher values\n",
  450. "\"\"\"\n",
  451. "cutoff_energy = 350.0 #MeV\n",
  452. "\n",
  453. "\"\"\"\n",
  454. "better statistics: cutoff=xxxMeV - sample size: xxx events and efficiency=xxxx\n",
  455. "\"\"\"\n",
  456. "up_nobrem_found = upstream_found[ak.sum(upstream_found[\"brem_photons_pe\"],axis=-1,keepdims=False)<cutoff_energy]\n",
  457. "up_nobrem_lost = upstream_lost[ak.sum(upstream_lost[\"brem_photons_pe\"],axis=-1,keepdims=False)<cutoff_energy]\n",
  458. "\n",
  459. "print(\"\\nupstream: cutoff energy = 350MeV, sample size:\",ak.num(up_nobrem_found,axis=0)+ak.num(up_nobrem_lost,axis=0))\n",
  460. "print(\"eff = \",np.round(t_eff(up_nobrem_found, up_nobrem_lost),4), \"+/-\", np.round(eff_err(up_nobrem_found, up_nobrem_lost),3))\n"
  461. ]
  462. },
  463. {
  464. "cell_type": "code",
  465. "execution_count": 14,
  466. "metadata": {},
  467. "outputs": [
  468. {
  469. "name": "stdout",
  470. "output_type": "stream",
  471. "text": [
  472. "nobrem_vertices\n",
  473. "upstream: cutoff energy = 350MeV, sample size: 1562\n",
  474. "eff = 0.9181 +/- 0.007\n",
  475. "\n",
  476. "downstream: cutoff energy = 350MeV, sample size: 5131\n",
  477. "eff = 0.8864 +/- 0.004\n"
  478. ]
  479. }
  480. ],
  481. "source": [
  482. "down_efficiencies = []\n",
  483. "down_deff = []\n",
  484. "\n",
  485. "for cutoff_energy in range(0,10050,200):\n",
  486. "\tdown_nobrem_f = downstream_found[ak.sum(downstream_found[\"brem_photons_pe\"],axis=-1,keepdims=False)<cutoff_energy]\n",
  487. "\tdown_nobrem_l = downstream_lost[ak.sum(downstream_lost[\"brem_photons_pe\"],axis=-1,keepdims=False)<cutoff_energy]\n",
  488. "\n",
  489. "\tif ak.num(down_nobrem_f,axis=0)+ak.num(down_nobrem_l,axis=0)==0:\n",
  490. "\t\tdown_efficiencies.append(0)\n",
  491. "\t\tdown_deff.append(0)\n",
  492. "\t\tcontinue\n",
  493. "\teff = t_eff(down_nobrem_f,down_nobrem_l)\n",
  494. "\tdeff = eff_err(down_nobrem_f,down_nobrem_l)\n",
  495. "\tdown_efficiencies.append(eff)\n",
  496. "\tdown_deff.append(deff)\n",
  497. "\n",
  498. "\n",
  499. "\t#print(\"\\ncutoff = \",str(cutoff_energy),\"MeV, sample size: \",ak.num(down_nobrem_f,axis=0)+ak.num(down_nobrem_l,axis=0))\n",
  500. "\t#print(\"eff = \",np.round(eff,4), \"+/-\", np.round(eff_err(down_nobrem_f, down_nobrem_l),4))\n",
  501. "\n",
  502. "\"\"\"\n",
  503. "we see that a cutoff energy of xxxMeV is ideal because the efficiency drops significantly for higher values\n",
  504. "\"\"\"\n",
  505. "cutoff_energy = 350.0 #MeV\n",
  506. "\n",
  507. "\"\"\"\n",
  508. "better statistics: cutoff=xxxMeV - sample size: xxx events and efficiency=xxxx\n",
  509. "\"\"\"\n",
  510. "down_nobrem_found = downstream_found[ak.sum(downstream_found[\"brem_photons_pe\"],axis=-1,keepdims=False)<cutoff_energy]\n",
  511. "down_nobrem_lost = downstream_lost[ak.sum(downstream_lost[\"brem_photons_pe\"],axis=-1,keepdims=False)<cutoff_energy]\n",
  512. "\n",
  513. "\n",
  514. "print(\"nobrem_vertices\\nupstream: cutoff energy = 350MeV, sample size:\",ak.num(up_nobrem_found,axis=0)+ak.num(up_nobrem_lost,axis=0))\n",
  515. "print(\"eff = \",np.round(t_eff(up_nobrem_found, up_nobrem_lost),4), \"+/-\", np.round(eff_err(up_nobrem_found, up_nobrem_lost),3))\n",
  516. "\n",
  517. "print(\"\\ndownstream: cutoff energy = 350MeV, sample size:\",ak.num(down_nobrem_found,axis=0)+ak.num(down_nobrem_lost,axis=0))\n",
  518. "print(\"eff = \",np.round(t_eff(down_nobrem_found, down_nobrem_lost),4), \"+/-\", np.round(eff_err(down_nobrem_found, down_nobrem_lost),3))\n"
  519. ]
  520. },
  521. {
  522. "cell_type": "code",
  523. "execution_count": 17,
  524. "metadata": {},
  525. "outputs": [
  526. {
  527. "data": {
  528. "image/png": "iVBORw0KGgoAAAANSUhEUgAABcMAAAJJCAYAAABmnC9/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAACW70lEQVR4nOzdeViU9f7/8RcCsii4hCIqCrZipiUuiZppiplLZiZWR7OwkwePa3WUzFxaPFqZbVpaZpqlJ7NsoZIsOxoabu2drG96KANJKyAJGHB+f/hjDhPb3MA9M9w8H9fllXPzmfv+3PMeppuXH963j91utwsAAAAAAAAAAAtr5OkJAAAAAAAAAABgNsJwAAAAAAAAAIDlEYYDAAAAAAAAACyPMBwAAAAAAAAAYHmE4QAAAAAAAAAAyyMMBwAAAAAAAABYHmE4AAAAAAAAAMDyCMMBAAAAAAAAAJZHGA4AAAAAAAAAsDzCcAAAAAAAAACA5RGGAwAAAAAAAAAsjzAcAADUa5MmTdK6des8PQ246LPPPlNGRoanpwEAAACgASIMBwAAgNusXr1a77//vqenAQAAAKABIgwHAAAN2ubNm3XhhRcqKChIPj4++uSTTzw9pSqlpKTIx8fH6U9oaKhiY2P1r3/9y23HX79+vdP2X3/9VcOGDVPjxo31xBNPOH3t3//+t2699VZlZ2c7tn311VcaM2aMjh8/bvqc65OdO3eWq2/pn71795Yb/9lnnykxMVFnn322goKCFBQUpHPPPVe33Xab9u/fX6M5XHPNNQoKCtJvv/1W6Zgbb7xR/v7+Hq/fwoUL5ePjoxMnTnh0Hu6ybt06+fj46OjRo6bsPy0tTQsXLqyy9gAAAPUZYTgAAKh3RowYoebNm6t58+Z68cUXlZSU5Hj8z3/+0+X9/Pzzz5owYYLOPvtsvfPOO9qzZ4/OO+88E2deewcPHpQkbdu2TXv27FFaWprWrFmjU6dO6frrr9dnn33mluPHxsY6tn322Wfq0aOHDh06pPfff19///vfnZ4TGxur8PBwde3aVTt27NATTzyhIUOG6Morr1SrVq1MnW999cADD2jPnj1Of7p06eI05umnn1ZsbKw+/vhjzZgxQ2+++abeeustzZw5U19++aV69uyp//u//zN87MTERBUUFOjFF1+s8Os5OTl69dVXNWLECIWHh9fo/OCd0tLStGjRIsJwAABgWX6engAAAIBRb775puPvkyZN0uWXX65JkyYZ3s/hw4dls9n0l7/8RQMGDKh0XH5+voKDg2sy1Tp38OBBNWvWTKNGjXJs69Onj4qLi/WXv/xFhw4dUteuXU09flBQkC644AJJ0qZNm5SYmKiuXbvqlVdeUdu2bcs9p0mTJlq8eLFOnz6tJUuWqFGjRtq4caPGjx9v2jw94ZdfftHp06cVFhZW632de+65uvTSSyv9+kcffaSkpCQNHz5cW7ZsUePGjR1fGzRokKZOnaqXX35ZQUFBho89bNgwtW3bVmvXrlVSUlK5r7/00kv6448/lJiYaHjf3sabvrfrG147AABQH7EyHAAAeIXdu3crPj5ezZo1U4sWLTR8+HB9++23ph1v0qRJ6tevnyQpISFBPj4+uvzyyx1tFw4ePKixY8eqRYsWOvvss53mecUVVygkJETBwcGKi4vTW2+95bTv0n189tlnuu6669SsWTO1bNlSs2fPVnFxsb755htdeeWVCgkJUVRUlJYtW+byvA8cOKCLL7643PYff/xRkhQTE1ODV0N65JFH9Nprrxk6/h133KHrr79eN954oz788MMKg3BJ+vzzz9W9e3cdPHhQV111lW655Rbdf//9Gj58eLXtLVx9X9T1+2fkyJHq0aOH1qxZo27duikoKEiRkZFasGCBTp8+XeFzPvvsM0VERGjYsGFav3698vLyanz86jzwwAPy9fXV008/7RSEl3XdddeVq8m3336rG264Qa1bt1ZAQIBiYmL05JNPOo3x9fXVTTfdpAMHDujzzz8vt9/nnnvOcZ5GlX5vfPnll7r++uvVrFkzhYeH65ZbblFOTo7TWFe+10r98MMPGjNmjEJDQ9WsWTP95S9/0c8//1zhsSv73nbltTHze9uV49f2uf/5z390/fXXKzw8XAEBAerQoYMmTpyowsJCLVy4UHfeeackKTo62tGeZ+fOnXXyuVj29auu/j///LP++te/KjIyUgEBAWrVqpX69u2r9957z+XXEwAAoCKE4QAAwOMWLlyoAQMGKDIyUi+99JKeeeYZ/fDDD7riiiv0+++/V/ncdevW1WhV+Pz58x1hUWlLipUrVzq+PmbMGJ1zzjl6+eWX9dRTT0mSPvzwQw0aNEg5OTl69tln9dJLLykkJEQjR47U5s2byx1j3Lhx6tatm1555RXdeuuteuSRRzRr1iyNHj1aw4cP16uvvqpBgwZpzpw52rp1a7VzPnnypDIyMtStWzcVFxeruLhY2dnZ2rBhg+6//35NnjxZvXr1MvxaSNL+/fs1bty4KgPx0uN36NBB8fHxeuKJJ7RmzRqtXr260kBWklq0aKH77rtP77zzjqKjo9W3b1998sknGjdunJo1a1bp81x9X9Tm/VOZAwcO6D//+Y8eeeQR3XnnnXr99dfVr18/LV68WGvXrq3wOZdeeqmef/55+fv769Zbb1Xr1q113XXXaevWrSosLDR0/KlTp8rPz0+hoaEaOnSodu/e7fhaSUmJPvjgA/Xo0UMREREu7/Orr75Sz5499cUXX+jhhx/Wm2++qeHDh2v69OlatGiR09hbbrlFPj4+5c71q6++Unp6um666Sb5+voaOqeyrr32Wp133nl65ZVXNHfuXL344ouaNWuW4+tGv9euueYanXPOOdqyZYsWLlyo1157TUOHDpXNZis3tqLvbSOvjVT339tGj1+T53766afq2bOn9u7dq8WLF+vtt9/WkiVLVFhYqKKiIk2ePFnTpk2TJG3dutXRnqd79+5VvnZGayVVX/8JEybotdde0z333KPt27frmWee0eDBg3Xy5MlqX0sAAIAq2QEAADzojTfesEuyL1u2zGn74cOH7ZLsL7zwQrnnXHnllfYmTZpU+Of+++93+dgffPCBXZL95ZdfdmxbsGCBXZL9nnvuKTf+0ksvtbdu3dqel5fn2FZcXGzv0qWLvX379vbTp0877ePhhx92ev7FF19sl2TfunWrY5vNZrO3atXKPmbMmGrnu337drukcn/8/Pzs9913n8vnXZHi4mL7DTfcYPf397e/+uqr1R4/MDDQvnfvXsPHmTp1qv25556rdpyr74uavH+q8+OPP9ol2Tt16mT/7bffHNuLiorsbdq0sY8YMaLaffz666/2tWvX2uPj4+1+fn72Zs2a2SdNmmR/99137cXFxZU+7+DBg/YZM2bYX331Vfu///1v+9q1a+0xMTF2X19f+zvvvGO32+32rKwsuyT7+PHjyz2/uLjYbrPZHH9K35N2u90+dOhQe/v27e05OTlOz/n73/9uDwwMtP/yyy9O2wcMGGAPCwuzFxUVObbdfvvtdkn2w4cPV/saVKT0e+PP9UpKSrIHBgY65mv0e23WrFlO+9u4cWO5+lf1ve3qa2PW97arx3/uuefskuxHjhwx/NxBgwbZmzdvbs/Ozq50Hg8++GC5/Zc979p8LpbdT3X1b9q0qX3mzJmVzhMAAKCmWBkOAAA86p577tHZZ5+tGTNmOFY7FxcXKzo6WkFBQfr+++/LPeftt9/W77//XuGfu+66q07mde211zo9PnXqlD7++GONHTtWTZs2dWz39fXVhAkT9OOPP+qbb75xes6IESOcHsfExMjHx8epvYSfn5/OOecc/fe//612TgcOHJB0ZtXmvn37tG/fPr3zzjsaPny47rnnnkpXoJ44ccLR8qCyP35+fnrxxRdls9k0btw4HT9+vNLjT5gwQQUFBfr3v/9d7Zz/7IknnnBpJb+r74uavH+qs2/fPklnVpyXXbnu7++vc845p9rWLpLUvHlz3XzzzXr33XeVmZmpf/7znzpy5IiuvPJKtW3bttJ5XXLJJVqxYoV
  529. "text/plain": [
  530. "<Figure size 1800x600 with 2 Axes>"
  531. ]
  532. },
  533. "metadata": {},
  534. "output_type": "display_data"
  535. }
  536. ],
  537. "source": [
  538. "#plot efficiencies wrt cutoff energy\n",
  539. "fig, ax = plt.subplots(nrows=1,ncols=2,figsize=(18,6))\n",
  540. "x_ = np.arange(0,10050,step=200)\n",
  541. "\n",
  542. "ax[0].errorbar(x_,up_efficiencies,yerr=up_deff, ls=\"\", capsize=1,fmt=\".\")\n",
  543. "ax[0].set(xlabel=\"cutoff energy [MeV]\",ylabel=r\"$\\epsilon$\",title=\"upstream\", ylim=[0.8,1.0])\n",
  544. "ax[0].set_yticks(np.arange(0.8,1.01,step=0.02),minor=False)\n",
  545. "ax[0].set_xticks(np.arange(0,10100,step=200),minor=True)\n",
  546. "ax[0].grid()\n",
  547. "\n",
  548. "ax[1].errorbar(x_,down_efficiencies,yerr=down_deff, ls=\"\", capsize=1,fmt=\".\")\n",
  549. "ax[1].set(xlabel=\"cutoff energy [MeV]\",ylabel=r\"$\\epsilon$\",title=\"downstream\", ylim=[0.8,1.0])\n",
  550. "ax[1].set_yticks(np.arange(0.8,1.01,step=0.02),minor=False)\n",
  551. "ax[1].set_xticks(np.arange(0,10100,step=200),minor=True)\n",
  552. "ax[1].grid(True)\n",
  553. "\n",
  554. "fig.suptitle(r\"$e^\\pm$ from $B\\rightarrow K^\\ast ee$, $p>5$GeV, nobrem electrons\")\n",
  555. "\n",
  556. "\n",
  557. "plt.show()"
  558. ]
  559. },
  560. {
  561. "cell_type": "code",
  562. "execution_count": null,
  563. "metadata": {},
  564. "outputs": [],
  565. "source": []
  566. },
  567. {
  568. "cell_type": "code",
  569. "execution_count": 13,
  570. "metadata": {},
  571. "outputs": [
  572. {
  573. "name": "stdout",
  574. "output_type": "stream",
  575. "text": [
  576. "brem vertices\n",
  577. "upstream eff = 0.851 +/- 0.004\n",
  578. "downstream eff = 0.836 +/- 0.005\n"
  579. ]
  580. }
  581. ],
  582. "source": [
  583. "cutoff_energy=350\n",
  584. "#possibly: instead of checking if any photons exceed the cutoff, use the sum of all photon energies to separate nobrem and brem\n",
  585. "\n",
  586. "upstream_brem_found = upstream_found[ak.sum(upstream_found[\"brem_photons_pe\"],axis=-1,keepdims=False)>=cutoff_energy]\n",
  587. "up_energy_found = ak.to_numpy(upstream_brem_found[\"energy\"])\n",
  588. "up_eph_found = ak.to_numpy(ak.sum(upstream_brem_found[\"brem_photons_pe\"], axis=-1, keepdims=False))\n",
  589. "up_residual_found = up_energy_found - up_eph_found\n",
  590. "up_energyloss_found = up_eph_found/up_energy_found\n",
  591. "\n",
  592. "\n",
  593. "upstream_brem_lost = upstream_lost[ak.sum(upstream_lost[\"brem_photons_pe\"],axis=-1,keepdims=False)>=cutoff_energy]\n",
  594. "up_energy_lost = ak.to_numpy(upstream_brem_lost[\"energy\"])\n",
  595. "up_eph_lost = ak.to_numpy(ak.sum(upstream_brem_lost[\"brem_photons_pe\"], axis=-1, keepdims=False))\n",
  596. "up_residual_lost = up_energy_lost - up_eph_lost\n",
  597. "up_energyloss_lost = up_eph_lost/up_energy_lost\n",
  598. "\n",
  599. "\n",
  600. "print(\"brem vertices\\nupstream eff = \", np.round(t_eff(upstream_brem_found,upstream_brem_lost),3), \"+/-\", np.round(eff_err(upstream_brem_found, upstream_brem_lost),3))\n",
  601. "\n",
  602. "\n",
  603. "downstream_brem_found = downstream_found[ak.sum(downstream_found[\"brem_photons_pe\"],axis=-1,keepdims=False)>=cutoff_energy]\n",
  604. "down_energy_found = ak.to_numpy(downstream_brem_found[\"energy\"])\n",
  605. "down_eph_found = ak.to_numpy(ak.sum(downstream_brem_found[\"brem_photons_pe\"], axis=-1, keepdims=False))\n",
  606. "down_residual_found = down_energy_found - down_eph_found\n",
  607. "down_energyloss_found = down_eph_found/down_energy_found\n",
  608. "\n",
  609. "\n",
  610. "downstream_brem_lost = downstream_lost[ak.sum(downstream_lost[\"brem_photons_pe\"],axis=-1,keepdims=False)>=cutoff_energy]\n",
  611. "down_energy_lost = ak.to_numpy(downstream_brem_lost[\"energy\"])\n",
  612. "down_eph_lost = ak.to_numpy(ak.sum(downstream_brem_lost[\"brem_photons_pe\"], axis=-1, keepdims=False))\n",
  613. "down_residual_lost = down_energy_lost - down_eph_lost\n",
  614. "down_energyloss_lost = down_eph_lost/down_energy_lost\n",
  615. "\n",
  616. "\n",
  617. "print(\"downstream eff = \", np.round(t_eff(downstream_brem_found,downstream_brem_lost),3), \"+/-\", np.round(eff_err(downstream_brem_found, downstream_brem_lost),3))"
  618. ]
  619. },
  620. {
  621. "cell_type": "code",
  622. "execution_count": 14,
  623. "metadata": {},
  624. "outputs": [
  625. {
  626. "name": "stdout",
  627. "output_type": "stream",
  628. "text": [
  629. "upstream:\n",
  630. "mean energyloss relative to initial energy (found): 0.33078325542598164\n",
  631. "mean energyloss relative to initial energy (lost): 0.5708618852236069\n",
  632. "downstream:\n",
  633. "mean energyloss relative to initial energy (found): 0.19104090843883118\n",
  634. "mean energyloss relative to initial energy (lost): 0.3051594568487781\n"
  635. ]
  636. }
  637. ],
  638. "source": [
  639. "print(\"upstream:\\nmean energyloss relative to initial energy (found): \",ak.mean(up_energyloss_found))\n",
  640. "print(\"mean energyloss relative to initial energy (lost): \", ak.mean(up_energyloss_lost))\n",
  641. "\n",
  642. "print(\"downstream:\\nmean energyloss relative to initial energy (found): \",ak.mean(down_energyloss_found))\n",
  643. "print(\"mean energyloss relative to initial energy (lost): \", ak.mean(down_energyloss_lost))"
  644. ]
  645. },
  646. {
  647. "cell_type": "code",
  648. "execution_count": 15,
  649. "metadata": {},
  650. "outputs": [
  651. {
  652. "data": {
  653. "image/png": "iVBORw0KGgoAAAANSUhEUgAABbYAAAJNCAYAAADtbSO1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAB1tklEQVR4nO3deZyN9f//8eeZfcYy1jH2nezKlqWyK6GVSklI9UHSIrQw+ijbT2khJdEnkcqSykdUhpCspeITRZmyb2MQZnn//ug7J8fMMNc423XN4367za3O+1znfb3OeR0zr3nN+7wvlzHGCAAAAAAAAAAAmwgJdAAAAAAAAAAAAFhBYxsAAAAAAAAAYCs0tgEAAAAAAAAAtkJjGwAAAAAAAABgKzS2AQAAAAAAAAC2QmMbAAAAAAAAAGArNLYBAAAAAAAAALZCYxsAAAAAAAAAYCs0tgEAAAAAAAAAtkJjGwAABI2tW7dqz549gQ4DAAAAABDkaGwDAICg8eabb+qrr74KdBgAAAAAgCBHYxsAAATUqlWr1L9/fx08eNA9tm3bNt166606cOBAACMDAAAAAAQrGtsAADjIkiVL5HK5PL4KFy6sRo0a6YMPPvDLuf/zn/94jB87dkw33HCDIiIi9Nprr2V5XKNGjVSqVCnVr19fX375pV577TV16NBB119/vUqWLOnTmO0kMTExS24zv9atW5ftY7Zu3ap+/fqpatWqio6OVnR0tKpXr64HH3xQGzdutBzDLbfcoujoaB0/fjzHY+6++26Fh4f79I8SCQkJcrlcOnz4sFfmW7t2rRISEi76vOwuPT1dcXFxeumll3I8xtuvq104Of/r169Xp06dVKhQIRUsWFBt2rTRmjVrshxn9fvLyZMnNWTIEJUpU0ZRUVFq2LCh3n///TwfBwAAkBdhgQ4AAAB4z+bNmyVJH3/8seLi4mSM0Z49ezRq1CjddddduuKKK1S/fn2fnrtRo0busa1bt+qWW27RqVOn9NVXX6lVq1ZZHlegQAE999xzysjI0NixYxUSEqL33ntPd955p0/itLsXXnhBbdq08RirW7duluPeeOMNDRo0SDVr1tQjjzyiOnXqyOVyafv27Zo7d66aNGmiX375RVWrVs31ufv166dFixZpzpw5GjBgQJb7k5OTtXDhQnXp0kWlSpWy/uQCZO3atRo9erTuu+8+FSlSJNDh+MSqVat06NAh3XrrrYEOJeg4Nf8bNmzQtddeq6ZNm+rdd9+VMUYTJkxQu3bttGLFCjVv3jzLY3L7/eXWW2/Vhg0bNG7cONWoUUNz5szRXXfdpYyMDPXs2dPycQAAAHlBYxsAAAfZvHmzYmNj1a1bN/dY8+bNlZaWpnvuuUdbtmzxaWM7OjpaV1xxhSTp/fffV79+/VS/fn3Nnz9fZcqUyfZxP/zwg3r16qX4+Hh17txZZcqU0fPPP693331X77zzjkqUKOGTeP3p6NGjysjI8MpzqV69uq6++uqLHrNmzRoNGDBAN954oz766CNFRES472vbtq0GDhyoDz/8UNHR0ZbOfcMNN6hMmTJ6++23s21sz507V3/99Zf69etnaV743kcffaTGjRurYsWKPjvH6dOnFRMT47P5kdW+fftUoEABFS5cOMt9zz77rIoUKaKlS5e689K+fXtVqVJFTzzxRLYrt3Pz/WXJkiVavny5u0ktSW3atNHvv/+uoUOH6o477lBoaGiujwMAAMgrtiIBAMBBNm3apIYNG2YZ/+OPPyRJtWrVsjznSy+9pEWLFlk69xNPPKG77rpLd999t1auXJljU1uSihYtqjFjxmjp0qWqXLmyWrZsqe+++049evRQbGxsjo9bvXq1OnbsqNjYWBUtWlQ33nijdu7cmefjcqtr165q3Lixpk+frgYNGig6Olrly5fXqFGjlJGRke1jtm7dqtKlS+uGG27Qf/7zH6WkpOT5/LnxwgsvKDQ0VG+88YZHU/t83bt3z5KXnTt3qmfPnoqLi1NkZKRq1aqlKVOmuO8PDQ1V7969tWnTJv3www9Z5pw5c6b7eVqVuQ3Gli1bdOutt6pw4cKKjY3VPffco0OHDmX7mAMHDuiuu+5SbGysSpUqpb59+yo5OdnjmNWrV6tdu3YqVKiQYmJi1KJFC3322Wce5x06dKgkqXLlyu7tFxITE3M9x/nx//TTT5eM6dChQ3rggQdUvnx5RUZGqmTJkmrZsqW++OKLHF+fn376SS6XSx9++KF7bNOmTXK5XKpTp47Hsd26dfP45IQxRgsXLtRtt92W4/znS0pKumQOMp/v5s2bdfvtt6to0aIeq/8v9V46f46tW7eqe/fuio2NVbFixfTYY48pLS1NP//8s66//noVKlRIlSpV0oQJE3IVf6ZFixbJ5XLpyy+/zHLf66+/LpfLpVtvvTXH/J85c0ZXXnmlqlWr5pHD/fv3Kz4+Xq1bt1Z6enquYslpmw+Xy6XffvvN0vM6duyYZsyYofbt26tcuXLatWtXtsetWbNGrVu39vhjQ6FChXTttddq7dq12rdvn6XzZlq4cKEKFiyo7t27e4z36dNHe/fu1bfffmvpOOny3wu+fi8BAIDgRGMbAACHOHLkiPbs2aMGDRooLS1NaWlpOnjwoN599109//zzuv/++9W0aVPL827cuFE9evS4aHM789wVKlRQx44d9dprr2n69Ol68803c2ysZipXrpy6dOniMZbZQA0PD8/2MQkJCbruuutUvnx5zZ07V2+99ZaSkpLUrl07nTx50vJxVmzatEn/+9//9NJLL2no0KFavHixWrVqpeeee05vv/12to+5+uqr9c477yg8PFz9+/dXXFycunfvrgULFujs2bOWzj9w4ECFhYWpcOHC6tSpk1avXu1xf3p6ulasWKHGjRurdOnSuZ5327ZtatKkiX788UdNmjRJn376qW688UYNHjxYo0ePdh/Xt29fuVyuLM9127ZtWr9+vXr37n1ZqzBvueUWVatWTR999JESEhK0aNEiderUSampqVmOve2221SjRg3Nnz9fw4cP15w5c/Too4+671+5cqXatm2r5ORkzZgxQ3PnzlWhQoXUtWtXzZs3T5J0//336+GHH5YkLViwQN98842++eYbXXXVVbmew0pMktSrVy8tWrRII0eO1LJly/TWW2+pffv2OnLkSI6vS506dVS6dGmP5vcXX3yh6Ohobdu2TXv37pUkpaWlaeXKlWrfvr37uMwmZm4b21ZycOutt6patWr68MMPNW3aNEm5fy9l6tGjhxo0aKD58+erf//+eumll/Too4/q5ptv1o033qiFCxeqbdu2GjZsmBYsWJCr5yBJXbp0UVxcnGbOnJnlvlmzZumqq67SK6+8kmP+o6Ki9MEHH+jgwYPq27evJCkjI0N33323jDGaO3durt/rmfNmfn311VcqW7as4uPjVaxYsUs+/vTp05o3b55uuukmxcfH6+GHH1aRIkU0b9481a5dO9vHnDt3TpGRkVnGM8ey++PUpb6/SNKPP/6oWrVqKSzM88O/mZ8G+vHHHy0dd77LfS/46r0EAACClAEAAI6wbNkyIynLV1hYmBkzZkye501LSzM9e/Y04eHhZuHChZc8d1RUlFm3bl2ez3cpn3zyiZFkJkyY4DG+Y8cOI8nMnj3b0nFW/PHHH0aSqVKlijl+/Lh7/Ny5cyY+Pt506dLlknMcO3bMvP3226Zjx44mLCzMxMbGmvvuu898/vnnJi0tLcfHbd682TzyyCNm4cKFZtWqVebtt982tWrVMqGhoWbp0qXu4/bv328kmTvvvDPLHGlpaSY1NdX9lZGR4b6vU6dOply5ciY5OdnjMYMGDTJRUVHm6NGj7rHrrrvOlChRwpw7d8499vjjjxtJZseOHZd8DbIzatQoI8k8+uijHuPvvfdelnxlHnthbgcMGGCioqLcz+vqq682cXFxJiUlxeM1qFu3rilXrpz7uIkTJxpJZvfu3Vniyu0cuY3JGGMKFixohgwZYuXlMcYYc88995gqVaq4b7dv397079/fFC1a1LzzzjvGGGPWrFljJJlly5a5jxs
  654. "text/plain": [
  655. "<Figure size 1800x600 with 2 Axes>"
  656. ]
  657. },
  658. "metadata": {},
  659. "output_type": "display_data"
  660. }
  661. ],
  662. "source": [
  663. "#in abhängigkeit von der energie der elektronen\n",
  664. "fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(18,6))\n",
  665. "\n",
  666. "\n",
  667. "ax[0].hist(up_energyloss_lost, bins=100, density=True, alpha=0.5, histtype='bar', color=\"darkorange\", label=\"lost\")\n",
  668. "ax[0].hist(up_energyloss_found, bins=100, density=True, alpha=0.5, histtype='bar', color=\"blue\", label=\"found\")\n",
  669. "ax[0].set_xticks(np.arange(0,1.1,0.1), minor=True,)\n",
  670. "ax[0].set_yticks(np.arange(0,11,1), minor=True)\n",
  671. "ax[0].set_xlabel(r\"$E_\\gamma/E_0$\")\n",
  672. "ax[0].set_ylabel(\"counts (normed)\")\n",
  673. "ax[0].set_title(\"Upstream\")\n",
  674. "ax[0].legend()\n",
  675. "ax[0].grid()\n",
  676. "\n",
  677. "ax[1].hist(down_energyloss_lost, bins=100, density=True, alpha=0.5, histtype='bar', color=\"darkorange\", label=\"lost\")\n",
  678. "ax[1].hist(down_energyloss_found, bins=100, density=True, alpha=0.5, histtype='bar', color=\"blue\", label=\"found\")\n",
  679. "ax[1].set_xticks(np.arange(0,1.1,0.1), minor=True,)\n",
  680. "ax[1].set_yticks(np.arange(0,11,1), minor=True)\n",
  681. "ax[1].set_xlabel(r\"$E_\\gamma/E_0$\")\n",
  682. "ax[1].set_ylabel(\"counts (normed)\")\n",
  683. "ax[1].set_title(\"Downstream\")\n",
  684. "ax[1].legend()\n",
  685. "ax[1].grid()\n",
  686. "\n",
  687. "\"\"\"\n",
  688. "most electrons lose little energy relative to their initial energy downstream\n",
  689. "\"\"\"\n",
  690. "fig.suptitle(r\"$B\\rightarrow K^\\ast ee$, $p>5$GeV, photons w/ brem_vtx_z$<9500$mm\")\n",
  691. "plt.show()"
  692. ]
  693. },
  694. {
  695. "cell_type": "code",
  696. "execution_count": 20,
  697. "metadata": {},
  698. "outputs": [
  699. {
  700. "data": {
  701. "image/png": "iVBORw0KGgoAAAANSUhEUgAABjYAAAJOCAYAAAAUHj4bAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAADIAElEQVR4nOzdd3hUVcLH8d8lIYWQhISSEKkqIAJBBUXKSm+CKBYUMAIquqIiCq5dYFcBwbYLtrXhihRdBQXdCEixoggii7iWFQWUAEoIRSAknPcP3swypNwz5M4wQ76f55lHuXPm9Hvm3pw59zjGGCMAAAAAAAAAAIAIUOl4ZwAAAAAAAAAAAMAWExsAAAAAAAAAACBiMLEBAAAAAAAAAAAiBhMbAAAAAAAAAAAgYjCxAQAAAAAAAAAAIgYTGwAAAAAAAAAAIGIwsQEAAAAAAAAAACIGExsAAAAAAAAAACBiMLEBAAAAAAAAAAAiBhMbAAAAx9HatWu1cePG450NAAAAAAAiBhMbAAAAx9Hf//53LVmy5HhnAwAAAACAiMHEBgAAQIi9//77Gj58uLZt2+Y7tn79el188cXaunXrccwZAAAAAADhj4kNAAAQ0d555x05juP3SkpKUqtWrfTqq6+GJO1//OMffsdzc3PVu3dvxcTEaNq0acU+16pVK6WlpSkzM1Pvvfeepk2bpu7du6tXr16qWbNmUPMcSZYtW1asbYteK1asKPEza9eu1TXXXKNTTjlF8fHxio+PV6NGjXT99dfr888/DzgP/fv3V3x8vHbu3FlqmMGDB6ty5coBT0qNGzdOjuPo119/LfH95s2bq1OnTgHFGYiPP/5Y48aNK7NsFYFbOwSqItRrYWGhatWqpccee6zUMF7Xa6Q4kdv/s88+U8+ePZWYmKiqVauqc+fO+uijj4qFC2Ts3rNnj0aNGqWMjAzFxcXpjDPO0OzZs0tMP5CwAADgxBd9vDMAAABQHqtXr5Ykvfnmm6pVq5aMMdq4caPGjh2rgQMH6rTTTlNmZmZQ027VqpXv2Nq1a9W/f3/t3btXS5YsUYcOHYp9LiEhQX/+85916NAhTZw4UZUqVdIrr7yiK664Iij5jHQTJkxQ586d/Y41b968WLhnnnlGN910k5o0aaJbbrlFzZo1k+M4+vrrrzVr1iydffbZ+v7773XKKadYp33NNddo3rx5mjlzpkaMGFHs/by8PM2dO1d9+/ZVWlpa4IU7jj7++GONHz9eQ4cOVbVq1Y53dk4YFaFe33//fW3fvl0XX3zx8c5K2DlR23/lypU677zzdM455+jll1+WMUaTJ09W165dtXTpUrVt27bYZ2zG7osvvlgrV67UpEmT1LhxY82cOVMDBw7UoUOHNGjQoGMOCwAATnxMbAAAgIi2evVqJScnq1+/fr5jbdu2VUFBga688kp98cUXQZ3YiI+P12mnnSZJmj17tq655hplZmbq9ddfV0ZGRomf+/e//62srCylp6fr/PPPV0ZGhh588EG9/PLLeumll1SjRo2g5DeUduzYoUOHDnlSlkaNGuncc88tM8xHH32kESNGqE+fPvrnP/+pmJgY33tdunTRjTfeqNdee03x8fEBpd27d29lZGTohRdeKHFiY9asWdq3b5+uueaagOKNNL///ruqVKlyvLOBMPHPf/5TrVu3Vv369YOWBn0u9LZs2aKEhAQlJSUVe+++++5TtWrVlJ2d7WuXbt266eSTT9aYMWNKXLnhNna/8847WrRokW+CQpI6d+6sn376Sbfffrsuv/xyRUVFBRwWAABUDDyKCgAARLRVq1bpjDPOKHZ88+bNkqSmTZsGHOdjjz2mefPmBZT2mDFjNHDgQA0ePFjLly8vdVJDklJSUvTAAw8oOztbDRs2VPv27bVmzRoNGDBAycnJpX7uww8/VI8ePZScnKyUlBT16dNH33333TGHs3XBBReodevWevbZZ9WyZUvFx8erbt26Gjt2rA4dOlTiZ9auXavatWurd+/e+sc//qHdu3cfc/o2JkyYoKioKD3zzDN+kxpHuuyyy4q1y3fffadBgwapVq1aio2NVdOmTfXEE0/43o+KitKQIUO0atUq/fvf/y4W54svvugrZ7AVPdrniy++0MUXX6ykpCQlJyfryiuv1Pbt2/3Cbt++Xdddd53q1q2r2NhY1axZU+3bt9fixYt9cd1+++2SpIYNG/oeE7Ns2TJfOqtXr9all16qlJQUv1UubnUmSd9//72GDRumRo0aqUqVKjrppJN0wQUXFKvDorTWrl2ryy67TMnJyUpNTdVtt92mgoICffPNN+rVq5cSExPVoEEDTZ482fO6KrJ161YNHDhQycnJSktL09VXX628vDy/MB9++KG6du2qxMREValSRe3atdPbb7/tl25p9Wobx5H5/+qrr1zz5NbWJfnqq6/kOI5ee+0137FVq1bJcRw1a9bML2y/fv38VqUZYzR37lxdcsklpcZ/pE2bNrm2gRd9Llh9SZLmzZsnx3H03nvvFXvvqaee8qVbVvvv379fZ555pk499VS/NszJyVF6ero6deqkwsJCq/yU9pgnx3H0448/WpdLOvzoxOeff17dunVTnTp19MMPP5QY7qOPPlKnTp38JpsSExN13nnn6eOPP9aWLVsCSleS5s6dq6pVq+qyyy7zOz5s2DD98ssv+vTTT48pbHn7QjD7EgAA8A4TGwAAIGL99ttv2rhxo1q2bKmCggIVFBRo27Ztevnll/Xggw/q2muv1TnnnBNwvJ9//rkGDBhQ5uRGUdr16tVTjx49NG3aND377LP6+9//Xuof1ovUqVNHffv29TtW9Af0ypUrl/iZcePGqWPHjqpbt65mzZql5557Tps2bVLXrl21Z8+egMMFYtWqVfrPf/6jxx57TLfffrveeustdejQQX/+85/1wgsvlPiZc889Vy+99JIqV66s4cOHq1atWrrsssv0xhtv6MCBAwGlf+ONNyo6OlpJSUnq2bOnPvzwQ7/3CwsLtXTpUrVu3Vq1a9e2jnf9+vU6++yztW7dOj3yyCNasGCB+vTpo5EjR2r8+PG+cFdffbUcxylW1vXr1+uzzz7TkCFDQvpL4f79++vUU0/VP//5T40bN07z5s1Tz549dfDgQV+YrKwszZs3T/fff78WLlyo5557Tt26ddNvv/0mSbr22mt18803S5LeeOMNffLJJ/rkk0901lln+eK4+OKLdeqpp+q1117T008/7SuzTZ398ssvql69uiZNmqTs7Gw98cQTio6OVps2bfTNN98UK9OAAQPUsmVLvf766xo+fLgee+wx3XrrrbrooovUp08fzZ07V126dNEdd9yhN954w9O6KnLJJZeocePGev3113XnnXdq5syZuvXWW33vL1++XF26dFFeXp6ef/55zZo1S4mJibrgggs0Z84cq3q1iSOQPEnubV2SZs2aqXbt2n6TH4sXL1Z8fLzWr1+vX375RZJUUFCg5cuXq1u3br5wRX/Etp3YCKQNytPnigSjL/Xt21e1atXSiy++WOy96dOn66yzzlJmZmaZ7R8XF6dXX31V27Zt09VXXy1JOnTokAYPHixjjGbNmmU9jhTFW/RasmSJTjrpJKWnpys1NdX187///rvmzJmjCy+8UOnp6br55ptVrVo1zZkzR6effnqJn8nPz1dsbGyx40XHSpr4dRu7161bp6ZNmyo62v9BEkWrLNetW3dMYYuUty8Ea1wCAAAeMQAAABFq4cKFRlKxV3R0tHnggQeOOd6CggIzaNAgU7lyZTN37lzXtOPi4syKFSuOOT038+fPN5LM5MmT/Y5/++23RpKZMWNGQOECsXnzZiPJnHzyyWbnzp2+4/n5+SY9Pd307dvXNY7c3FzzwgsvmB49epjo6GiTnJxshg4dat59911TUFBQ6udWr15tbrnlFjN37lzz/vvvmxdeeME0bdrUREVFmezsbF+4nJwcI8lcccUVxeIoKCgwBw8e9L0OHTrke69nz56mTp06Ji8vz+8zN910k4mLizM7duzwHevYsaOpUaOGyc/P9x0bPXq0kWS+/fZb1zo
  702. "text/plain": [
  703. "<Figure size 2000x600 with 3 Axes>"
  704. ]
  705. },
  706. "metadata": {},
  707. "output_type": "display_data"
  708. }
  709. ],
  710. "source": [
  711. "#energyloss in abh von der energie der elektronen\n",
  712. "#upstream\n",
  713. "fig, ((ax0, ax1)) = plt.subplots(nrows=1, ncols=2, figsize=(20,6))\n",
  714. "\n",
  715. "a0=ax0.hist2d(up_energyloss_found, up_energy_found, bins=(np.linspace(0,1,80), np.linspace(0,5e4,80)), cmap=plt.cm.jet, cmin=1, vmax=15)\n",
  716. "ax0.set_ylim(0,5e4)\n",
  717. "ax0.set_xlim(0,1)\n",
  718. "ax0.set_xlabel(r\"energyloss $E_\\gamma/E_0$\")\n",
  719. "ax0.set_ylabel(r\"$E_0$\")\n",
  720. "ax0.set_title(\"found energyloss wrt electron energy\")\n",
  721. "\n",
  722. "a1=ax1.hist2d(up_energyloss_lost, up_energy_lost, bins=(np.linspace(0,1,50), np.linspace(0,5e4,50)), cmap=plt.cm.jet, cmin=1, vmax=15)\n",
  723. "ax1.set_ylim(0,5e4)\n",
  724. "ax1.set_xlim(0,1)\n",
  725. "ax1.set_xlabel(r\"energyloss $E_\\gamma/E_0$\")\n",
  726. "ax1.set_ylabel(r\"$E_0$\")\n",
  727. "ax1.set_title(\"lost energyloss wrt electron energy\")\n",
  728. "\n",
  729. "fig.colorbar(a1[3],ax=ax1)\n",
  730. "fig.suptitle(r\"$B\\rightarrow K^\\ast ee$, $p>5$GeV, Upstream photons w/ brem_vtx_z$<9500$mm\")\n",
  731. "\n",
  732. "\n",
  733. "plt.show()"
  734. ]
  735. },
  736. {
  737. "cell_type": "code",
  738. "execution_count": 17,
  739. "metadata": {},
  740. "outputs": [
  741. {
  742. "data": {
  743. "image/png": "iVBORw0KGgoAAAANSUhEUgAABj8AAAJOCAYAAADoCxXRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAADGZ0lEQVR4nOzdeXwURf7/8feQm5AMCUdCkFO5BAIKyrlyX3IoCChoBFRkRUUEPPACXDnk3gURZRWU2wNYATcCgiByyCGyiF/QFQWEAEoIh0BIUr8//GWWIcf0QCdMMq/n4zEPSPenq6qreibdqakqhzHGCAAAAAAAAAAAoJAocr0LAAAAAAAAAAAAYCc6PwAAAAAAAAAAQKFC5wcAAAAAAAAAAChU6PwAAAAAAAAAAACFCp0fAAAAAAAAAACgUKHzAwAAAAAAAAAAFCp0fgAAAAAAAAAAgEKFzg8AAAAAAAAAAFCo0PkBAAAAAAAAAAAKFTo/AAAA8tnu3bt18ODB610MAAAAAAAKLTo/AAAA8tnbb7+ttWvXXu9iAAAAAABQaNH5AQAAkA82bNig/v376/jx465te/fuVbdu3XTs2LHrWDIAAAAAAAofOj8AAECB8+mnn8rhcLi9IiMjVa9ePX3wwQf5kvf777/vtj05OVkdOnRQcHCwpk+fnuW4evXqKSYmRvHx8fr88881ffp0tWnTRu3bt1epUqXytMwFyRdffJGlbTNfW7ZsyfaY3bt36+GHH9aNN96osLAwhYWFqUqVKhowYIC2b9/udRm6du2qsLAwnTp1KseY+++/X0FBQV53XM2ZM8ftnEJDQxUbG6sWLVpo7Nixbp1jBdGmTZs0cuTIXOvOH4wcOVIOh0O//fabbWn6Q92mp6erdOnSmjJlSo4xeVG3BUFhbv+vv/5a7dq1U0REhIoVK6YWLVroq6++covx9nfD2bNnNXjwYMXFxSk0NFR169bVokWLrjoOAAAUTIHXuwAAAADe2rlzpyTpX//6l0qXLi1jjA4ePKgRI0aoV69eql69uuLj4/M073r16rm27d69W127dtW5c+e0du1aNW3aNMtx4eHhevXVV5WRkaGxY8eqSJEimj9/vu677748KWdBN2bMGLVo0cJtW61atbLEvfXWW3riiSdUrVo1PfXUU6pZs6YcDoe+//57LVy4ULfddpt+/PFH3XjjjZbzfvjhh7Vs2TItWLBAAwcOzLI/JSVFS5cuVadOnRQTE+P9yUmaPXu2qlevrkuXLun48ePauHGjXn/9dU2cOFGLFy9W69atryrd623Tpk0aNWqU+vbtq+LFi1/v4hQq/lC3GzZs0IkTJ9StW7frXRSfU1jbf9u2bbrjjjt0++23a+7cuTLGaPz48WrVqpXWrVunRo0aucVb/d3QrVs3bdu2TePGjVPVqlW1YMEC9erVSxkZGerdu7fXcQAAoGCi8wMAABQ4O3fulNPpVJcuXVzbGjVqpLS0ND3wwAP65ptv8rTzIywsTNWrV5ckLVq0SA8//LDi4+P18ccfKy4uLtvj/vOf/yghIUGxsbG68847FRcXp9GjR2vu3Ll67733VLJkyTwpb346efKkMjIybDmXKlWqqGHDhrnGfPXVVxo4cKA6duyojz76SMHBwa59LVu21OOPP64PP/xQYWFhXuXdoUMHxcXF6d13382282PhwoU6f/68Hn74Ya/SvVytWrVUv35918/33HOPnn76aTVt2lTdunXTDz/8cNUdKwXJH3/8oaJFi17vYsBHfPTRR6pfv74qVKiQJ+lzveW/o0ePKjw8XJGRkdnuf/nll1W8eHElJia62qZ169aqXLmyhg0blmUEiJXfDZ9++qlWr17t6siQpBYtWuiXX37RM888o3vvvVcBAQGW4wAAQMHFtFcAAKDA2bFjh+rWrZtl++HDhyVJNWrU8DrNKVOmaNmyZV7lPWzYMPXq1Uv333+/1q9fn2PHhyRFRUXptddeU2JioipVqqQmTZpo165d6tmzp5xOZ47Hbdy4UW3btpXT6VRUVJQ6duyoH3744arjrOrcubPq16+vWbNmqU6dOgoLC1O5cuU0YsQIZWRkZHvM7t27VaZMGXXo0EHvv/++zpw5c9X5WzFmzBgFBATorbfecuv4uFyPHj2ytMsPP/yg3r17q3Tp0goJCVGNGjX0xhtvuPYHBASoT58+2rFjh/7zn/9kSXP27Nmu87RT+fLlNWnSJJ05c0ZvvfWW276NGzeqVatWioiIUNGiRdW4cWOtXLnStf+7776Tw+HQhx9+6Nq2Y8cOORwO1axZ0y2tLl26uI1cypxG6LvvvlOvXr3kdDoVExOjhx56SCkpKa64EydO6NFHH1W5cuUUEhKiUqVKqUmTJlqzZo0rnWeeeUaSVKlSJdeUNF988YVbPjt37lT37t0VFRXlGpHjqU0y/fjjj+rXr5+qVKmiokWLqmzZsurcuXOWdsrMa/fu3erRo4ecTqeio6M1ZMgQpaWlad++fWrfvr0iIiJUsWJFjR8/3mP7ZKb5zTffqFu3boqMjJTT6dQDDzygEydOZHvMsWPHcq3TTJ7a11Pdejr+ynO41rbOzrVcg5JkjNHSpUt1zz335JjH5Q4dOpRrO+R2vUnWrrm8uI4ut2zZMjkcDn3++edZ9r355puuvHNq/8TERN1yyy266aab3NovKSlJsbGxat68udLT0y2VJacppRwOh37++Wevzis5OVnvvPOOWrdurRtuuEE//fRTjrFfffWVmjdv7tYpFRERoTvuuEObNm3S0aNHvcpbkpYuXapixYqpR48ebtv79eunI0eOaOvWrV7FSdd+LeT1tQQAALJH5wcAAChQfv/9dx08eFB16tRRWlqa0tLSdPz4cc2dO1ejR4/WI488ottvv93rdLdv366ePXvm2gGSmXf58uXVtm1bTZ8+XbNmzdLbb7+d4x/fM91www3q1KmT27bMP7IHBQVle8zIkSPVrFkzlStXTgsXLtQ///lPHTp0SK1atdLZs2e9jvPGjh079H//93+aMmWKnnnmGX3yySdq2rSpXn31Vb377rvZHtOwYUO99957CgoKUv/+/VW6dGn16NFDS5Ys0cWLF73K//HHH1dgYKAiIyPVrl07bdy40W1/enq61q1bp/r166tMmTKW0927d69uu+027dmzR5MmTdKKFSvUsWNHDRo0SKNGjXLFPfTQQ3I4HFnOde/evfr666/Vp0+fPPlG8J133qmAgABt2LDBtW39+vVq2bKlUlJS9M4772jhwoWKiIhQ586dtXjxYklSzZo1VaZMGbc/Tq9Zs0ZhYWHau3evjhw5IklKS0vT+vXrs51W65577lHVqlX18ccf6/nnn9eCBQv09NNPu/YnJCRo2bJleuWVV7Rq1Sr985//VOvWrfX7779Lkh555BE9+eSTkqQlS5Zo8+bN2rx5s2699Va3fLp166abbrpJH374oWbOnGm5TSTpyJEjKlGihMaNG6fExES98cYbCgwMVIMGDbRv374s59SzZ0/VqVNHH3/8sfr3768pU6bo6aef1t13362OHTtq6dKlatmypZ577jktWbLEUht17dpVN910kz766CONHDlSy5YtU7t27XTp0iWv61Sy1r651a2V470tl6e2zs61XoOZf+i22vlhtR2uvN4k658DmfLiOpKkTp06qXTp0po9e3aWfXPmzNGtt96q+Pj4HNu/cePG+uCDD3T8+HE99NBDkqSMjAzdf//9MsZo4cKFlj+nMtPMfK1du1Zly5ZVbGysoqOjPR7/xx9/aPHixbrrrrsUGxurJ598UsWLF9fixYt1880353hcamqqQkJCsmzP3HZlx6an3w2StGfPHtWoUUOBge4TXWSOCN2zZ49XcZe71mshr64lAACQAwMAAFCArFq1ykjK8goMDDSvvfbaVaeblpZmevfubYKCgszSpUs95h0aGmq2bNly1fl5snz5ciPJjB8/3m37/v37jSQzb948r+K8cfjwYSPJVK5c2Zw6dcq1PTU11cTGxppOnTp5TCM5Odm8++67pm3btiYwMNA4nU7Tt29f89lnn5m0tLQ
  744. "text/plain": [
  745. "<Figure size 2000x600 with 3 Axes>"
  746. ]
  747. },
  748. "metadata": {},
  749. "output_type": "display_data"
  750. }
  751. ],
  752. "source": [
  753. "#downstream\n",
  754. "fig, ((ax0, ax1)) = plt.subplots(nrows=1, ncols=2, figsize=(20,6))\n",
  755. "\n",
  756. "a0=ax0.hist2d(down_energyloss_found, down_energy_found, bins=(np.linspace(0,1,80), np.linspace(0,5e4,80)), cmap=plt.cm.jet, cmin=1, vmax=15)\n",
  757. "ax0.set_ylim(0,5e4)\n",
  758. "ax0.set_xlim(0,1)\n",
  759. "ax0.set_xlabel(r\"energyloss $E_\\gamma/E_0$\")\n",
  760. "ax0.set_ylabel(r\"$E_0$\")\n",
  761. "ax0.set_title(\"found energyloss wrt electron energy\")\n",
  762. "\n",
  763. "a1=ax1.hist2d(down_energyloss_lost, down_energy_lost, bins=(np.linspace(0,1,50), np.linspace(0,5e4,50)), cmap=plt.cm.jet, cmin=1, vmax=15)\n",
  764. "ax1.set_ylim(0,5e4)\n",
  765. "ax1.set_xlim(0,1)\n",
  766. "ax1.set_xlabel(r\"energyloss $E_\\gamma/E_0$\")\n",
  767. "ax1.set_ylabel(r\"$E_0$\")\n",
  768. "ax1.set_title(\"lost energyloss wrt electron energy\")\n",
  769. "\n",
  770. "fig.colorbar(a1[3],ax=ax1)\n",
  771. "fig.suptitle(r\"$B\\rightarrow K^\\ast ee$, $p>5$GeV, Downstream photons w/ brem_vtx_z$<9500$mm\")\n",
  772. "\n",
  773. "\n",
  774. "plt.show()"
  775. ]
  776. },
  777. {
  778. "cell_type": "code",
  779. "execution_count": 18,
  780. "metadata": {},
  781. "outputs": [
  782. {
  783. "data": {
  784. "image/png": "iVBORw0KGgoAAAANSUhEUgAABk4AAAJOCAYAAADxgPt3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAADT4UlEQVR4nOzdd3gU1f7H8c+mh5AsoSWEDiK9KB3UgHRpFoqgkSZyBUQE7NIsoIKAgu3nVUG6Si8iBBCkdxH0ClxBUAigQgJIIOX8/oDsdUnbCbumvV/Psw/szHfOnDkzmZ2zZ885NmOMEQAAAAAAAAAAAOSV3RkAAAAAAAAAAADIKWg4AQAAAAAAAAAAuI6GEwAAAAAAAAAAgOtoOAEAAAAAAAAAALiOhhMAAAAAAAAAAIDraDgBAAAAAAAAAAC4joYTAAAAAAAAAACA62g4AQAAAAAAAAAAuI6GEwAAAAAAAAAAgOtoOAEAAAAAAAAAALiOhhMAAJCj9O7dW9OnT8/ubMBF+/fv1/Hjx7M7GwAAAAAAuA0NJwAAAMiy//u//9O6deuyOxsAAAAAALgNDScAACBPmT9/vqpXr67AwEDZbDbt27cvu7OUoZUrV8pmszm9QkJCVLduXX3++ef/2P4/++wzp+Xnzp1Tu3bt5Ofnp2nTpjmt27hxo/r3768zZ844lv3www+6//77dfr0aY/nOTf55ptvUp3flNe2bdtSxe/fv1/9+vVTxYoVFRgYqMDAQFWqVEkDBgzQrl27spSH++67T4GBgTp//ny6MQ899JB8fX0tn78xY8bIZrPp999/T3N9jRo11KxZM0tpumrLli0aM2ZMhseVH2R2DrIiP5RtUlKSihcvrsmTJ6cb44myzQ3y8vnfsWOH2rRpo+DgYBUsWFDNmzfX5s2bnWKs3rcvXryooUOHKiIiQgEBAapTp47mzZuX5v6txAIAgOxFwwkAAMh2HTp0UKFChVSoUCHNmTNHAwcOdLx//fXXXU7n7NmzioqKUsWKFbVq1Spt3bpVt956qwdzfvP27NkjSVqyZIm2bt2qLVu26KOPPtKlS5fUo0cP7d+//x/Zf926dR3L9u/fr3r16mnv3r1at26dBg8e7LRN3bp1FRYWplq1amnt2rWaNm2aWrVqpbZt26pYsWIezW9uNW7cOG3dutXpVaNGDaeYDz/8UHXr1tX27dv15JNPavny5VqxYoWGDh2qgwcPqn79+vrvf/9red/9+vVTfHy85syZk+b62NhYLVq0SB06dFBYWFiWji87bNmyRWPHjs2TX+5mt/xQths3btTZs2d1//33Z3dWcpy8ev537typu+66S5cvX9bMmTM1c+ZMxcfHq0WLFtq6dWuqeFfu25J0//33a8aMGRo9erS++uor1a9fXz169EjznmslFgAAZC+f7M4AAADA8uXLHf/v3bu3mjVrpt69e1tO59ChQ0pISNDDDz+syMjIdOP++usvFShQICtZdbs9e/bIbrerU6dOjmWNGzdWYmKiHn74Ye3du1e1atXy6P4DAwNVpUoVSdK8efPUr18/1apVSwsWLFBERESqbYKCgvTyyy8rOTlZ48ePl5eXl2bPnq0HH3zQY/nMDn/++aeSk5NVtGjRm06rUqVKatSoUbrrN2/erIEDB6p9+/b68ssv5efn51h39913a9CgQfriiy8UGBhoed/t2rVTRESEPvnkEw0cODDV+rlz5+ry5cvq16+f5bRzk5z0d4/s9+WXX6pevXoqW7asR9Lnessep06dUlBQkEJCQlKtGzlypAoVKqRVq1Y5zk3Lli1VoUIFjRgxIlXPk8zu29K1Xptr1qzRnDlz1KNHD0lS8+bN9csvv+jpp59W9+7d5e3tbTkWAABkP3qcAAAAj9i0aZNat24tu92u0NBQtW/fXocPH/bY/nr37q077rhDktS9e3fZbDY1a9bMMdTKnj171KVLF4WGhqpixYpO+WzRooWCg4NVoEABNWnSRCtWrHBKOyWN/fv3q2vXrrLb7SpcuLCGDRumxMRE/fTTT2rbtq2Cg4NVrlw5vfnmmy7ne/fu3apTp06q5b/++qskqWrVqlkoDWny5MlavHixpf2PGDFCPXr00EMPPaQNGzak2WgiSd9//71uv/127dmzR/fcc4/69u2r1157Te3bt890SBtXrwt3Xz8dO3ZUvXr19NFHH6l27doKDAxU6dKlNXr0aCUnJ6e5zf79+1WiRAm1a9dOn332mS5cuJDl/Wdm3Lhx8vb21ocffujUaPJ3Xbt2TXVODh8+rJ49e6p48eLy9/dX1apV9e677zrFeHt7q1evXtq9e7e+//77VOl++umnjuP0tJS/pb179+r+++9XSEiI7Ha7Hn74YZ09e9YRd/bsWT322GMqXbq0/P39VaxYMTVt2lTR0dGOdJ5++mlJUvny5R3D6HzzzTdO+0nv796Vcjty5Ij69OmjSpUqqUCBAipZsqQ6duyYqgw9dX9wtaxSnD59Wj169JDdbldYWJj69u2r2NjYVHGZ3fMyK1tX7pl/z//BgwczzFdm5zotBw8elM1m0xdffOFYtnv3btlsNlWvXt0ptlOnTk496iTJGKNFixbpgQceSHcff3fixIkMz4E7rjdPfs5I0uLFi2Wz2bR27dpU695//33HvtM7/6tWrdJtt92mW265xen8xcTEKDw8XM2aNVNSUpLL+UlvKCybzaZjx45ZOrZz587p448/VsuWLVWqVCn9/PPPacZt3rxZzZo1c2rQCg4O1l133aUtW7bo1KlTlvYrSYsWLVLBggXVtWtXp+V9+vTRyZMntX37dsuxN3stePpaAgAg3zAAAABuNnr0aOPl5WX69u1rVqxYYb788ktTs2ZNU7p0aXPhwgWP7PPIkSPm3XffNZLMuHHjzNatW83BgwfN6NGjjSRTtmxZ8+yzz5o1a9aYxYsXG2OM+eabb4yvr6+pW7eumT9/vlm8eLFp3bq1sdlsZt68eU7HI8lUrlzZvPLKK2bNmjXmmWeeMZLM4MGDTZUqVcw777xj1qxZY/r06WMkmQULFmSa599//91IMkOGDDEJCQkmISHBnD592nz22WcmODjYPProo1kuj549expfX1+zaNGiTPffvXt3c/fddxt/f3/z0UcfZZr2iRMnzLJly4wxxgwaNMh8+umnJjEx0UyfPt1cvXo13e1cvS48cf2UKFHCBAUFmapVq5qZM2ea1atXmwcffNBISveYL1++bGbPnm06duxo/Pz8TEBAgOnSpYtZsGCBiY+Pd2m/69evN5JM8eLFjbe3twkODjatW7c23377rSMmMTHRBAYGmsaNG1s6poMHDxq73W5q1qxpPvvsM7N69WozfPhw4+XlZcaMGeMUe/jwYWOz2czQoUNTpSHJPPfcc5b2nSLlb+Ps2bNprq9evbqJjIxMFV+2bFnz9NNPm6+//tpMmjTJBAUFmdtuu81x/bRp08YUK1bM/N///Z/55ptvzOLFi82oUaMcf5cnTpwwTzzxhJFkFi5caLZu3Wq2bt1qYmNjU+3nxr97V8ttw4YNZvjw4ebLL780GzZsMIsWLTL33nuvCQwMNP/5z39SHZO77w+ultXf9z9q1CizZs0aM2nSJOPv72/69OnjlKYr97yMytbVe6aVfGV2rtNTokQJ89hjjznev/766yYwMNBIMr/99psxxpiEhAQTEhJinnnmGadtN23aZCSZQ4cOufUc3Mz15qnrKEVCQoIpXry4eeihh1Kta9Cggbn99tuNMRmf/0OHDpng4GBz//33G2OMSUpKMnfffbcpXry4OXnypMt5McY40k15rVu3zpQsWdKEh4c7/o4zcunSJTNv3jzTqVMn4+fnZwIDA80DDzxgvvjiC3PlypU0t/Hz8zOPPPJIquU9evQwkszXX39tjHHtvp2iUaNGpn79+qmWHzhwwEgyH374oeXYm70WPH0tAQCQX9BwAgAA3GrZsmVGknnzzTedlh86dMhIMrNmzUq1Tdu2bU1QUFCar9dee83lfad82fHFF184lqV8gTBq1KhU8Y0aNTLFixd3+jI+MTHR1KhRw5QqVco
  785. "text/plain": [
  786. "<Figure size 2000x600 with 3 Axes>"
  787. ]
  788. },
  789. "metadata": {},
  790. "output_type": "display_data"
  791. }
  792. ],
  793. "source": [
  794. "#plot residual energy against energyloss and try to find a good split (eg energyloss before and after the magnet)\n",
  795. "#upstream\n",
  796. "nbins=60\n",
  797. "\n",
  798. "fig, ((ax0, ax1)) = plt.subplots(nrows=1, ncols=2, figsize=(20,6))\n",
  799. "\n",
  800. "a0=ax0.hist2d(up_energyloss_found, up_residual_found, bins=(np.linspace(0,1,nbins), np.linspace(0,5e4,nbins)), cmap=plt.cm.jet, cmin=1, vmax=20)\n",
  801. "ax0.set_ylim(0,5e4)\n",
  802. "ax0.set_xlim(0,1)\n",
  803. "ax0.set_xlabel(r\"energyloss $E_\\gamma/E_0$\")\n",
  804. "ax0.set_ylabel(r\"$E_0-E_\\gamma$\")\n",
  805. "ax0.set_title(\"found energyloss wrt residual electron energy\")\n",
  806. "\n",
  807. "a1=ax1.hist2d(up_energyloss_lost, up_residual_lost, bins=(np.linspace(0,1,nbins), np.linspace(0,5e4,nbins)), cmap=plt.cm.jet, cmin=1, vmax=20) \n",
  808. "ax1.set_ylim(0,5e4)\n",
  809. "ax1.set_xlim(0,1)\n",
  810. "ax1.set_xlabel(r\"energyloss $E_\\gamma/E_0$\")\n",
  811. "ax1.set_ylabel(r\"$E_0-E_\\gamma$\")\n",
  812. "ax1.set_title(\"lost energyloss wrt residual electron energy\")\n",
  813. "\n",
  814. "fig.colorbar(a1[3],ax=ax1)\n",
  815. "fig.suptitle(r\"$e^\\pm$ from $B\\rightarrow K^\\ast ee$, $p>5$GeV, Upstream photons w/ brem_vtx_z$<9500$mm\")\n",
  816. "\n",
  817. "\"\"\"\n",
  818. "\"\"\"\n",
  819. "plt.show()"
  820. ]
  821. },
  822. {
  823. "cell_type": "code",
  824. "execution_count": 19,
  825. "metadata": {},
  826. "outputs": [
  827. {
  828. "data": {
  829. "image/png": "iVBORw0KGgoAAAANSUhEUgAABkEAAAJOCAYAAAAAi6D6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAADEXklEQVR4nOzdd3hUxf7H8c+mh5AsoYQQOkivCggIGpAuTRFU0EgTuQIiAnYlgAoqRb1gu14FpAgqTYp0BJHeRNAreAVBIYAYQpFAyvz+8Je9LGl7kl3YJO/X8+wDe/Z75syZc7I7s7MzYzPGGAEAAAAAAAAAAOQzPjc6AwAAAAAAAAAAAJ5AJwgAAAAAAAAAAMiX6AQBAAAAAAAAAAD5Ep0gAAAAAAAAAAAgX6ITBAAAAAAAAAAA5Et0ggAAAAAAAAAAgHyJThAAAAAAAAAAAJAv0QkCAAAAAAAAAADyJTpBAAAAAAAAAABAvkQnCAAAAAAAAAAAyJfoBAEAAB7Vp08fTZ8+/UZnAy7at2+fjh49eqOzAQAAAACAW9AJAgAAAId//etfWrdu3Y3OBgAAAAAAbkEnCAAA8Grz5s1TrVq1FBwcLJvNpr17997oLGVp+fLlstlsTo+wsDA1aNBAn3322XU7/ieffOK0PT4+Xh06dFBAQICmTp3q9NrGjRs1YMAAnTp1yrHthx9+ULdu3XTy5EmP5zkv+frrr9Nd37TH1q1b08Xv27dP/fv3V+XKlRUcHKzg4GBVqVJFAwcO1M6dO3OUh3vuuUfBwcE6e/ZspjEPPvig/P39LV+/6dOnO51TUFCQIiMj1bJlS40fP97pHsmLNm/erNGjR2dZdgXB6NGjZbPZ9Mcff7glvYJQrikpKYqIiNCbb76ZaYy7yzWvyM/Xf/v27WrXrp1CQ0NVuHBhtWzZUt9++226OCufDRcuXNCwYcMUFRWloKAg1a9fX3Pnzs3w+FZiAQBA5ugEAQAAbtepUycVKVJERYoU0Zw5czRo0CDH89dee83ldE6fPq2YmBhVrlxZK1as0JYtW1S1alUP5jz3du/eLUlavHixtmzZos2bN+vDDz/UxYsX1bNnT+3bt++6HL9BgwaObfv27VPDhg21Z88erVu3TkOGDHHap0GDBipZsqTq1q2rtWvXaurUqWrTpo3at2+vEiVKeDS/edW4ceO0ZcsWp0ft2rWdYj744AM1aNBA27Zt0xNPPKGlS5dq2bJlGjZsmA4cOKBGjRrpv//9r+Vj9+/fX4mJiZozZ06GryckJGjhwoXq1KmTSpYsmaPzmzZtmrZs2aLVq1frnXfeUf369fX666+rRo0aWrNmTY7S9AabN2/WmDFj8uWXtTdSQSjXjRs36vTp0+rWrduNzorXya/Xf8eOHbrjjjt06dIlzZw5UzNnzlRiYqJatWqlLVu2ZLiPK58N3bp104wZMxQbG6uvvvpKjRo1Us+ePTN8T7cSCwAAMud3ozMAAADyn6VLlzr+36dPH7Vo0UJ9+vSxnM7BgweVlJSkhx56SNHR0ZnG/fXXXypUqFBOsup2u3fvlt1uV5cuXRzbmjZtquTkZD300EPas2eP6tat69HjBwcHq3r16pKkuXPnqn///qpbt67mz5+vqKiodPuEhIRo7NixSk1N1fjx4+Xj46PZs2frgQce8Fg+b4Q///xTqampKl68eK7TqlKlipo0aZLp699++60GDRqkjh076osvvlBAQIDjtTvvvFODBw/W559/ruDgYMvH7tChg6KiovTxxx9r0KBB6V7/9NNPdenSJfXv399y2mlq166thg0bOp7fe++9evLJJ9W8eXN169ZNhw4dynEHS17hTe8ruPG++OILNWzYUOXLl/fYMbjnrr8TJ04oJCREYWFh6V576aWXVKRIEa1YscJxXVq3bq1KlSpp5MiRGY4Iye6zYfny5Vq9erXmzJmjnj17SpJatmypX3/9VU899ZTuv/9++fr6Wo4FAABZYyQIAABwyaZNm9S2bVvZ7XaFh4erY8eOOnTokMeO16dPHzVv3lySdP/998tms6lFixaO6UZ2796t7t27Kzw8XJUrV3bKZ6tWrRQaGqpChQrptttu07Jly5zSTktj37596tGjh+x2u4oWLarhw4crOTlZP/30k9q3b6/Q0FBVqFBBb7zxhsv53rVrl+rXr59u+2+//SZJqlGjRg5KQ3rzzTe1aNEiS8cfOXKkevbsqQcffFAbNmzIsANEkr7//nvdcsst2r17t+666y7169dPr776qjp27JjttC6u3hfuvn86d+6shg0b6sMPP1S9evUUHByssmXLKjY2VqmpqRnus2/fPpUqVUodOnTQJ598ovPnz+f4+NkZN26cfH199cEHHzh1gFytR48e6a7JoUOH1KtXL0VERCgwMFA1atTQO++84xTj6+ur3r17a9euXfr+++/TpTtt2jTHebpTuXLlNGnSJJ0/f14ffPCBY7srf3MHDhyQzWbT559/7ti2a9cu2Ww21apVyym2S5cujpFMaX+rBw4cUM+ePWW321WyZEn169dPCQkJTvudPn1ajz76qMqWLavAwECVKFFCzZo1c4xcGT16tJ566ilJUsWKFR1T1Xz99dfZvq+4cl1+/vln9e3bV1WqVFGhQoVUunRpde7cOd018tT7T1q6e/bsUbdu3RQWFia73a6HHnpIp0+fznCfkydPZluu2V3frMrV1TSuzr87rnVGcnoPSpIxRgsXLtS9996bafpXO3bsWLbXwB33nCc/yxYtWiSbzaa1a9eme+29995zHDer65+YmKibb75ZN910k9M1jIuLU2RkpFq0aKGUlBSX8pPZVFM2m01Hjhxx+bykv6eH/Oijj9S6dWuVKVNGv/zyS4Zx3377rVq0aOHUMRUaGqo77rhDmzdv1okTJywdV5IWLlyowoULq0ePHk7b+/btq+PHj2vbtm05is3tveDJewkAAK9gAAAAshEbG2t8fHxMv379zLJly8wXX3xh6tSpY8qWLWvOnz/vkWP+/PPP5p133jGSzLhx48yWLVvMgQMHTGxsrJFkypcvb5555hmzevVqs2jRImOMMV9//bXx9/c3DRo0MPPmzTOLFi0ybdu2NTabzcydO9fpfCSZatWqmZdfftmsXr3aPP3000aSGTJkiKlevbr55z//aVavXm369u1rJJn58+dnm+c//vjDSDJDhw41SUlJJikpyZw8edJ88sknJjQ01DzyyCM5Lo9evXoZf39/s3DhwmyPf//995s777zTBAYGmg8//DDbtI8dO2aWLFlijDFm8ODBZtq0aSY5OdlMnz7dXLlyJdP9XL0vPHH/lCpVyoSEhJgaNWqYmTNnmlWrVpkHHnjASMr0nC9dumRmz55tOnfubAICAkxQUJDp3r27mT9/vklMTHTpuOvXrzeSTEREhPH19TWhoaGmbdu25ptvvnHEJCcnm+DgYNO0aVNL53TgwAFjt9tNnTp1zCeffGJWrVplRowYYXx8fMzo0aOdYg8dOmRsNpsZNmxYujQkmWeffdbSsdNMmzbNSDI7duzI8PULFy4YX19f06pVK2OM639zxvx9zR599FHH89dee80EBwcbSeb33383xhiTlJRkwsLCzNNPP22Mcf5bHTVqlFm9erWZPHmyCQwMNH379nVKv127dqZEiRLmX//6l/n666/NokWLzKhRoxz5OHbsmHn88ceNJLNgwQKzZcsWs2XLFpOQkJDl+4qr12XDhg1mxIgR5osvvjAbNmwwCxcuNHfffbcJDg42//nPfxxxnnr/ufocnnrqKbNy5UozefJkExISYm6++Wanv2VXy9WV65tVuVq5R9x5rTOTk3vQGGM2bdpkJJmDBw+6/Rrk5p7z1L2UVg4RERHmwQcfTPfarbfeam655RZjTPbX/+DBgyY0NNR069bNGGNMSkqKufPOO01ERIQ5fvy4S3kxxjjSTXusW7fOlC5d2kRGRjqOlZWLFy+auXPnmi5dupiAgAATHBxs7r33XvP555+by5cvZ7hPQECAefjhh9Nt79mzp5FkVq5c6djmymeDMcY0adLENGrUKF2a+/f
  830. "text/plain": [
  831. "<Figure size 2000x600 with 3 Axes>"
  832. ]
  833. },
  834. "metadata": {},
  835. "output_type": "display_data"
  836. }
  837. ],
  838. "source": [
  839. "#downstream\n",
  840. "fig, ((ax0, ax1)) = plt.subplots(nrows=1, ncols=2, figsize=(20,6))\n",
  841. "\n",
  842. "a0=ax0.hist2d(down_energyloss_found, down_residual_found, bins=(np.linspace(0,1,60), np.linspace(0,5e4,60)), cmap=plt.cm.jet, cmin=1, vmax=25)\n",
  843. "ax0.set_ylim(0,5e4)\n",
  844. "ax0.set_xlim(0,1)\n",
  845. "ax0.set_xlabel(r\"energyloss $E_\\gamma/E_0$\")\n",
  846. "ax0.set_ylabel(r\"$E_0-E_\\gamma$\")\n",
  847. "ax0.set_title(\"found energyloss wrt residual electron energy\")\n",
  848. "\n",
  849. "a1=ax1.hist2d(down_energyloss_lost, down_residual_lost, bins=(np.linspace(0,1,60), np.linspace(0,5e4,60)), cmap=plt.cm.jet, cmin=1, vmax=20) \n",
  850. "ax1.set_ylim(0,5e4)\n",
  851. "ax1.set_xlim(0,1)\n",
  852. "ax1.set_xlabel(r\"energyloss $E_\\gamma/E_0$\")\n",
  853. "ax1.set_ylabel(r\"$E_0-E_\\gamma$\")\n",
  854. "ax1.set_title(\"lost energyloss wrt residual electron energy\")\n",
  855. "\n",
  856. "fig.colorbar(a0[3],ax=ax1)\n",
  857. "fig.suptitle(r\"$e^\\pm$ from $B\\rightarrow K^\\ast ee$, $p>5$GeV, Downstream photons w/ brem_vtx_z$<9500$mm\")\n",
  858. "\n",
  859. "\"\"\"\n",
  860. "\"\"\"\n",
  861. "plt.show()"
  862. ]
  863. },
  864. {
  865. "cell_type": "code",
  866. "execution_count": null,
  867. "metadata": {},
  868. "outputs": [],
  869. "source": []
  870. }
  871. ],
  872. "metadata": {
  873. "kernelspec": {
  874. "display_name": "env1",
  875. "language": "python",
  876. "name": "python3"
  877. },
  878. "language_info": {
  879. "codemirror_mode": {
  880. "name": "ipython",
  881. "version": 3
  882. },
  883. "file_extension": ".py",
  884. "mimetype": "text/x-python",
  885. "name": "python",
  886. "nbconvert_exporter": "python",
  887. "pygments_lexer": "ipython3",
  888. "version": "3.11.5"
  889. }
  890. },
  891. "nbformat": 4,
  892. "nbformat_minor": 2
  893. }