Projektpraktikum/B_rework.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [],
   "source": [
    "import uproot\t\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",
    "from mpl_toolkits.axes_grid1 import ImageGrid\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10522"
      ]
     },
     "execution_count": 107,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "file = uproot.open(\"tracking_losses_ntuple_Bd2KstEE.root:PrDebugTrackingLosses.PrDebugTrackingTool/Tuple;1\")\n",
    "\n",
    "#selektiere nur elektronen von B->K*ee und nur solche mit einem momentum von ueber 5 GeV \n",
    "allcolumns = file.arrays()\n",
    "found = allcolumns[(allcolumns.isElectron) & (~allcolumns.lost) & (allcolumns.fromSignal) & (allcolumns.p > 5e3)] #B: 9056\n",
    "lost = allcolumns[(allcolumns.isElectron) & (allcolumns.lost) & (allcolumns.fromSignal) & (allcolumns.p > 5e3)] #B: 1466\n",
    "\n",
    "ak.num(found, axis=0) + ak.num(lost, axis=0)\n",
    "#ak.count(found, axis=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "eff all =  0.8606728758791105 +/- 0.003375885792719708\n"
     ]
    }
   ],
   "source": [
    "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",
    "\n",
    "print(\"eff all = \", t_eff(found, lost), \"+/-\", eff_err(found, lost))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'energy': 46180.704276008204,\n",
       " 'brem_photons_pe': [3264.454345703125,\n",
       "  4453.86376953125,\n",
       "  178.029052734375,\n",
       "  14471.001953125,\n",
       "  1095.5640869140625,\n",
       "  3793.752685546875,\n",
       "  357.2669982910156,\n",
       "  825.275634765625,\n",
       "  8990.57421875,\n",
       "  3479.052490234375],\n",
       " 'brem_vtx_z': [161.9427032470703,\n",
       "  186.9705047607422,\n",
       "  387.3406982421875,\n",
       "  486.6791076660156,\n",
       "  1340.39501953125,\n",
       "  2322.552490234375,\n",
       "  9494.216796875,\n",
       "  12068.0263671875,\n",
       "  12118.072265625,\n",
       "  12129.564453125]}"
      ]
     },
     "execution_count": 115,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#try excluding all photons that originate from a vtx @ z>9500mm\n",
    "\n",
    "brem_e_f = found[\"brem_photons_pe\"]\n",
    "brem_z_f = found[\"brem_vtx_z\"]\n",
    "e_f = found[\"energy\"]\n",
    "\n",
    "\n",
    "\n",
    "brem_f = ak.ArrayBuilder()\n",
    "\n",
    "for itr in range(4):#range(ak.num(found,axis=0)):\n",
    "    brem_f.begin_record()\n",
    "    #[:,0] energy\n",
    "    brem_f.field(\"energy\").append(e_f[itr])\n",
    "    #[:,1] photon energy \n",
    "    brem_f.field(\"brem_photons_pe\").append(brem_e_f[itr])\n",
    "    #[:,2] 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",
    "#this is a event cut! not suitable\n",
    "cut = (brem_f[\"brem_vtx_z\"]<9500)\n",
    "brem_f[0].tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[{'energy': 134321.90632702777,\n",
       "  'brem_photons_pe': [1881.1756591796875,\n",
       "   17914.712890625,\n",
       "   479.9449768066406,\n",
       "   3987.47119140625,\n",
       "   4792.82421875,\n",
       "   19725.302734375,\n",
       "   2376.974853515625,\n",
       "   6870.6201171875,\n",
       "   19237.05859375,\n",
       "   3409.642822265625],\n",
       "  'brem_vtx_z': [184.35940551757812,\n",
       "   190.19970703125,\n",
       "   929.2517700195312,\n",
       "   5855.25390625,\n",
       "   8506.958984375,\n",
       "   11310.3232421875,\n",
       "   12051.8232421875,\n",
       "   12121.9033203125,\n",
       "   12132.9287109375,\n",
       "   12201.8427734375]},\n",
       " {'energy': 134321.90632702777,\n",
       "  'brem_photons_pe': [1881.1756591796875,\n",
       "   17914.712890625,\n",
       "   479.9449768066406,\n",
       "   3987.47119140625,\n",
       "   4792.82421875,\n",
       "   19725.302734375,\n",
       "   2376.974853515625,\n",
       "   6870.6201171875,\n",
       "   19237.05859375,\n",
       "   3409.642822265625],\n",
       "  'brem_vtx_z': [184.35940551757812,\n",
       "   190.19970703125,\n",
       "   929.2517700195312,\n",
       "   5855.25390625,\n",
       "   8506.958984375,\n",
       "   11310.3232421875,\n",
       "   12051.8232421875,\n",
       "   12121.9033203125,\n",
       "   12132.9287109375,\n",
       "   12201.8427734375]},\n",
       " {'energy': 134321.90632702777,\n",
       "  'brem_photons_pe': [1881.1756591796875,\n",
       "   17914.712890625,\n",
       "   479.9449768066406,\n",
       "   3987.47119140625,\n",
       "   4792.82421875,\n",
       "   19725.302734375,\n",
       "   2376.974853515625,\n",
       "   6870.6201171875,\n",
       "   19237.05859375,\n",
       "   3409.642822265625],\n",
       "  'brem_vtx_z': [184.35940551757812,\n",
       "   190.19970703125,\n",
       "   929.2517700195312,\n",
       "   5855.25390625,\n",
       "   8506.958984375,\n",
       "   11310.3232421875,\n",
       "   12051.8232421875,\n",
       "   12121.9033203125,\n",
       "   12132.9287109375,\n",
       "   12201.8427734375]},\n",
       " {'energy': 134321.90632702777,\n",
       "  'brem_photons_pe': [1881.1756591796875,\n",
       "   17914.712890625,\n",
       "   479.9449768066406,\n",
       "   3987.47119140625,\n",
       "   4792.82421875,\n",
       "   19725.302734375,\n",
       "   2376.974853515625,\n",
       "   6870.6201171875,\n",
       "   19237.05859375,\n",
       "   3409.642822265625],\n",
       "  'brem_vtx_z': [184.35940551757812,\n",
       "   190.19970703125,\n",
       "   929.2517700195312,\n",
       "   5855.25390625,\n",
       "   8506.958984375,\n",
       "   11310.3232421875,\n",
       "   12051.8232421875,\n",
       "   12121.9033203125,\n",
       "   12132.9287109375,\n",
       "   12201.8427734375]},\n",
       " {'energy': 134321.90632702777,\n",
       "  'brem_photons_pe': [1881.1756591796875,\n",
       "   17914.712890625,\n",
       "   479.9449768066406,\n",
       "   3987.47119140625,\n",
       "   4792.82421875,\n",
       "   19725.302734375,\n",
       "   2376.974853515625,\n",
       "   6870.6201171875,\n",
       "   19237.05859375,\n",
       "   3409.642822265625],\n",
       "  'brem_vtx_z': [184.35940551757812,\n",
       "   190.19970703125,\n",
       "   929.2517700195312,\n",
       "   5855.25390625,\n",
       "   8506.958984375,\n",
       "   11310.3232421875,\n",
       "   12051.8232421875,\n",
       "   12121.9033203125,\n",
       "   12132.9287109375,\n",
       "   12201.8427734375]},\n",
       " None,\n",
       " None,\n",
       " None,\n",
       " None,\n",
       " None]"
      ]
     },
     "execution_count": 116,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msa = ak.mask(brem_f,cut)\n",
    "msa[2].tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#ignore all brem vertices @ z>9500mm \n",
    "\n",
    "\"\"\"\n",
    "ph_e = found[\"brem_photons_pe\"]\n",
    "event_cut = ak.all(ph_e<cutoff_energy,axis=1)\n",
    "ph_e = ph_e[event_cut]\n",
    "\"\"\"\n",
    "\n",
    "brem_f = found[\"brem_photons_pe\"]\n",
    "brem_vtx_f = found[\"brem_vtx_z\"]\n",
    "event_cut = ak.any(brem_vtx_f<9500,axis=1)\n",
    "brems = found[event_cut]\n",
    "brems[0].tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#finden wir die elektronen die keine bremsstrahlung gemacht haben mit hoher effizienz?\n",
    "#von energie der photonen abmachen\n",
    "#scan ab welcher energie der photonen die effizienz abfällt\n",
    "\n",
    "#abhängigkeit vom ort der emission untersuchen   <- noch nicht gemacht\n",
    "\n",
    "\n",
    "\n",
    "#idea: we make an event cut st all events that contain a photon of energy > cutoff_energy are not included\n",
    "\"\"\"\n",
    "ph_e = found[\"brem_photons_pe\"]\n",
    "event_cut = ak.all(ph_e<cutoff_energy,axis=1)\n",
    "ph_e = ph_e[event_cut]\n",
    "\"\"\"\n",
    "\n",
    "efficiencies_found = ak.ArrayBuilder()\n",
    "\n",
    "\n",
    "\n",
    "for cutoff_energy in range(0,1050,25):\n",
    "\tnobrem_f = found[ak.all(found[\"brem_photons_pe\"]<cutoff_energy,axis=1)]\n",
    "\tnobrem_l = lost[ak.all(lost[\"brem_photons_pe\"]<cutoff_energy,axis=1)]\n",
    "\n",
    "\n",
    "\n",
    "\tprint(\"sample size: \",ak.num(nobrem_f,axis=0)+ak.num(nobrem_l,axis=0))\n",
    "\tprint(\"eff (cutoff = \",str(cutoff_energy),\") = \",np.round(t_eff(nobrem_f,nobrem_l),4), \"+/-\", np.round(eff_err(nobrem_f, nobrem_l),4))\n",
    "\n",
    "\"\"\"\n",
    "we see that a cutoff energy of 350MeV is ideal because the efficiency drops significantly for higher values\n",
    "\"\"\"\n",
    "cutoff_energy = 350.0 #MeV\n",
    "\n",
    "\"\"\"\n",
    "better statistics: cutoff=350MeV - sample size: 150 events and efficiency=0.9533\n",
    "\"\"\"\n",
    "nobrem_found = found[ak.all(found[\"brem_photons_pe\"]<cutoff_energy,axis=1)]\n",
    "nobrem_lost = lost[ak.all(lost[\"brem_photons_pe\"]<cutoff_energy,axis=1)]\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",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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",
    "#if any photon of an electron has an energy higher the cutoff then it is included\n",
    "cutoff_energy=350\n",
    "\n",
    "brem_found = found[ak.any(found[\"brem_photons_pe\"]>=cutoff_energy,axis=1)]\n",
    "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",
    "residual_found = energy_found - eph_found\n",
    "energyloss_found = eph_found/energy_found\n",
    "\n",
    "brem_lost = lost[ak.any(lost[\"brem_photons_pe\"]>=cutoff_energy,axis=1)]\n",
    "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",
    "residual_lost = energy_lost - eph_lost\n",
    "energyloss_lost = eph_lost/energy_lost\n",
    "\n",
    "print(\"eff = \", np.round(t_eff(brem_found,brem_lost),4), \"+/-\", np.round(eff_err(brem_found, brem_lost),4))\n",
    "brem_lost"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_energyloss_found = ak.mean(energyloss_found)\n",
    "mean_energyloss_lost = ak.mean(energyloss_lost)\n",
    "print(\"mean energyloss relative to initial energy (found): \", mean_energyloss_found)\n",
    "print(\"mean energyloss relative to initial energy (lost): \", mean_energyloss_lost)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#in abhängigkeit von der energie der elektronen\n",
    "plt.hist(energyloss_lost, bins=200, 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",
    "plt.xticks(np.arange(0,1.1,0.1), minor=True,)\n",
    "plt.yticks(np.arange(0,10,1), minor=True)\n",
    "plt.xlabel(r\"$E_\\gamma/E_0$\")\n",
    "plt.ylabel(\"counts (normed)\")\n",
    "plt.title(r'$E_{ph}/E_0$')\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#energyloss in abh von der energie der elektronen\n",
    "fig, ((ax0, ax1)) = plt.subplots(nrows=1, ncols=2, figsize=(20,6))\n",
    "\n",
    "a0=ax0.hist2d(energyloss_found, energy_found, bins=[100,500], cmap=plt.cm.jet, cmin=1, vmax=10)\n",
    "ax0.set_ylim(0,1e5)\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=[100,500], cmap=plt.cm.jet, cmin=1, vmax=10) \n",
    "ax1.set_ylim(0,1e5)\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\"$e^\\pm$ from $B\\rightarrow K^\\ast ee$, $p>5$GeV\")\n",
    "\n",
    "\"\"\"\n",
    "we can see that high energy electrons are often found even though they emit a lot of their energy through bremsstrahlung\n",
    "\"\"\"\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot residual energy against energyloss and try to find a good split (eg energyloss before and after the magnet)\n",
    "fig, ((ax0, ax1)) = plt.subplots(nrows=1, ncols=2, figsize=(20,6))\n",
    "\n",
    "a0=ax0.hist2d(energyloss_found, residual_found, bins=[100,500], cmap=plt.cm.jet, cmin=1, vmax=40)\n",
    "ax0.set_ylim(0,0.5e5)\n",
    "ax0.set_xlim(0,1)\n",
    "ax0.set_xlabel(r\"energyloss $E_\\gamma/E_0$\")\n",
    "ax0.set_ylabel(r\"$E_0-E_\\gamma$\")\n",
    "ax0.set_title(\"found energyloss wrt residual electron energy\")\n",
    "\n",
    "a1=ax1.hist2d(energyloss_lost, residual_lost, bins=[100,500], cmap=plt.cm.jet, cmin=1, vmax=40) \n",
    "ax1.set_ylim(0,0.5e5)\n",
    "ax1.set_xlim(0,1)\n",
    "ax1.set_xlabel(r\"energyloss $E_\\gamma/E_0$\")\n",
    "ax1.set_ylabel(r\"$E_0-E_\\gamma$\")\n",
    "ax1.set_title(\"lost energyloss wrt residual electron energy\")\n",
    "\n",
    "fig.colorbar(a1[3],ax=ax1)\n",
    "fig.suptitle(r\"$e^\\pm$ from $B\\rightarrow K^\\ast ee$, $p>5$GeV\")\n",
    "\n",
    "\"\"\"\n",
    "\"\"\"\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#try to find a split between energy lost before and after the magnet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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",
   "execution_count": null,
   "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",
    "#should be fulfilled by all candidates\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",
    "sf_energy_found = ak.to_numpy(scifi_found[\"energy\"])\n",
    "sf_eph_found = ak.to_numpy(ak.sum(scifi_found[\"brem_photons_pe\"], axis=-1, keepdims=False))\n",
    "sf_vtx_type_found = scifi_found[\"all_endvtx_types\"]\n",
    "\n",
    "\n",
    "brem_vtx_type_found = scifi_found[scifi_found[\"endvtx_type\"]==101]\n",
    "\n",
    "sf_energy_lost = ak.to_numpy(scifi_lost[\"energy\"])\n",
    "sf_eph_lost = ak.to_numpy(ak.sum(scifi_lost[\"brem_photons_pe\"], axis=-1, keepdims=False))\n",
    "sf_vtx_type_lost = scifi_lost[\"all_endvtx_types\"]\n",
    "brem_vtx_type_lost = scifi_lost[scifi_lost[\"endvtx_type\"]==101]\n",
    "\n",
    "\n",
    "\n",
    "#ak.num(scifi_found[\"energy\"], axis=0)\n",
    "#scifi_found.snapshot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ak.num(scifi_found[\"energy\"], axis=0)\n",
    "scifi_found[\"all_endvtx_types\"][1,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scifi_fitpars_found = ak.ArrayBuilder()\n",
    "vtx_types_found = ak.ArrayBuilder()\n",
    "\n",
    "for i in range(0,ak.num(scifi_found, 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",
    "    #[:,4] -> energy \n",
    "    scifi_fitpars_found.real(sf_energy_found[i])\n",
    "    #[:,5] -> photon energy\n",
    "    scifi_fitpars_found.real(sf_eph_found[i])\n",
    "    scifi_fitpars_found.end_list()\n",
    "    \n",
    "    vtx_types_found.begin_list()\n",
    "    #[:,0] -> endvtx_type\n",
    "    vtx_types_found.extend(sf_vtx_type_found[i,:])\n",
    "    vtx_types_found.end_list()\n",
    "    \n",
    "\n",
    "scifi_fitpars_lost = ak.ArrayBuilder()\n",
    "vtx_types_lost = ak.ArrayBuilder()\n",
    "\n",
    "for i in range(0,ak.num(scifi_lost, 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",
    "    #[:,4] -> energy \n",
    "    scifi_fitpars_lost.real(sf_energy_lost[i])\n",
    "    #[:,5] -> photon energy\n",
    "    scifi_fitpars_lost.real(sf_eph_lost[i])\n",
    "    scifi_fitpars_lost.end_list()\n",
    "    \n",
    "    vtx_types_lost.begin_list()\n",
    "    #endvtx_type\n",
    "    vtx_types_lost.extend(sf_vtx_type_lost[i,:])\n",
    "    vtx_types_lost.end_list()\n",
    "    \n",
    "\n",
    "\n",
    "scifi_fitpars_lost = ak.to_numpy(scifi_fitpars_lost)\n",
    "scifi_fitpars_found = ak.to_numpy(scifi_fitpars_found)\n",
    "\n",
    "vtx_types_lost = ak.Array(vtx_types_lost)\n",
    "vtx_types_found = ak.Array(vtx_types_found)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vtx_types_found[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#b parameter des fits [:,1] hat für lost eine breitere Verteilung. Warum?\n",
    "#evtl multiple scattering candidates (lost); findet man einen gewissen endvtx_type (mult scattering)\n",
    "#steiler velo winkel (eta)? vertex type? evtl bremsstrahlung?\n",
    "\n",
    "#isolate b parameters for analysis\n",
    "b_found = scifi_fitpars_found[:,1]\n",
    "b_lost = scifi_fitpars_lost[:,1]\n",
    "\n",
    "brem_energy_found = scifi_fitpars_found[:,5]\n",
    "brem_energy_lost = scifi_fitpars_lost[:,5]\n",
    "\n",
    "\n",
    "bs_found, vtx_types_found = ak.broadcast_arrays(b_found, vtx_types_found)\n",
    "bs_found = ak.to_numpy(ak.ravel(bs_found))\n",
    "vtx_types_found = ak.to_numpy(ak.ravel(vtx_types_found))\n",
    "\n",
    "bs_lost, vtx_types_lost = ak.broadcast_arrays(b_lost, vtx_types_lost)\n",
    "bs_lost = ak.to_numpy(ak.ravel(bs_lost))\n",
    "vtx_types_lost = ak.to_numpy(ak.ravel(vtx_types_lost))\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "#Erste Annahme ist Bremsstrahlung\n",
    "\n",
    "fig = plt.figure(figsize=(18,6))\n",
    "axes = ImageGrid(fig, 111,  # similar to subplot(111)\n",
    "                 nrows_ncols=(1, 2),  # creates 2x2 grid of axes\n",
    "                 axes_pad=1,  # pad between axes in inch.\n",
    "                 cbar_mode=\"single\",\n",
    "                 cbar_location=\"right\",\n",
    "                 cbar_pad=0.1,\n",
    "                 aspect=False\n",
    "                )\n",
    "\n",
    "\n",
    "h0 = axes[0].hist2d(b_found, brem_energy_found, bins=200, cmap=plt.cm.jet, cmin=1,vmax=30)\n",
    "axes[0].set_xlim(-1,1)\n",
    "axes[0].set_xlabel(\"b parameter [mm]\")\n",
    "axes[0].set_ylabel(r\"$E_{ph}$\")\n",
    "axes[0].set_title(\"found photon energy wrt b parameter\")\n",
    "\n",
    "h1 = axes[1].hist2d(b_lost, brem_energy_lost, bins=200, cmap=plt.cm.jet, cmin=1,vmax=30)\n",
    "axes[1].set_xlim(-1,1)\n",
    "axes[1].set_xlabel(\"b parameter [mm]\")\n",
    "axes[1].set_ylabel(r\"$E_{ph}$\")\n",
    "axes[1].set_title(\"lost photon energy wrt b parameter\")\n",
    "\n",
    "fig.colorbar(h0[3], cax=axes.cbar_axes[0], orientation='vertical')\n",
    "\n",
    "\"\"\"\n",
    "\"\"\"\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(18,6))\n",
    "\n",
    "a0=ax[0].hist2d(bs_found, vtx_types_found, bins=110, density=True, cmap=plt.cm.jet, cmin=1e-20,vmax=2)\n",
    "ax[0].set_ylim(0,110)\n",
    "ax[0].set_xlim(-1,1)\n",
    "ax[0].set_xlabel(\"b\")\n",
    "ax[0].set_ylabel(\"endvtx id\")\n",
    "ax[0].set_title(\"found endvtx id wrt b parameter\")\n",
    "ax[0].set_yticks(np.arange(0,110,1),minor=True)\n",
    "\n",
    "a1=ax[1].hist2d(bs_lost, vtx_types_lost, bins=110, density=True, cmap=plt.cm.jet, cmin=1e-20,vmax=2)\n",
    "ax[1].set_ylim(0,110)\n",
    "ax[1].set_xlim(-1,1)\n",
    "ax[1].set_xlabel(\"b\")\n",
    "ax[1].set_ylabel(\"endvtx id\")\n",
    "ax[1].set_title(\"lost endvtx id wrt b paraneter\")\n",
    "ax[1].set_yticks(np.arange(0,110,1), minor=True)\n",
    "\n",
    "\"\"\"\n",
    "vtx_id: 101 - Bremsstrahlung\n",
    "B:\n",
    "wir können nicht wirklich sagen dass bei den lost teilchen jegliche endvertex types überwiegen, im gegensatz zu den found \n",
    "\"\"\"\n",
    "fig.colorbar(a0[3], ax=ax, orientation='vertical')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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_xticks(np.arange(-1,1,0.1),minor=True)\n",
    "ax1.set_xlabel(\"b\")\n",
    "ax1.set_ylabel(\"normed\")\n",
    "ax1.set_title(\"fitparameter b der scifi track\")\n",
    "ax1.legend()\n",
    "#evtl multiple scattering candidates (lost); findet man einen gewissen endvtx_type (mult scattering)\n",
    "#steiler velo winkel (eta)? vertex type? evtl bremsstrahlung?\n",
    "\n",
    "\n",
    "ax2.hist(scifi_fitpars_found[:,2], bins=500, density=True, alpha=0.5, histtype='bar', color=\"blue\", label=r\"$c_x$ found\")\n",
    "ax2.hist(scifi_fitpars_lost[:,2], bins=500, density=True, alpha=0.5, histtype='bar', color=\"darkorange\", label=r\"$c_x$ lost\")\n",
    "ax2.set_xlim([-3e-5,3e-5])\n",
    "ax2.set_xticks(np.arange(-3e-5,3.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=500, density=True, alpha=0.5, histtype='bar', color=\"blue\", label=r\"$d_x$ found\")\n",
    "ax3.hist(scifi_fitpars_lost[:,3], bins=500, 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",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": []
  },
  {
   "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
}