Projektpraktikum/notebooks/B_updown.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 8,
   "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": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10522"
      ]
     },
     "execution_count": 9,
     "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": 10,
   "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": 11,
   "metadata": {},
   "outputs": [],
   "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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "cut_brem_found = ak.ArrayBuilder()\n",
    "\n",
    "for itr in range(ak.num(brem_f, axis=0)):\n",
    "    cut_brem_found.begin_record()\n",
    "    cut_brem_found.field(\"energy\").real(brem_f[itr,\"energy\"])\n",
    "    \n",
    "    cut_brem_found.field(\"brem_photons_pe\")\n",
    "    cut_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",
    "            cut_brem_found.real(brem_f[itr,\"brem_photons_pe\", jentry])\n",
    "            \n",
    "            #cut_brem_found.field(\"brem_vtx_z\").real(brem_f[itr, \"brem_vtx_z\",jentry])\n",
    "    cut_brem_found.end_list()\n",
    "    \n",
    "    cut_brem_found.field(\"brem_vtx_z\")\n",
    "    cut_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",
    "            cut_brem_found.real(brem_f[itr, \"brem_vtx_z\",jentry])\n",
    "    cut_brem_found.end_list()\n",
    "    \n",
    "\n",
    "    \n",
    "    cut_brem_found.end_record()\n",
    "\n",
    "cut_brem_found = ak.Array(cut_brem_found)\n",
    "\n",
    "\n",
    "\n",
    "cut_brem_lost = ak.ArrayBuilder()\n",
    "\n",
    "for itr in range(ak.num(brem_l, axis=0)):\n",
    "    cut_brem_lost.begin_record()\n",
    "    cut_brem_lost.field(\"energy\").real(brem_l[itr,\"energy\"])\n",
    "    \n",
    "    \n",
    "    cut_brem_lost.field(\"brem_photons_pe\")\n",
    "    cut_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",
    "            cut_brem_lost.real(brem_l[itr,\"brem_photons_pe\", jentry])\n",
    "            \n",
    "            #cut_brem_found.field(\"brem_vtx_z\").real(brem_f[itr, \"brem_vtx_z\",jentry])\n",
    "    cut_brem_lost.end_list()\n",
    "    \n",
    "    cut_brem_lost.field(\"brem_vtx_z\")\n",
    "    cut_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",
    "            cut_brem_lost.real(brem_l[itr, \"brem_vtx_z\",jentry])\n",
    "    cut_brem_lost.end_list()\n",
    "        \n",
    "    cut_brem_lost.end_record()\n",
    "\n",
    "cut_brem_lost = ak.Array(cut_brem_lost)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre>{energy: 9.36e+03,\n",
       " brem_photons_pe: [2.47e+03, 170, 224, 388, 3.23e+03, 809, 172, 224],\n",
       " brem_vtx_z: [400, 501, 638, 667, 677, 709, 8.58e+03, 9.28e+03]}\n",
       "---------------------------------------------------------------------\n",
       "type: {\n",
       "    energy: float64,\n",
       "    brem_photons_pe: var * float64,\n",
       "    brem_vtx_z: var * float64\n",
       "}</pre>"
      ],
      "text/plain": [
       "<Record {energy: 9.36e+03, ...} type='{energy: float64, brem_photons_pe: va...'>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#data in cut_brem_found and cut_brem_lost\n",
    "\n",
    "cut_length_found = ak.num(cut_brem_found[\"brem_photons_pe\"],axis=-1)\n",
    "cut_length_lost = ak.num(cut_brem_lost[\"brem_photons_pe\"], axis=-1)\n",
    "\n",
    "cut_brem_found[1]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### in magnet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "inmagnet_found = ak.ArrayBuilder()\n",
    "\n",
    "for itr in range(ak.num(cut_brem_found, axis=0)):\n",
    "    \n",
    "    inmagnet_found.begin_record()\n",
    "    inmagnet_found.field(\"energy\").real(cut_brem_found[itr,\"energy\"])\n",
    "    \n",
    "    inmagnet_found.field(\"brem_photons_pe\")\n",
    "    inmagnet_found.begin_list()\n",
    "    for jentry in range(cut_length_found[itr]):\n",
    "        if (cut_brem_found[itr, \"brem_vtx_z\", jentry]>1500):\n",
    "            if cut_brem_found[itr, \"brem_vtx_z\", jentry]<=9500:\n",
    "                inmagnet_found.real(cut_brem_found[itr,\"brem_photons_pe\",jentry])\n",
    "            else:\n",
    "                continue\n",
    "    inmagnet_found.end_list()\n",
    "    \n",
    "    inmagnet_found.field(\"brem_vtx_z\")\n",
    "    inmagnet_found.begin_list()\n",
    "    for jentry in range(cut_length_found[itr]):\n",
    "        if cut_brem_found[itr, \"brem_vtx_z\", jentry]>1500:\n",
    "            if cut_brem_found[itr,\"brem_vtx_z\",jentry]<=9500:\n",
    "                inmagnet_found.real(cut_brem_found[itr,\"brem_vtx_z\",jentry])\n",
    "            else:\n",
    "                continue\n",
    "    inmagnet_found.end_list()\n",
    "    inmagnet_found.end_record()\n",
    "    \n",
    "\n",
    "inmagnet_found = ak.Array(inmagnet_found)\n",
    "\n",
    "\n",
    "inmagnet_lost = ak.ArrayBuilder()\n",
    "\n",
    "for itr in range(ak.num(cut_brem_lost, axis=0)):\n",
    "    \n",
    "    inmagnet_lost.begin_record()\n",
    "    inmagnet_lost.field(\"energy\").real(cut_brem_lost[itr,\"energy\"])\n",
    "    \n",
    "    inmagnet_lost.field(\"brem_photons_pe\")\n",
    "    inmagnet_lost.begin_list()\n",
    "    for jentry in range(cut_length_lost[itr]):\n",
    "        if (cut_brem_lost[itr, \"brem_vtx_z\", jentry]>1500):\n",
    "            if cut_brem_lost[itr, \"brem_vtx_z\", jentry]<=9500:\n",
    "                inmagnet_lost.real(cut_brem_lost[itr,\"brem_photons_pe\",jentry])\n",
    "            else:\n",
    "                continue\n",
    "    inmagnet_lost.end_list()\n",
    "    \n",
    "    inmagnet_lost.field(\"brem_vtx_z\")\n",
    "    inmagnet_lost.begin_list()\n",
    "    for jentry in range(cut_length_lost[itr]):\n",
    "        if cut_brem_lost[itr, \"brem_vtx_z\", jentry]>1500:\n",
    "            if cut_brem_lost[itr,\"brem_vtx_z\",jentry]<=9500:\n",
    "                inmagnet_lost.real(cut_brem_lost[itr,\"brem_vtx_z\",jentry])\n",
    "            else:\n",
    "                continue\n",
    "    inmagnet_lost.end_list()\n",
    "    inmagnet_lost.end_record()\n",
    "    \n",
    "\n",
    "inmagnet_lost = ak.Array(inmagnet_lost)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "cutoff_energy=350\n",
    "#possibly: instead of checking if any photons exceed the cutoff, use the sum of all photon energies to separate nobrem and brem\n",
    "\n",
    "inmagnet_brem_found = inmagnet_found[ak.sum(inmagnet_found[\"brem_photons_pe\"],axis=-1,keepdims=False)>=cutoff_energy]\n",
    "magnet_energy_found = ak.to_numpy(inmagnet_brem_found[\"energy\"])\n",
    "magnet_eph_found = ak.to_numpy(ak.sum(inmagnet_brem_found[\"brem_photons_pe\"], axis=-1, keepdims=False))\n",
    "magnet_residual_found = magnet_energy_found - magnet_eph_found\n",
    "magnet_energyloss_found = magnet_eph_found/magnet_energy_found\n",
    "\n",
    "\n",
    "inmagnet_brem_lost = inmagnet_lost[ak.sum(inmagnet_lost[\"brem_photons_pe\"],axis=-1,keepdims=False)>=cutoff_energy]\n",
    "magnet_energy_lost = ak.to_numpy(inmagnet_brem_lost[\"energy\"])\n",
    "magnet_eph_lost = ak.to_numpy(ak.sum(inmagnet_brem_lost[\"brem_photons_pe\"], axis=-1, keepdims=False))\n",
    "magnet_residual_lost = magnet_energy_lost - magnet_eph_lost\n",
    "magnet_energyloss_lost = magnet_eph_lost/magnet_energy_lost"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "24784.620206013704"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ak.mean(magnet_eph_lost)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjQAAAHRCAYAAAB96iOvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABNJ0lEQVR4nO3dT4zb2J0v+m8FF/E1nmGzytnMarqobLIamypv7saTFJlep0tyDTABZjFdYjdq10iLXdnY3qSaGm+Nbqp6HuYB2bgo98VbJqRf2rjAYHBdZDubGbxJi+UsHuYt0iW24fs81ZNpvkXlsCmJ+kdRKlH1/QCCXRJFHh5R5E+Hv3POShRFEYiIiIgK7DvnXQAiIiKiaTGgISIiosJjQENERESFx4CGiIiICo8BDRERERUeAxoiIiIqPAY0REREVHj/5bwLMEu6rqNUKuHLL7/ErVu3UKlUzrtIM3cR95mIiGhpA5pqtQpZllGv1wEA5XIZkiRBVdVzLtnsXMR9JiIiAoCVZRwpOAgClEoltNttyLIMAGg0GnAcB47jnHPpZmNW++y6LtbW1qAoSl5FJSIiyt1S5tD4vg8A8YUdABRFgeu6CMPwnEo1W7Pa5zAMcXJyMm3xiIiIZmopA5pnz55BkqSu59bW1gBgKS7OhmEgCIKu55Z9n5eJCD4vijAM0Ww2+45Zml4QBGg2m0v7Q43Gd9HOK2kKE9CEYYhGo4GVlRWsrKxgdXUV1WoVmqahVCqh0Wh0LSsu5r0W5aTabDZRKpXi/dE0Da7rAjg7MDVNi18T+9ZsNgEA169fj5cT7ynCPl9E1Wo1/hzFwzCM1GV930e1WkWpVMLq6irK5TI0TYNhGGg2myiXyxNt23XdrmOsXC6j1Wr1LddqtVAul+Nl8jwxNptNrK+vQ9d1XnQz0HV94OfeaDRQKpWg6zp/tCyBIAig6zp0XUe1Wh36nRn3vGIYBqrVanytHHQtGHe5hRcVjKIoEYDIsqz4Odu2IwBRpVKJoiiKTNOMJEnqep/neRGAqN1uz7W8w5imGQGIZFnue01V1UhRlKjT6cTPeZ4X1Wq1qFarRaZpRvV6PX59Vvts23bkOE7m919knU4nkmU5qlQqXQ/P8/qWrdfr8TGc/Lw6nU5Uq9UiAFGWr2un04nfm/zO9LIsq+/4yYvYt7T9nka73e76fiwjSZKGfm6VSmXhzmuztoyfuzhXJ78j9Xo9kiSp77Md97yiKEp8TYyiKHIcJ5IkKfNyRVC4Xk4ickz23BFdk8WvT1mW+yJb8QsmmWNy3trtNgD0Rdait1JvMq+iKCiXy7BtG5IkYXt7O77NlNc+67re9ffR0RHW1tZg23bX85Zljb1O4OyXuud5Q5fRNG2pupnv7+/DNM2R+1StVtFqtWCaZtxDTZAkCZZloVwu930245AkCfV6HY1GA7Zto1arpS7neR729vYmXv84RIti3qrVavxdWEYi/+3OnTsDlxnUKrvMlvFz39nZgaIoXZ0vTNNEo9GAYRhd599xziuGYcD3fTx58iR+TlVVyLKMnZ2d+Fw87nKFcd4R1SREFJvWooE//QrtdDpRu93u+9VimmakKEqu5TFNc6pfRrIsd5Wz3W5HiqIM/EVWqVQix3Ei0zSjTqcT1ev1eNlZ7fNFaKHpdDpD6z0rSZIi0zSH/tIRrYvjfE5ZP8tkK82gskiSNLNfvaIlMs9ffBehZaJSqYz8zEXr3TLXQ1LRPnfbtiNFUYaWV5y7a7Va32uqqvbt7zjnFUmSUq+Tvd/FcZcrikIFNKKSez/4tItCpVKJTNOM/1YUJfcLs23bqQfhOMRFRjTz27YdybI81gFUr9dTl5vFPl+EgCaKojhAlGW5qw6zEseqeEiSFNm23becCGrHCabSToqe50WVSiVSVTWSZTmq1+up7xUXgmTTsmDbdurzw3Q6nciyrPgYE/+XJCmqVCpdwZGoi+Ryg8oibrHV6/VIVdVIVdW+4098VwBEqqr2NbePWkey7OL4HlSm5LpqtVoky/LIz8pxnPhzF+vzPC8us6IoXeUddJEe57gQAY24HS1JUiRJUuo5UtSDuL2YXGbYcSTqS1XVyLKsqN1uR6qqRpIkRaqqdt32lmU5kiRp4HGYxvO8SJKkuG6S9SBuV5qmOfBzF+URdW7bdtTpdOIyTvJ9lmU5XnetVos/m0kv7pZlRbIsR7VabeQPBXH9Sqsz8fmKc8c45xXxwz/t+yWOTXENGWe5KJr+GJj1MSQUKqARB23yA0wGM70nhGSuSd6/voW06HYclmV1fTGTH+oow6L9cfdZBGP1en3oL4hJA5px15s3ccIWX/BRF85hxJcqmaM0qU6nE9m2HQdJ4gSUvIiIX2biYj8pz/MiVVXjv8V3IS3ITm6rd5+yBL7JvB5VVaN6vR4HRr2tqOIkrKpq/MtSvD95sREXtuSFQ3xPei9K4kLXe3yNs452ux2XU5R9UJkqlUrXidWyrLEukKJ8yXrtvUgk19n7nPgsRx1/osyqqka1Wq3ruBfHRjIQSH43xQ/AUcdRu93u2o6oL7E/YtuO43TV7SQBQG8epCAufMKgzz2K+lu8VVWduIWh97sz6Ef0IFnOHaIe085TyYAuisY7r4jjPa3MySBm3OWiaPpjYB7HUBQVLKARH574MooI9TxbEMSFdFLiAxO/TPJoFRiX+OUhDEs8nCSgmWS9syC+MMkWLHFCyBLtT/IraxRxwk5e5JK/5Aet3/O8yDTNrod4f+8v/SiK4uMpbX3iB0GyLtrtduZkYHGy7/2Me394DFouedEV+5P2XRLf9bRf7r0XtnHXIT6P3u9db5nSfimO810VAWTvRSqtib83EVw8N855pfcXvCAuduJYEZ9B2vdgnONIXOB6L35pn41YdtJzmihzkmiZEYYFNL0X6yzf+d6WPHGOHvb9n7Z1N3lLuJfY32Hn597zyrDPWhyXiqKMvZww7TEwj2OoMN22RfdkRVHgeV78ANDV5XneKpUKNE2DpmkTdUsV5RXJWIZhzKVbq6ZpODw8xMHBAYCzROowDAdOj6Cq6lhTJ0y63lk4OjoCAGxvb8fJddvb2wCydV2v1Wpot9vQNA2bm5vQdT1zd8ZKpRIneU+SUK0oCmq1GgzDiJPHVVVFEATwfR/7+/txd8tqtRq/T9RFkni/6P4PnCUeTpsM3Jt0LpKXe5PaNzY2uv6WJCmuT7E/aSNSi/WNqrcs60hLLE12gZZlGY1Go2tYiN7E7TSyLENRlNRu8kEQxPsdhiHCMOyrw1ar1fV5jtK7z72fgdjPW7du9ZVlkuOot75EuZPJyeI50elhXKLMyTp79OjR2B0FxHdFJNibpjnR9oHuziY7OzsAzr4jaceJ6GZdLpfjUdrHOTZ6SZIUJ+sn69513fi7OqhjR9p5ZZxk6TAMx14urbxJkx4DszyGChPQiA8tecApihJfQCf58qcRkzpmeViWBdd1sbq6OlZg5ft+fBJTFCX+EgwanyQvzWYTruvi4OAAkiSh2WzCMAw4jjPwCzPOQZ9lvbPg+z4qlUrXyV1cOKbpDVKpVOB5XtcJP0vwqaoqFEWJy5S8wA8LlCRJivdJ/CvGirFtu+vR6XQQRVFqICl6L4iB7gDg8PBwYM+nrEQZJwn+ho19I+pp1PryWEcv0ZvGMAyUSqWJxujpvUA3m824x5K4+BweHvadu8Tyw3o3jSI+/9797f0+ZzmOZkUch6JufN/vC4JHyRLEpHFdF61WKw6SeiXHbGm321N/hyzLgmma8H0/Hg8qCIL4vDWsHnrPK+Kcm3aOEs/Jsjz2ckVSmIBGBAqapnU9L75w4pdOVpZlod1uZ3pYlgVFUeA4zlgnALEvYlnxC2DWo6mKE2wQBHHw1G63pz5pzWq9kxAXgd6uzc+ePQMwfcCbNM1xJstyfJJKXlzSWlTSiPcmWzYmIT4f0zTRarWwsbGRe/dXUcYsJ8NhvwjHDUrzWIcgyzKOj4/jVrFyudzVwjVM7wVaXLRUVY3XYVlW38Xw0aNHUFV1qs8lOZzDMFmPo1mQJAmVSgWu6yIIAjx69GjioQqCIIhb/pK
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "plt.hist(magnet_energyloss_found, alpha=0.5, bins=80, density=True, histtype='bar', color=\"blue\", label=\"found\")\n",
    "plt.hist(magnet_energyloss_lost, alpha=0.5, bins=80, density=True, histtype='bar', color=\"darkorange\", label=\"lost\")\n",
    "\n",
    "#plt.vlines(ak.mean(both_eloss),0,3,colors=\"red\", label=\"mean\")\n",
    "plt.xlabel(r\"Energyloss Ratio $E_\\gamma/E_0$\")\n",
    "plt.ylabel(\"counts (normed)\")\n",
    "plt.title(r'$B^0\\rightarrow K^{\\ast 0} e^+e^-$, $p>5$GeV, photons w/ brem_vtx_z$<9500$mm')\n",
    "plt.legend(title=\"LHCb Simulation\", title_fontsize=15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjQAAAHhCAYAAAB5kCQCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABN2UlEQVR4nO3dT4wb+Z3f/U8ngRUBwqi6tZecMl30HjaXzBRbF1/Ga5GZc9ykdMgCuUyTHvTN2CGncxnNxT3F1bVhF3sC7AP4IpEa4zk6pJ54AjxG8KhZoz0FWD8syYdgn8NOkyNoocixXc+h51cm2SSbLP6t7vcLIKRmFX/1q18Vq7781e/PRhiGoQAAABLsn606AwAAALMioAEAAIlHQAMAABKPgAYAACQeAQ0AAEg8AhoAAJB4BDQAACDx/sWqM7BIxWJRqVRKX3/9tW7fvq1cLrfqLC3cVdxnAAAubUCTz+dl27ZKpZIkKZ1Oy7IsZTKZFedsca7iPgMAIEkbl3Gk4CAIlEql1G63Zdu2JKlSqajRaKjRaKw4d4uxqH1uNpva2tqS4zjzyioAAHN3KdvQ+L4vSdGNXZIcx1Gz2VS3211RrhZrUfvc7XZ1eno6a/YAAFioSxnQPH36VJZl9b23tbUlSZfi5lwulxUEQd97l32fLxMTfF4V3W5X1Wr13DmL2QVBoGq1eml/qGFyV+26MkxiApput6tKpaKNjQ1tbGxoc3NT+Xxe2WxWqVRKlUqlb11zMx+0LhfVarWqVCoV7U82m1Wz2ZR0dmJms9lomdm3arUqSbp161a0nvlMEvb5Ksrn89FxNK9yuTx0Xd/3lc/nlUqltLm5qXQ6rWw2q3K5rGq1qnQ6PdW2m81m3zmWTqdVr9fPrVev15VOp6N15nlhrFar2t7eVrFY5KYbQ7FYHHncK5WKUqmUisUiP1ougSAIVCwWVSwWlc/nx35nJr2ulMtl5fP56F456l4w6XprL0wYx3FCSaHnedF7tVotlBTmcrkwDMPQdd3Qsqy+z7VarVBS2G63l5rfcVzXDSWFtm2fW5bJZELHccJOpxO912q1wkKhEBYKhdB13bBUKkXLF7XPtVotbDQasT9/lXU6ndC27TCXy/W9Wq3WuXVLpVJ0Dvcer06nExYKhVBSGOfr2ul0os/2fmcGeZ537vyZF7Nvw/Z7Fu12u+/7cRlZljX2uOVyubW7ri3aZTzu5lrd+x0plUqhZVnnju2k1xXHcaJ7YhiGYaPRCC3Lir1eEiSul5OJHHt77piuyebXp23b5yJb8wumt43JqrXbbUk6F1mb3kqDjXkdx1E6nVatVpNlWbp37170mGle+1wsFvv+Pjk50dbWlmq1Wt/7nudNnKZ09ku91WqNXSebzV6qbuaHh4dyXffCfcrn86rX63JdN+qhZliWJc/zlE6nzx2bSViWpVKppEqlolqtpkKhMHS9Vqulg4ODqdOfhKlRnLd8Ph99Fy4j0/7t7t27I9cZVSt7mV3G4763tyfHcfo6X7iuq0qlonK53Hf9neS6Ui6X5fu+njx5Er2XyWRk27b29vaia/Gk6yXGqiOqaZgodliNhr79FdrpdMJ2u33uV4vruqHjOHPNj+u6M/0ysm27L5/tdjt0HGfkL7JcLhc2Go3Qdd2w0+mEpVIpWndR+3wVamg6nc7Yco/LsqzQdd2xv3RM7eIkxynuseytpRmVF8uyFvar19REzvMX31Womcjlchcec1N7d5nLoVfSjnutVgsdxxmbX3PtLhQK55ZlMplz+zvJdcWyrKH3ycHv4qTrJUWiAhpTyIMHfthNIZfLha7rRn87jjP3G3OtVht6Ek7C3GRMNX+tVgtt257oBCqVSkPXW8Q+X4WAJgzDKEC0bbuvDOMy56p5WZYV1mq1c+uZoHaSYGrYRbHVaoW5XC7MZDKhbdthqVQa+llzI+itWjZqtdrQ98fpdDqh53nROWb+b1lWmMvl+oIjUxa9643Ki3nEViqVwkwmE2YymXPnn/muSAozmcy56vaL0ujNuzm/R+WpN61CoRDatn3hsWo0GtFxN+m1Wq0oz47j9OV31E16kvPCBDTmcbRlWaFlWUOvkaYczOPF3nXGnUemvDKZTOh5Xthut8NMJhNalhVmMpm+x962bYeWZY08D4dptVqhZVlR2fSWg3lc6bruyONu8mPKvFarhZ1OJ8rjNN9n27ajtAuFQnRspr25e54X2rYdFgqFC38omPvXsDIzx9dcOya5rpgf/sO+X+bcNPeQSdYLw9nPgUWfQ0aiAhpz0vYewN5gZvCC0NvWZN6/vo1h0e0kPM/r+2L2HtSLjIv2J91nE4yVSqWxvyCmDWgmTXfezAXbfMEvunGOY75UvW2UptXpdMJarRYFSeYC1HsTMb/MzM1+Wq1WK8xkMtHf5rswLMju3dbgPsUJfHvb9WQymbBUKkWB0WAtqrkIZzKZ6Jel+Xzvzcbc2HpvHOZ7MnhTMje6wfNrkjTa7XaUT5P3UXnK5XJ9F1bP8ya6QZr89Zbr4E2iN83B98yxvOj8M3nOZDJhoVDoO+/NudEbCPR+N80PwIvOo3a73bcdU15mf8y2G41GX9lOEwAMtoM0zI3PGHXcw/B8jXcmk5m6hmHwuzPqR/Qoca4dphyHXad6A7ownOy6Ys73YXnuDWImXS8MZz8HlnEOhWHCAhpz8MyX0USoq6xBMDfSaZkDZn6ZzKNWYFLml4cxruHhNAHNNOkugvnC9NZgmQtCnGh/ml9ZFzEX7N6bXO8v+VHpt1qt0HXdvpf5/OAv/TAMo/NpWHrmB0FvWbTb7diNgc3FfvAYD/7wGLVe703X7M+w75L5rg/75T54Y5s0DXM8Br93g3ka9ktxku+qCSAHb1LDqvgHG4Kb9ya5rgz+gjfMzc6cK+YYDPseTHIemRvc4M1v2LEx6057TTN57mVqZoxxAc3gzTrOd36wJs9co8d9/2et3e19JDzI7O+46/PgdWXcsTbnpeM4E69nzHoOLOMcSky3bdM92XEctVqt6CWpr8vzsuVyOWWzWWWz2am6pZr8msZY5XJ5Kd1as9msHj16pOPjY0lnDam73e7I6REymcxEUydMm+4inJycSJLu3bsXNa67d++epHhd1wuFgtrttrLZrO7cuaNisRi7O2Mul4saeU/ToNpxHBUKBZXL5ajxeCaTURAE8n1fh4eHUXfLfD4ffc6URS/zedP9XzpreDhrY+DBRuem8fJgo/adnZ2+vy3LisrT7M+wEalNeheVW5w0hjUs7e0Cbdu2KpVK37AQgw23h7FtW47jDO0mHwRBtN/dblfdbvdcGdbr9b7jeZHBfR48BmY/b9++fS4v05xHg+Vl8t3bONm8Zzo9TMrkubfMHj58OHFHAfNdMQ3sXdedavtSf2eTvb09SWffkWHnielmnU6no1HaJzk3BlmWFTXW7y37ZrMZfVdHdewYdl2ZpLF0t9udeL1h+e017TmwyHMoMQGNOWi9J5zjONENdJov/zBmUsc4L8/z1Gw2tbm5OVFg5ft+dBFzHCf6Eowan2ReqtWqms2mjo+PZVmWqtWqyuWyGo3GyC/MJCd9nHQXwfd95XK5vou7uXHM0hskl8up1Wr1XfDjBJ+ZTEaO40R56r3BjwuULMuK9sn8a8aKqdVqfa9Op6MwDIcGkqb3ghnoTpIePXo0sudTXCaP0wR/48a+MeV0UXrzSGOQ6U1TLpeVSqWmGqNn8AZdrVajHkvm5vPo0aNz1y6z/rjeTRcxx39wfwe/z3HOo0Ux56EpG9/3zwXBF4kTxAzTbDZVr9ejIGlQ75gt7XZ75u+Q53lyXVe+70fjQQVBEF23xpXD4HXFXHOHXaPMe7ZtT7xekiQmoDGBQjab7XvffOHML524PM9Tu92O9fI8T47jqNFoTHQBMPti1jW/ABY9mqq5wAZBEAVP7XZ75ovWotKdhrkJDHZtfvr0qaTZA95es5xntm1HF6nem8uwGpVhzGd7azamYY6P67qq1+va2dmZe/dXk8c4F8NxvwgnDUrnkYZh27aeP38e1Yql0+m+Gq5xBm/Q5qaVyWSiNDzPO3c
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nstart = 0\n",
    "nend = 5e4\n",
    "plt.hist(magnet_residual_found, alpha=0.5, bins=70, density=True, histtype='bar', color=\"blue\", label=\"found\", range=[nstart, nend])\n",
    "plt.hist(magnet_residual_lost, alpha=0.5, bins=70, density=True, histtype='bar', color=\"darkorange\", label=\"lost\", range=[nstart, nend])\n",
    "\n",
    "#plt.vlines(ak.mean(both_eloss),0,3,colors=\"red\", label=\"mean\")\n",
    "#plt.xlim(0,50000)\n",
    "plt.xlabel(r\"Residual Energy in magnet $E_\\gamma$\")\n",
    "plt.ylabel(\"counts (normed)\")\n",
    "plt.title(r'$B^0\\rightarrow K^{\\ast 0} e^+e^-$, $p>5$GeV, photons w/ brem_vtx_z$<9500$mm')\n",
    "plt.legend(title=\"LHCb Simulation\", title_fontsize=15)\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": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cut_length_found[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Split in Upstream and Downstream Events and analyse separately"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#try to find a split between energy lost before and after the magnet (z~5000mm)\n",
    "\n",
    "upstream_found = ak.ArrayBuilder()\n",
    "downstream_found = ak.ArrayBuilder()\n",
    "\n",
    "for itr in range(ak.num(cut_brem_found, axis=0)):\n",
    "    upstream_found.begin_record()\n",
    "    upstream_found.field(\"energy\").real(cut_brem_found[itr,\"energy\"])\n",
    "    \n",
    "    downstream_found.begin_record()\n",
    "    downstream_found.field(\"energy\").real(cut_brem_found[itr,\"energy\"])\n",
    "    \n",
    "    upstream_found.field(\"brem_photons_pe\")\n",
    "    downstream_found.field(\"brem_photons_pe\")\n",
    "    upstream_found.begin_list()\n",
    "    downstream_found.begin_list()\n",
    "    for jentry in range(cut_length_found[itr]):\n",
    "        if (cut_brem_found[itr, \"brem_vtx_z\", jentry]>5000):\n",
    "            if cut_brem_found[itr, \"brem_vtx_z\", jentry]<=9500:\n",
    "                downstream_found.real(cut_brem_found[itr,\"brem_photons_pe\",jentry])\n",
    "            else:\n",
    "                continue\n",
    "        else:\n",
    "            upstream_found.real(cut_brem_found[itr,\"brem_photons_pe\", jentry])            \n",
    "    upstream_found.end_list()\n",
    "    downstream_found.end_list()\n",
    "    \n",
    "    upstream_found.field(\"brem_vtx_z\")\n",
    "    downstream_found.field(\"brem_vtx_z\")\n",
    "    upstream_found.begin_list()\n",
    "    downstream_found.begin_list()\n",
    "    for jentry in range(cut_length_found[itr]):\n",
    "        if cut_brem_found[itr, \"brem_vtx_z\", jentry]>5000:\n",
    "            if cut_brem_found[itr,\"brem_vtx_z\",jentry]<=9500:\n",
    "                downstream_found.real(cut_brem_found[itr,\"brem_vtx_z\",jentry])\n",
    "            else:\n",
    "                continue\n",
    "        else:\n",
    "            upstream_found.real(cut_brem_found[itr, \"brem_vtx_z\",jentry])\n",
    "    upstream_found.end_list()\n",
    "    downstream_found.end_list()\n",
    "    upstream_found.end_record()\n",
    "    downstream_found.end_record()\n",
    "    \n",
    "\n",
    "upstream_found = ak.Array(upstream_found)\n",
    "downstream_found = ak.Array(downstream_found)\n",
    "\n",
    "\n",
    "upstream_lost = ak.ArrayBuilder()\n",
    "downstream_lost = ak.ArrayBuilder()\n",
    "\n",
    "for itr in range(ak.num(cut_brem_lost, axis=0)):\n",
    "    upstream_lost.begin_record()\n",
    "    upstream_lost.field(\"energy\").real(cut_brem_lost[itr,\"energy\"])\n",
    "    \n",
    "    downstream_lost.begin_record()\n",
    "    downstream_lost.field(\"energy\").real(cut_brem_lost[itr,\"energy\"])\n",
    "    \n",
    "    upstream_lost.field(\"brem_photons_pe\")\n",
    "    downstream_lost.field(\"brem_photons_pe\")\n",
    "    upstream_lost.begin_list()\n",
    "    downstream_lost.begin_list()\n",
    "    for jentry in range(cut_length_lost[itr]):\n",
    "        if (cut_brem_lost[itr, \"brem_vtx_z\", jentry]>5000):\n",
    "            if cut_brem_lost[itr, \"brem_vtx_z\", jentry]<=9500:\n",
    "                downstream_lost.real(cut_brem_lost[itr,\"brem_photons_pe\",jentry])\n",
    "            else:\n",
    "                continue\n",
    "        else:\n",
    "            upstream_lost.real(cut_brem_lost[itr,\"brem_photons_pe\", jentry])            \n",
    "    upstream_lost.end_list()\n",
    "    downstream_lost.end_list()\n",
    "    \n",
    "    upstream_lost.field(\"brem_vtx_z\")\n",
    "    downstream_lost.field(\"brem_vtx_z\")\n",
    "    upstream_lost.begin_list()\n",
    "    downstream_lost.begin_list()\n",
    "    for jentry in range(cut_length_lost[itr]):\n",
    "        if cut_brem_lost[itr, \"brem_vtx_z\", jentry]>5000:\n",
    "            if cut_brem_lost[itr,\"brem_vtx_z\",jentry]<=9500:\n",
    "                downstream_lost.real(cut_brem_lost[itr,\"brem_vtx_z\",jentry])\n",
    "            else:\n",
    "                continue\n",
    "        else:\n",
    "            upstream_lost.real(cut_brem_lost[itr, \"brem_vtx_z\",jentry])\n",
    "    upstream_lost.end_list()\n",
    "    downstream_lost.end_list()\n",
    "    upstream_lost.end_record()\n",
    "    downstream_lost.end_record()\n",
    "    \n",
    "\n",
    "upstream_lost = ak.Array(upstream_lost)\n",
    "downstream_lost = ak.Array(downstream_lost)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre>{energy: 4.62e+04,\n",
       " brem_photons_pe: [3.26e+03, 4.45e+03, 178, 1.45e+04, 1.1e+03, 3.79e+03],\n",
       " brem_vtx_z: [162, 187, 387, 487, 1.34e+03, 2.32e+03]}\n",
       "-------------------------------------------------------------------------\n",
       "type: {\n",
       "    energy: float64,\n",
       "    brem_photons_pe: var * float64,\n",
       "    brem_vtx_z: var * float64\n",
       "}</pre>"
      ],
      "text/plain": [
       "<Record {energy: 4.62e+04, ...} type='{energy: float64, brem_photons_pe: va...'>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "upstream_found[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "upstream: cutoff energy = 350MeV, sample size: 1562\n",
      "eff =  0.9181 +/- 0.007\n"
     ]
    }
   ],
   "source": [
    "#plot efficiency against cutoff energy \n",
    "up_efficiencies = []\n",
    "up_deff = []\n",
    "\n",
    "\n",
    "for cutoff_energy in range(0,10050,200):\n",
    "\tup_nobrem_f = upstream_found[ak.sum(upstream_found[\"brem_photons_pe\"],axis=-1,keepdims=False)<cutoff_energy]\n",
    "\tup_nobrem_l = upstream_lost[ak.sum(upstream_lost[\"brem_photons_pe\"],axis=-1,keepdims=False)<cutoff_energy]\n",
    "\t\n",
    "\tif ak.num(up_nobrem_f,axis=0)+ak.num(up_nobrem_l,axis=0)==0:\n",
    "\t\tup_efficiencies.append(0)\n",
    "\t\tup_deff.append(0)\n",
    "\t\tcontinue\n",
    "\n",
    "\teff = t_eff(up_nobrem_f,up_nobrem_l)\n",
    "\tdeff = eff_err(up_nobrem_f,up_nobrem_l)\n",
    "\tup_efficiencies.append(eff)\n",
    "\tup_deff.append(deff)\n",
    "\n",
    "\t#print(\"\\ncutoff = \",str(cutoff_energy),\"MeV, sample size: \",ak.num(up_nobrem_f,axis=0)+ak.num(up_nobrem_l,axis=0))\n",
    "\t#print(\"eff = \",np.round(eff,4), \"+/-\", np.round(eff_err(up_nobrem_f, up_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",
    "up_nobrem_found = upstream_found[ak.sum(upstream_found[\"brem_photons_pe\"],axis=-1,keepdims=False)<cutoff_energy]\n",
    "up_nobrem_lost = upstream_lost[ak.sum(upstream_lost[\"brem_photons_pe\"],axis=-1,keepdims=False)<cutoff_energy]\n",
    "\n",
    "print(\"\\nupstream: cutoff energy = 350MeV, sample size:\",ak.num(up_nobrem_found,axis=0)+ak.num(up_nobrem_lost,axis=0))\n",
    "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nobrem_vertices\n",
      "upstream: cutoff energy = 350MeV, sample size: 1562\n",
      "eff =  0.9181 +/- 0.007\n",
      "\n",
      "downstream: cutoff energy = 350MeV, sample size: 5131\n",
      "eff =  0.8864 +/- 0.004\n"
     ]
    }
   ],
   "source": [
    "down_efficiencies = []\n",
    "down_deff = []\n",
    "\n",
    "for cutoff_energy in range(0,10050,200):\n",
    "\tdown_nobrem_f = downstream_found[ak.sum(downstream_found[\"brem_photons_pe\"],axis=-1,keepdims=False)<cutoff_energy]\n",
    "\tdown_nobrem_l = downstream_lost[ak.sum(downstream_lost[\"brem_photons_pe\"],axis=-1,keepdims=False)<cutoff_energy]\n",
    "\n",
    "\tif ak.num(down_nobrem_f,axis=0)+ak.num(down_nobrem_l,axis=0)==0:\n",
    "\t\tdown_efficiencies.append(0)\n",
    "\t\tdown_deff.append(0)\n",
    "\t\tcontinue\n",
    "\teff = t_eff(down_nobrem_f,down_nobrem_l)\n",
    "\tdeff = eff_err(down_nobrem_f,down_nobrem_l)\n",
    "\tdown_efficiencies.append(eff)\n",
    "\tdown_deff.append(deff)\n",
    "\n",
    "\n",
    "\t#print(\"\\ncutoff = \",str(cutoff_energy),\"MeV, sample size: \",ak.num(down_nobrem_f,axis=0)+ak.num(down_nobrem_l,axis=0))\n",
    "\t#print(\"eff = \",np.round(eff,4), \"+/-\", np.round(eff_err(down_nobrem_f, down_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",
    "down_nobrem_found = downstream_found[ak.sum(downstream_found[\"brem_photons_pe\"],axis=-1,keepdims=False)<cutoff_energy]\n",
    "down_nobrem_lost = downstream_lost[ak.sum(downstream_lost[\"brem_photons_pe\"],axis=-1,keepdims=False)<cutoff_energy]\n",
    "\n",
    "\n",
    "print(\"nobrem_vertices\\nupstream: cutoff energy = 350MeV, sample size:\",ak.num(up_nobrem_found,axis=0)+ak.num(up_nobrem_lost,axis=0))\n",
    "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",
    "\n",
    "print(\"\\ndownstream: cutoff energy = 350MeV, sample size:\",ak.num(down_nobrem_found,axis=0)+ak.num(down_nobrem_lost,axis=0))\n",
    "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABc4AAAJLCAYAAADeqgXEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB9vElEQVR4nO3dQYhjaX4g+H+YorMTV1cro8Y1DA1mSgENY3xoS2lmG0wVTEnMwT4YUqo67GH34JJYgrn0IUSCh+7GsIHiaoJFkY3Zgy+VUt2mD4vULN202V0qpekFY4MhlHsxg2u2I9TlWrKqD/32kPNUUoYiQi9CipAUvx8EmaH3ve99ep9C+c9/fPp/W0mSJAEAAAAAAERExG/d9gAAAAAAAGCVSJwDAAAAAMAEiXMAAAAAAJggcQ4AAAAAABMkzgEAAAAAYILEOQAAAAAATJA4B4ANNRwOYzQa3fYwAAAAYO1InAPAhup0OtHr9W57GAAAALB2JM4BAAAAAGCCxDkAAAAAAEyQOAcAYGFGo1EMh8PbHgYAAMC1SJwDALAwT58+jWazedvD4BKdTue2h3AnjUaj6HQ6Ua1W48GDBzEYDG57SCtvOBzG0dFRlMvl2x4KAHDHSJwDAJkcHByMkz4PHjyIer0eo9Hotoe1cIPBIOr1ejx48CC2tramvtLnXiwWo9ForPXzbzQaUSwWzzzHarU63lx2NBpFo9E4cy/q9XpEvLxX6Srz7e3tyOVy4/NsUHt7Dg4Ozsxr+vXRRx9den6v14t6vR47OzvjuZ983fd6vRgOh1EsFq+VAO50OjNfg+VyOQ4ODmaeMxgMolqtnnnNrvqnHY6OjqLVakWn01nr942bMBgMolwux87OTtTrde8lAMCN20qSJLntQQAAV3deQuHk5CQiXiYyZ52zt7eX6Tqj0Sjee++9KJVK0Ww2YzQaRbFYjOFwGM1mM3N/6yJdHRoRUSgUot/vj481Go1xYq/f70ehULiVMS7CgwcPxom8brcbpVLpTJtyuRy9Xi9KpVK02+1xgrzT6cSHH34YzWYzHj58GM+ePYt8Ph/VajUeP368sa+NVTc5p6+66PU6GAziww8/HCfD9/b2olwuRz6fj9FoFM+ePYtWqzWVLF/E6z99fUWc/Vk7z8HBQTQajbV6DxoOh7GzsxMR6/++cROOjo7Gv6S7rf+6pj9H6XseAHA3vHbbAwAArqfVas18/ODgIPL5fFQqlYVcp9FoxGAwiHa7HREvEwjHx8fjROqmmlzBmiZvUs1mM3q9XgwGg3jvvffi9PT0poe3MGliqFAozJzPdIXxrARlpVKJSqUyTmKenJzEw4cP4/nz5xJNt+To6Ci2t7fPTT7n8/mZj0/+oqhUKkWr1TrTtlAoRK1Wm0popr+ou452ux0PHjyIiJh7Bfsvf/nLyOVya5M0j5j9y0zOd95r9SZVq9VoNpt+yQEAd4xSLQDAXI6OjiLibNJnk5PmERGffPLJ+O/vv//+meMPHz6MiJeJ53UtvTCZpPzggw/OHK/X63F0dBT9fv/CBGWaJJ8s18JL9Xp9akX1sjWbzWg2m5HP52d+zZKWP4l4+XPd7XYvTFrWarXxL9IW8drP5XJRq9XG359XpmXS0dHR1DmwaOkvDQGAu0fiHAC41KrXDV6mNGGSz+dnJoMn780iVt3Oo16vL3ROJpNCk78ISetXD4fDeP78+bmrLXu93rgsSLPZjGq1GuVyOd5+++25kp93QavVinK5HNVqNYrF4lI35+x0OjEcDuPk5CRT7fE0aR4R44T4ZSqVykJ/eTa5sez+/v6FbXu9XoxGo3j8+PHCrg+TOp2O9zAAuMMkzgGAcx0dHY03A0y9/fbb480xJxMKvV4vqtXqeGV6vV4fb/D3qoODgyiXy1EsFmNnZyfK5fL4vFnS8hGNRiMivto0Lt2scHIcw+FwalPPcrl85dWwg8FgfO55JW+ePXsWES9Xy95USYHhcLjQxOvkRpFpcjzdsPGDDz6Ibrd76Qryn/zkJ7G3txfD4TCOj4+jVCrF8+fPV6LMwqrY29uL09PTqNfr0Wg0Ymdn58LX/VWlCed6vT7edLNarV6YRD86Ohr/MqZSqWT6xEC6sn2W0Wg03mB0a2srisXihYnIXC43/lkbjUYXvs6bzWbmsZ43xqOjo6n3oaOjo6kxX/YLiKzvaZOGw+HUhssXbXI67/tslvt+G++vWV8Xy+gvnbN089tyuTw1z+neDan33ntvPEeTlvVv33Vel6PRKKrV6vhrZ2dnvMkqAJBBAgCsvHa7nVQqlaRSqST5fD6p1WqXntNsNpN2u72Q6x8fHycRkUREcnx8fOY6+Xx+fLzVaiWlUinJ5XJnzun3+0k+n08qlcr4/NPT06RSqSQRkZRKpeT09HR8rNvtJqVSadxPrVZL9vb2kkKhkOzt7SW1Wm3quu12e3x/9vb2xuMqFApXet7NZnPcf7fbPXM8Hfd5x5clfZ6Lkj6HSqUyno9cLpf0+/3MfbVarblen5Oyvr6v8vOQRb/fT/b29sb9p6+fdrudFAqFJCKSfD6fNJvNa10nncdcLnftvlKnp6dJqVSa+pmc/Nrb25t53uTP2aLeN7rdblIoFMY/G8fHx+PrXPQzOfl+c97rPG1z3Z+7dIzp9ZrNZlKpVJJCoZDUarWp+zj53pTK+p6WHp+8Xvo8J98zX31uWd5n573vt/X+mvV10e12x+NYRH/Hx8fj55HOzeRrbvL1PzlXr74fLvvfvuu8Ll99fzo9PT1zfQDgchLnALDC0v9sT/6nem9vbyqpUqvVkkKhcOYrl8sluVxu5rGsSbqLEudJMp1cyOfz47E1m82ppGY6plnSREepVDr32KykXprcyeVyZ5ICl437MpOJi8k+W63W+FilUrlS39eVz+cXkuBst9tTifNcLpfk8/mZyZh5nJ6ezp1wn+f1fZ3215Ump/b29pK9vb2k2Wwm/X4/abVaUwmt60oT8mkC/ar3/lWnp6dTyf6LkueTxxfxek5/9l59LUy+V5yXxE+S6Z/5WXObJm4XJU125nK5qeudnp6Oj82a66u8p736fjl5vyeTxK/2O8/77FXu+02+v15lfBclzq/SX6FQmJlATu/D5OvqosT5q8eX8W/fVV6X6fvTq3PS7XYlzgEgI4lzAFhRabLg1f9Mp//Jvspq4KuaJ0GS/if+vNW/aYLzvGRZv98/N3kzuSrvVZNJlVkJx3RcV0kyTyYTX10Nms/nb3QOXpU+7+uOYXJV6WWrGBcp6+v7Nn4e0vtRqVTO9J++JheZvJ1cAby3t7fQOeh2u1Ov4Vefz+Rr+7Lrpivya7Xama80uZf+EmaWNJF/3vF0vOl4zksoLmqV/uSYZr0/pXPy6nvbVd/TJpOts96XLvrFzGXvs1e57zf5/nqV8V2UOM/aXzpns94v+v1+ksvlpu7tZYnzJFnuv33XeV3O+vm46JdVAMBZapwDwApKa8zmcrmpTfp6vV70er0olUrnbtR42ybroU9KaxX/4R/+4czjhUJhXCd5suZ2RMT29nZExMxaxumx86THs26mOblhZqvVitPT00iSJE5PT6NWq403zrytjeNKpVI0m8147733rrVR6OTmp5OvtcnavouW9fV9Gz8Po9FofF/L5fKZ/tPX6iI3hC2VStHtdqPf78dwOIwHDx4sbCPYUqkU/X5//P3k6ztrjepCoRDNZnNco/3o6CiePn0ajUZjvFFop9OJ0Wg0rq08+ZXWZR6NRudee3JOe73emdrTo9EoarVapnFfVTrXr87Ddd7TXu17Uq1WG7/XffLJJzPPu+h9Nut9v8n31+u+Lq7bX1orfNb7RaFQiNPT02i1WnNd+1XL+LfvIue9Lnd2diIiotFojOvWpyY33wUALidxDgArqFqtRkTE48ePx8mMdAO3NLm2buZJrKSJgMs24sviqhsHTiZo02Rg2l+r1RpvYNhoNKaSkJdJkziL+KrValGr1WJnZ+dKCfzhcDielzRRmT6vTqez0HmYlPX1fRs/D+mc5nK5mQna9L4tY/PTQqEQ7XY7jo+P4+TkJHZ2dsb34Dr
      "text/plain": [
       "<Figure size 1800x600 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#plot efficiencies wrt cutoff energy\n",
    "fig, ax = plt.subplots(nrows=1,ncols=2,figsize=(18,6))\n",
    "x_ = np.arange(0,10050,step=200)\n",
    "\n",
    "ax[0].errorbar(x_,up_efficiencies,yerr=up_deff, ls=\"\", capsize=1,fmt=\".\")\n",
    "ax[0].set(xlabel=\"cutoff energy [MeV]\",ylabel=r\"$\\epsilon$\",title=\"upstream\", ylim=[0.8,1.0])\n",
    "#ax[0].set_yticks(np.arange(0.8,1.01,step=0.02),minor=False)\n",
    "#ax[0].set_xticks(np.arange(0,10100,step=200),minor=True)\n",
    "#ax[0].grid()\n",
    "\n",
    "ax[1].errorbar(x_,down_efficiencies,yerr=down_deff, ls=\"\", capsize=1,fmt=\".\")\n",
    "ax[1].set(xlabel=\"cutoff energy [MeV]\",ylabel=r\"$\\epsilon$\",title=\"downstream\", ylim=[0.8,1.0])\n",
    "#ax[1].set_yticks(np.arange(0.8,1.01,step=0.02),minor=False)\n",
    "#ax[1].set_xticks(np.arange(0,10100,step=200),minor=True)\n",
    "#ax[1].grid(True)\n",
    "\n",
    "fig.suptitle(r\"$e^\\pm$ from $B\\rightarrow K^\\ast ee$, $p>5$GeV, nobrem electrons\")\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "brem vertices\n",
      "upstream eff =  0.851 +/- 0.004\n",
      "downstream eff =  0.836 +/- 0.005\n"
     ]
    }
   ],
   "source": [
    "cutoff_energy=350\n",
    "#possibly: instead of checking if any photons exceed the cutoff, use the sum of all photon energies to separate nobrem and brem\n",
    "\n",
    "upstream_brem_found = upstream_found[ak.sum(upstream_found[\"brem_photons_pe\"],axis=-1,keepdims=False)>=cutoff_energy]\n",
    "up_energy_found = ak.to_numpy(upstream_brem_found[\"energy\"])\n",
    "up_eph_found = ak.to_numpy(ak.sum(upstream_brem_found[\"brem_photons_pe\"], axis=-1, keepdims=False))\n",
    "up_residual_found = up_energy_found - up_eph_found\n",
    "up_energyloss_found = up_eph_found/up_energy_found\n",
    "\n",
    "\n",
    "upstream_brem_lost = upstream_lost[ak.sum(upstream_lost[\"brem_photons_pe\"],axis=-1,keepdims=False)>=cutoff_energy]\n",
    "up_energy_lost = ak.to_numpy(upstream_brem_lost[\"energy\"])\n",
    "up_eph_lost = ak.to_numpy(ak.sum(upstream_brem_lost[\"brem_photons_pe\"], axis=-1, keepdims=False))\n",
    "up_residual_lost = up_energy_lost - up_eph_lost\n",
    "up_energyloss_lost = up_eph_lost/up_energy_lost\n",
    "\n",
    "\n",
    "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",
    "\n",
    "\n",
    "downstream_brem_found = downstream_found[ak.sum(downstream_found[\"brem_photons_pe\"],axis=-1,keepdims=False)>=cutoff_energy]\n",
    "down_energy_found = ak.to_numpy(downstream_brem_found[\"energy\"])\n",
    "down_eph_found = ak.to_numpy(ak.sum(downstream_brem_found[\"brem_photons_pe\"], axis=-1, keepdims=False))\n",
    "down_residual_found = down_energy_found - down_eph_found\n",
    "down_energyloss_found = down_eph_found/down_energy_found\n",
    "\n",
    "\n",
    "downstream_brem_lost = downstream_lost[ak.sum(downstream_lost[\"brem_photons_pe\"],axis=-1,keepdims=False)>=cutoff_energy]\n",
    "down_energy_lost = ak.to_numpy(downstream_brem_lost[\"energy\"])\n",
    "down_eph_lost = ak.to_numpy(ak.sum(downstream_brem_lost[\"brem_photons_pe\"], axis=-1, keepdims=False))\n",
    "down_residual_lost = down_energy_lost - down_eph_lost\n",
    "down_energyloss_lost = down_eph_lost/down_energy_lost\n",
    "\n",
    "\n",
    "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))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "upstream:\n",
      "mean energyloss relative to initial energy (found):  0.33078325542598164\n",
      "mean energyloss relative to initial energy (lost):  0.5708618852236069\n",
      "downstream:\n",
      "mean energyloss relative to initial energy (found):  0.19104090843883118\n",
      "mean energyloss relative to initial energy (lost):  0.3051594568487781\n"
     ]
    }
   ],
   "source": [
    "print(\"upstream:\\nmean energyloss relative to initial energy (found): \",ak.mean(up_energyloss_found))\n",
    "print(\"mean energyloss relative to initial energy (lost): \", ak.mean(up_energyloss_lost))\n",
    "\n",
    "print(\"downstream:\\nmean energyloss relative to initial energy (found): \",ak.mean(down_energyloss_found))\n",
    "print(\"mean energyloss relative to initial energy (lost): \", ak.mean(down_energyloss_lost))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABbsAAAJPCAYAAABVWwkOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB/t0lEQVR4nOz9T2wjaX7nf35U0+sc4ZfIilT5MocGOoOewxgLTFZQtcAUFsiyk5w+u0VmDjAG5tAt0g3eCt1kazBAqS6tCrmuhB3M6oV34UuKVC8W2IubLIwKC9Thl2J0+rKzO92MLCMX8AI/lxSdyIVGPXZzD+mIIqUgRQb/Bcn3CyCqFH+e+JJBRj7x5ZfPs9HtdrsCAAAAAAAAAGCJvbXoAAAAAAAAAAAAmBTJbgAAAAAAAADA0iPZDQAAAAAAAABYeiS7AQAAAAAAAABLj2Q3AAAAAAAAAGDpkewGAAAAAAAAACw9kt0AAAAAAAAAgKVHshsAAAAAAAAAsPRIdgMAAAAAAAAAlh7JbgAAVpDv+/I8b9FhAAAAAAAwNyS7AQBYQUdHR7Jte9Fh4AaNRmPRISy1VqulYrGou3fv8louMd/31Wg0lM/ndffuXbmuu+iQgLFw/QEAIDlIdgMAEs113TCZtbGx0fe4e/eu7t69q3Q6rUqlIt/3Fx1ubJVKRel0+tpzzOfzarVakt4khCqVyrXXolgsSnrzWgXV3FtbWzIMI9wvaAPzd3h4eO28Bo+nT5/euH+Q0E2lUuG5733ft1oteZ6ndDo9UZKw0WhEvgez2awODw8j93FdV/l8/tp7dh6/Kmg0GnIcR7Vabak/+6vO9/3wfTRIrVaT4zhqNBqcS8xcpVJRNptVOp3W3bt3VSwWR3rfxb2WNxoNZbPZ8JjZbHakf5Pj7gcAwNrrAgCwBOr1eldSV1LXsqy+deVyOVzXbrcXFOF0GIYRPpdmsxm5TSaT6UrqZjKZ7vn5ebi8Xq93DcPoOo7TbbfbXcdxus1ms2sYRte27Tk9A1zVe06vPoa9X9vtdteyrHDbcrncbTab3U6nE57f3vXTev8H76+oz9ogtm13JS3kfRbEWq/X53bM8/Pzvs8eBgveG47jDN2u0+mszHV80Xh/Rmu3213TNLuFQiFcdn5+3rUsq2sYxo3vuzjX8lwud2190J8pl8sDjxV3PwAA0O2S7AYALIUgYTIoaRIk/QzDWEB003NTkjFI7A9LKtq23c1kMl3LsrqFQoGkxwI5jtM1TbPb6XQiH4P0frmTyWSGbus4zo1fkIzj/Py8L4kzinK5vLDP3iKS3ZlMhoTsiIJr803Xod73Ha/tZHh/RjMMI/I6dX5+HiayB11r41zLC4XCwD5LsC7quhV3PwAA8AbDmAAAlsKzZ8/C/3/06NG19dvb25Le/GR+WX8G3zsExePHj6+tLxaLqtVqarfbKpfLA9sJhi/pHcoEb16/ef8M3LZt2bYt0zQjH1GCoUEkKZPJqNlsDtxWkgqFgur1uiRN5b1vGIYKhUL496AhTHrVarW+fVZZMHQMbuZ5nlzXVSaT4Vo0J+vw/uwd0mvUa14w1Fkmk7m2zjCMsF9RqVQi9x/3Wu55nmq1miRFXhuD4ceuHi/ufgAA4BskuwEASyG4eTdNMzJp0jtG8NnZ2VxiKhaLUx2buDdB0XtDHozH7HmeXrx4IcuyBu4f3Pzbtq18Pq9sNqt79+6NlLBcdY7jKJvNKp/PK51Oz3xCsUajIc/zdHZ2NtZY2kGiW1KYxL5JLpeLTOLE1Tu56cHBwdBtW62WfN/X3t7e1I6fVI1Gg8/SGBzHkfRNgg6ztervT8/zwjk8PM/T559/PvKXKMH1/r333otcH1x3g+v21X3HvZYH19BB12XLsmQYhjzP6/u3P+5+AADgGyS7AQCJ57puWL2Vy+Uitzk9PZX0pkJrWBXsNHmeN9WEae8EV0FCO5g08PHjx2o2mzfe2H/++ecql8vyPE+dTkeZTEYvXryY22uSdOVyWefn5yoWi6pUKkqlUmEV3bQFSeJisRhO/JjP54cmS2q1WphoyeVyY1XDBlWHUXzfDye53NjYUDqdHpoUMwwj/Kz5vj/0fW7b9tixRgkmYwvOR+9krONU5PfuN8qknYeHh+EEcKlUqi+Gq/Ht7u6Gfz98+DCcJHeSdqU3r3GtVuvbplar9Z2vQc/D933l8/nwkUqllEqlRkowu67bN+FtKpUKX+fDw8Nr63rj931f2Ww2XBdVYRu8bwZdt4fxPE/5fD58jYdNfNpqtZTP58P4isVi5KSY43wOGo2G8vl8WEHrum74fO/evdu3X28SNjhu3F9ZuK57baLYqzE2Go3wOQTHuun9GbyeVyefDarvg2PevXs3VtVwMHnj3bt3w/dg8Og9ZpzXJYg9lUrp7OxMnU5H9Xp94Be/g9qQNPA61XvtvPpZi3MtDz5Hw2IMjhl8KTTJftLs3rOL+iwAABDbosdRAQDgJr3jdUeNSRxM5DRo/azU6/WuaZpTay94Drlcrnt+ft7N5XIjTZoVxXGcvkm4RlGv17u5XK6by+WuTeI1je3H0W63u+VyOWw7GMO8Xq+HYwCbpjnxhIjBOZz2JJ7n5+fdTCbTNU0zcjKzQROM9U4OOa0xWZvNZteyrPCz0el0wuMMm4Cyd8LAQe/zYJtJPnftdrvvedu2HZ7zqxNwRp2jYF0wYWcmk+mWy+W+1z5qvOhgsrpcLhcuCz53ipgANlgftDnoczluu8H56X2OuVwuHHP/pudx9XNwfn5+7fjDNJvNge+5drs99NoavB5Rr0Ww76hx9L62wTU/+Gz2vgd64wjeK73vgUwm07dPMJ7yqJ+DZrPZ934sFArdcrnctSyrWy6XwzGTg+MF15BguyCeUSd3HfRaBM8hk8lEbhM890Gv4aD3Z+8Y/73nO2gv7hwPtm1HXtd65x8Y9xrbe22YZP6J3mvZoBiuvv96l497LR/U1lXBNSE4j3H3m9V7NgmfBQAA4iDZDQBIvN5EUKDT6YTJrSChMmwSv1kxTXMqScnehECQ5J4k8XB+fj5ykjxIzvUm4YKJMKMSXONuP4ngZrlcLnfL5XLXtu1uu93uS9hMI0kdJNGDpPc0J/U8Pz/vS9IPS5L0rp/G+zlI8lx9L/QmVQYl3rvd/uR71LkNEhrTEBzn6hc8QfJ20OvSm5C/+jyDhGHUcxw0WV23+83zvppoHCWZGKfd3lgNw+h7rXsTn1ff68Hn4Opr0mw2R04yd7v9n7NB66I+Z0ECOcq414Te1zaYCLD3OL3vj2H7BcezbTv8Ai7O52DYF09Bks8wjGuvc29idZLPcO+XvFFyudy1CQxHneQzuBYF5y54j00yqWW5XI68bgbv3XESnr1fAA1qd1y9Cdubthn02Rn1Wt77JVHUJJOB3mTxJPsFZvWeXfRnAQCAcZHsBgAkXu9N5dUqv6gE1zwFSZhJY+i9eb2pknOagvivJt+Cm9urz2vc7ScVvBa5XO5a21er26aht5JtWkmWq+33voevPqfe9/ZNxw6q3wuFwrVHkPALvjiJEiRsBq0P4g3iGZSgnVZF/LAvL3qTJleTVYOSMN3uN++Rq0mYIBE7KNHfm3TqbfemZGLcdrvdbl9y76re6tao40W9ZsO+xLgqSKxGvRd6q6yvikq4BoYl/aP0vrZR53LYF1zBZ2pQIjPO56C3Ev+q3s/FsATvJF+E9r4eUUntqOczarK79/MUVMOP836JEnW83n/XRvm34eqXjtMU/Fsy7D05SkI8MOxaPuzXEr2Cz2/wPoq7X2BW79lFfxYAABgXY3YDABKtd5xex3F0fn6ubrer8/NzFQqFcPLGRU3KlclkZNu2Hj58ONFklb0TcPZOStg7Buu0BeNuGobRd8xWq6VWq6VMJtM3bui420/K9/3wNc1ms9faDsYtneaEpJlMRs1mU+12W57n6e7du1OdiDSTyajdbod/976/xx3X1LIs2bYdjqVcq9V0dHSkSqUSTm7WaDTk+/61MXRTqVQ43qzv+wOP3XtOW61W3xi1QduFQmGsuAcJxtKNGlPXNM0wjmB8/lEE75G
      "text/plain": [
       "<Figure size 1800x600 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#in abhängigkeit von der energie der elektronen\n",
    "fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(18,6))\n",
    "\n",
    "\n",
    "ax[0].hist(up_energyloss_lost, bins=100, density=True, alpha=0.5, histtype='bar', color=\"darkorange\", label=\"lost\")\n",
    "ax[0].hist(up_energyloss_found, bins=100, density=True, alpha=0.5, histtype='bar', color=\"blue\", label=\"found\")\n",
    "#ax[0].set_xticks(np.arange(0,1.1,0.1), minor=True,)\n",
    "#ax[0].set_yticks(np.arange(0,11,1), minor=True)\n",
    "ax[0].set_xlabel(r\"$E_\\gamma/E_0$\")\n",
    "ax[0].set_ylabel(\"counts (normed)\")\n",
    "ax[0].set_title(\"Upstream\")\n",
    "ax[0].legend()\n",
    "#ax[0].grid()\n",
    "\n",
    "ax[1].hist(down_energyloss_lost, bins=100, density=True, alpha=0.5, histtype='bar', color=\"darkorange\", label=\"lost\")\n",
    "ax[1].hist(down_energyloss_found, bins=100, density=True, alpha=0.5, histtype='bar', color=\"blue\", label=\"found\")\n",
    "#ax[1].set_xticks(np.arange(0,1.1,0.1), minor=True,)\n",
    "#ax[1].set_yticks(np.arange(0,11,1), minor=True)\n",
    "ax[1].set_xlabel(r\"$E_\\gamma/E_0$\")\n",
    "ax[1].set_ylabel(\"counts (normed)\")\n",
    "ax[1].set_title(\"Downstream\")\n",
    "ax[1].legend()\n",
    "#ax[1].grid()\n",
    "\n",
    "\"\"\"\n",
    "most electrons lose little energy relative to their initial energy downstream\n",
    "\"\"\"\n",
    "fig.suptitle(r\"$B\\rightarrow K^\\ast ee$, $p>5$GeV, photons w/ brem_vtx_z$<9500$mm\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAHOCAYAAACVVnNvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABP5UlEQVR4nO3dTYzbWHo3+n9N7h3HgGGzyrPJ4gJd1Mwis4lNlTez6UybRK+nS3ItMkAW6RK74V1jWpwKArR70zIVbwszVDlALjBAUCW6g7uciEbcm1lcF9meTXLvOyOWB14ki7TENgw4NZlpvovqw9YHJVHUN/X/AYJdEj+OKIp8dM5zztkIwzAEERERUYZ8a9EFICIiIpo2BjhERESUOQxwiIiIKHMY4BAREVHmMMAhIiKizGGAQ0RERJnDAIeIiIgyhwEOERERZQ4DHKIZ8Txv0UWYqyAIUKvV4Pv+oouSOb7vo1arIQiCRReFFmzdriuTYICTQhAEqFar2NjYwMbGBjY3N1EsFqFpGnK5HKrV6qKLGKnVasjlclFZNU2D4zgALr4omqZFr/WW2zAM3qwSKhaL0XEUD8MwYpf1PA/FYhG5XA6bm5vI5/PQNA2GYaBWqyGfz4+1b8dxuj7jfD4P27b7lrNtG/l8PlpmmhfKWq2G7e1t6LrOm3AKuq4P/Nyr1SpyuRx0XUer1ZpzyWjafN+HruvQdR3FYnHodybpdcUwDBSLxeg+NOi6nXS5zAgpNUVRQgChZVnRc/V6PQQQFgqFBZasm2maIYBQluW+11RVDRVFCdvtdvSceD+maYbNZjN0XTdsNBrzKu7KabfboSzLYaFQ6Hq4rtu3bLlcjs6PZrPZtY1SqRQCCNN8LdvtdrRu5/nYy7KsUJKksbefhHhvce87rWaz2XVuZpUkSUM/t0KhEALoOmeyLoufveu6fd+RcrkcSpLU99kmva4oitJ1v2k0GqEkSamXyxIGOBOQJCn2opP2JjUr4sbZewEtFAphqVTqW9513bBUKoWlUik0TTMsl8uZu9BMU7lcDuv1+sjlxE3KNM2By1iWlfrcEQGGqqoDlxGf6SyIQHqaF0xFUTJ/U280GiGAod8x8R3O+rHolMXPXlGUUFGUvufjfhQnua6I73zvudO7n6TLZc3y3IVXjIjE42pFRIAzzaBA1KakIcty18Wx2WyGiqKM/KWvqmpYLpczE+G32+2R7zsNSZJC0zSHHidRs5fkYpL2gtNZizOoLJIkzSxYnXaAsy61FoVCYeRnvm4Bzqp99vV6fWRA1mw2QwCxPypVVe17v0muK5Ikxd6Der+LSZfLGubgpCTyWFRV7Xpe5D4oigJJkqa2P1mWYZrm2OsFQQDf9yFJEmRZhm3b0DQNR0dHKJVKsesUi0XIsgxN03BwcIDj42PUarVJ38LCSZKEx48fo9lsTi1XqlqtIggCGIaBfD6Pzc3N2PwX0W6u6/rIbdbr9b7nRN6OyPOKa4eXJAmFQgEAUKlU+l63bRuqqo51XorE4Xw+D8dxov+LvLO43IFWqxUtt7GxgWKxGLtdXddhGAY0TevKDRNlFTlCIlehM2do1Pqd5bZtG47jJC6PruvI5XIjz3nHcaK8CLFNz/OifKjePCeRT9Gb92DbdqLzorOcm5ub2Nzc7FtPfL/FZ9W7zLDzSBwzTdOiZHFN07C5uQlN06LPWuQEbW5uDswzi+N5HjY3N6Nj03kcDMOI8gAHffaiPOKY27aNIAiiMo7zfc7lctA0LcqB6cx1GSc3TeQ4NhoNPH78GLIsD33/AGK/f2I9sUyS64rneQiCAIqi9G1PPHd8fJx4OWDyc2DW59DYFh1hrSoRcXdWIXb+Sp/FL4+4CHwU0eShqmpYKBRCVVUT/4JPWntTr9fDUqkUlsvlge89yTJJ1Ov1UFXVUJKksF6vh5ZlRblQ4+Y9maYZyrI8URNcu90O6/V6WC6Xo5oy9PxKE7/cAKTKZXJdt6vZSZxncb8EO/cVVx097v4784JEjV69Xo9+YXeek+LXoKqq0S9PsX5ns5jrun1t/+I87VxOVKv3nitJ1m82m1EZO2si48oThhc1BuVyuWt7SZryRBk7j6tocurcnthm73Pisxx1/olyq6oalkqlrvNenBv1ej06Bzu/a6J2aNR51Gw2u/Yjjpl4P2LfjUaj6/iO8+t/UI5is9nsKtugzz4M+2ukVVUduwai97sjzt2471ScNNcOcRzjrlPi/YpzLsl1RZzzg9IMxL6SLheGk58D8ziHxsEAJyVxwokLCIBQkqSZJuOKm/s4xAkk8oXGyb9IEoSUSqWuL2xcsmSSZcYhvkCdAZi4QPTeQJKwLCuUZTkslUoTN9+IC3jnTU98uYfdyFzXDU3T7HqI9RVF6bsAiM8zbnsi+O48Fs1mM3Vysbj4935mvUH+oOU6b8Li/cSdx+J7JM67QTe5pOuLz6L3nO8tTxheHM/ecyfJd0UElL03rbgmgd7EcvFcku+0OOd7czLEzU+cK+IziPseJDmPxA2v92bYe2w7lx03p0uUuZNpml3vbViA03vzTvOd77xOi6bdUc237XY7CjjS5LF1NiH3Eu930HUx7roy7LMW56WiKImXEyY9B+ZxDiXFJqoURFW4oihwXTd6AOirKp+mQqEQVccn7YoryvL48WMAF1XBSdcdVt0KXLzXk5MTHB0dAUBUZdzZbJdkmXGdnp4CAPb29qIq1r29PQBI1e2xVCqh2WxC0zTcvn0buq6n7j5ZKBTQaDQAAJZlJV5PURSUSiUYhhFV2aqqCt/34XkeKpVK1L2zs4lFHItOYv3OJhbTNHFwcJDqPQm954No+hDvV9jZ2en6W5Kk6HiK9xNXXS62N+y4pVk/rkmgt7u1LMuoVqtdzRzlcnlgOTrXUxQltlnS9/3ofQdBgCAI+o6hbduxTWaD9L7v3s9AvNdbt271lWWc86j3mIlyb21t9T3XbDYTl7+zzJ3H7Pj4OGpeHUV8V2zbhmmaqZruO68/+/v7AC6+I3HniujWnc/nkcvl0Gw2E50bvSRJitICOo+9aE4EBl9z464rSZqagyBIvFxceTuNew7M8hxKigFOCuJE6/ySKIoS3cRHXbBEG3+ah2VZcBwHm5ubIwMp0fYqLsLiSzmNNs9arQbHcXB0dARJklCr1WAYBhqNRnTSJlkmDc/zUCgUui724kbS+eUZV6FQgOu6XTeANGO6qKoKRVGiMnXe8IcFTpIkRe9J/Cva5Ov1etej3W4jDMPYQFFVVciyHLWHA8DJycnAnKu0RBnHCQaH5TeI4zRse5OuP0i9XockSTAMA7lcbqw8jN4bdq1Ww507dwB8czM6OTnpuy6I5cWyaYjPv/c9995c0pxHsyLOQ3FsPM/rC4pHSRPUxHEcB7ZtR0FTr84xY5rN5sTfIcuyYJomPM/rykUS161hx6H3uiKuoXHXKPGcLMuJl8siBjgpiMBC07Su58VFQvxaG8SyLDSbzVQPy7KgKAoajcbIi1JvIrT4hTKN0WbFRd33/ShgajabXWVKssy4xE2hN7ny6dOnAEYHl+OYZMA6WZaji1bnzSauxiWOWLez5mMc4nibpgnbtrGzszPVpPfOMqa5OA77xZgkSJ10/V6yLOPs7CyqNcvn84kT63tv2OImpqpqtA3LsvpujsfHx2MnffcS6476DNKeR7MgkuEdx4Hv+zg+Pk6cZC2IjhO+70/UWUDsNy6xf1bK5TKazSbCMITrutjZ2YHv+ygUCiPPhc7rivjM4wZ/FM8pipJ4uSxigDOmIAiiX0O9N+rOi8e0bybAxS8dwzDw+PHjREGCqGnqDMTEL59xLyidROBUKpVQLpdhmmbfxTvJMmnE1Z4BF7+aVVWdOHjK5/OwLAv1eh2NRiP159jba0Ec93F/eYqLU1wTCICBtXilUim6Aezv70/0eQ8iLo7jjLwsjklcuUXQksvlZrb+IOKG2Wg0opvdOMdM3LA7g0kxQm21Wo39ZT5u81Qc8Z57m6R6pT2PZkU0l1qWNbDJcZj9/X24rgtZllOPuC7WK5fLAwNEUdPVaDQS9axLo/PHyCid1xV
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(up_energyloss_lost, bins=100, density=True, alpha=0.5, histtype='bar', color=\"darkorange\", label=\"lost\")\n",
    "plt.hist(up_energyloss_found, bins=100, density=True, alpha=0.5, histtype='bar', color=\"blue\", label=\"found\")\n",
    "plt.xlabel(r\"$E_\\gamma/E_0$\")\n",
    "plt.ylabel(\"counts (normed)\")\n",
    "plt.title(r\"$B\\rightarrow K^\\ast ee$, $p>5$GeV, photons w/ brem_vtx_z$<9500$mm\")\n",
    "plt.legend(title=\"LHCb Simulation\", title_fontsize=15)\n",
    "plt.text(0.35,2.0, \"Upstream\", size=15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAHRCAYAAACGvdZwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABSTElEQVR4nO3dTYzc5p0/+G9nZqMxVpDY7bnksrBYySGnyKzWWUhEwrcBki6qDxNgLu6iHd2McdE9wMI2FnCbNbpqYrI1wBxyUZPynHaBDKkZeE8LqEkrwOL/ByYoShksNv9DVEULWniUyQz30HloVhXrjfXSVc3vByhIXcWXp55i8fnV87qVpmkKIiIiogr61nkngIiIiOi8MBAiIiKiymIgRERERJXFQIiIiIgqi4EQERERVRYDISIiIqosBkJERERUWX963glYR4ZhoFar4fnz57hx4wYajcZ5J2npqvieiYiIGAgN0HUdsiyj1WoBAOr1OiRJgqqq55yy5anieyYiIgKALc4s/Y04jlGr1dDpdCDLMgCg3W7D9334vn/OqVuOZb3nIAiws7MDRVEWlVQiIqKFYx+hnCiKACALCABAURQEQYAkSc4pVcu1rPecJAm63e68ySMiIloqBkI5jx8/hiRJfc/t7OwAwIUo1E3TRBzHfc9d9Pd8kYigtSqSJIHjOEPXLM0vjmM4jnNhf+DR9Kp2Xyly4QOhJEnQbrextbWFra0tbG9vQ9d1aJqGWq2Gdrvdt60IAgaty83YcRzUarXs/WiahiAIAJxd0JqmZa+J9+Y4DgDg9ddfz7YT+2zCe64iXdezz1E8TNMs3DaKIui6jlqthu3tbdTrdWiaBtM04TgO6vX6TOcOgqDvGqvX6/A8b2g7z/NQr9ezbRZ5Q3UcB9euXYNhGCysSzAMY+Tn3m63UavVYBgGf+xcAHEcwzAMGIYBXdfHfmemva+Ypgld17OyclRZMO12ay+tCEVRUgCpbdvZc67rpgDSRqORpmmaWpaVSpLUt18YhimAtNPprDS941iWlQJIZVkeek1V1VRRlLTX62XPhWGYNpvNtNlsppZlpa1WK3t9We/Zdd3U9/3S+1dZr9dLZVlOG41G3yMMw6FtW61Wdg3nP69er5c2m80UQFrma97r9bJ989+ZQbZtD10/iyLeW9H7nken0+n7flxEkiSN/dwajcba3deW7SJ+7uJenf+OtFqtVJKkoc922vuKoihZmZimaer7fipJUuntNkFlRo2JSDU/EkoMERe/dmVZHoqkxS+mfB+a89bpdABgKJIXo78GOzkrioJ6vQ7XdSFJEvb397PmsEW9Z8Mw+v4+PT3Fzs4OXNfte9627amPCZzVDIRhOHYbTdMu1HD/o6MjWJY18T3pug7P82BZVjbiT5AkCbZto16vD30205AkCa1WC+12G67rotlsFm4XhiEODw9nPv40RA3moum6nn0XLiLRv+/27dsjtxlVC3yRXcTP/eDgAIqi9A1KsSwL7XYbpmn23X+nua+YpokoivDo0aPsOVVVIcsyDg4OsnvxtNttjPOOxFZBRM1FNSj446/eXq+XdjqdoV9JlmWliqIsND2WZc31S0yW5b50djqdVFGUkb8AG41G6vt+allW2uv10larlW27rPdchRqhXq83Nt/LkiQptSxr7C8rUZs5zedU9rPM1wqNSoskSUv7lS1qPhf5C7MKNSGNRmPiZy5qCy9yPuRt2ufuum6qKMrY9Ip7d7PZHHpNVdWh9zvNfUWSpMJycvC7OO12m6ISgZD4cAYvmKLCpNFopJZlZX8rirLwAt113cKLdxqicBLNEa7rprIsT3XhtVqtwu2W8Z6rEAilaZoFlrIs9+VhWeJaFQ9JklLXdYe2E8HwNEFY0c00DMO00WikqqqmsiynrVarcF9RgOSrwAXXdQufH6fX66W2bWfXmPi/JElpo9HoC6pEXuS3G5UW0RTYarVSVVVTVVWHrj/xXQGQqqo61Cww6Rj5tIvre1Sa8sdqNpupLMsTPyvf97PPXRwvDMMszYqi9KV3VOE+zXUhAiHRbC5JUipJUuE9UuSDaAbNbzPuOhL5papqatt22ul0UlVVU0mSUlVV+5rnZVlOJUkaeR0WCcMwlSQpy5t8PohmVcuyRn7uIj0iz13XTXu9XpbGWb7Psixnx242m9lnM2tQYNt2Ksty2mw2J/7AEOVXUZ6Jz1fcO6a5r4gKg6Lvl7g2RRkyzXZpOv81sOxrSKhEICQu9vwHnw+CBm8k+b40i/61LxRF09OwbbvvC52/GCYZ9+ti2vcsgrhWqzX2F8usgdC0x100caMXN4ZJBe444suY74M1q16vl7qumwVX4saVL3zEL0ERJMwqDMNUVdXsb/FdKArO8+cafE9lAuZ8vyVVVdNWq5UFVIO1tuLmrapq9ktW7J8vpESBmC9wxPdksDATBeTg9TXNMTqdTpZOkfZRaWo0Gn03ZNu2pypYRfry+TpYuOSPOfic+CwnXX8izaqqps1ms++6F9dGPoDIfzfFD8dJ11Gn0+k7j8gv8X7EuX3f78vbWQKHwX6egigwhVGfe5oO17Crqjpzjcbgd2fUj+9Rytw7RD4W3afygWCaTndfEdd7UZrzwc+026Xp/NfAKq6hNK1IICQ+dPElFhHxedZYiAJ4VuKDFr+EFlELMS3xS0cY1yFzlkBoluMug/ii5WvMxI2kzK+LWX7VTSJu9PnCMV9zMOr4YRimlmX1PcT+gzULaZpm11PR8cQPiXxedDqd0p2kRSEx+BkP/mAZtV2+sBbvp+i7JL7rRTUFgwXitMcQn8fg924wTUW/TKf5rorAc7BwK2qKGOwgL56b5r4yWGMgiEJSXCviMyj6HkxzHYmCcbDQLPpsxLaz3tNEmvNETZAwLhAaLOTLfOcHaw7FPXrc93/e2uR80/Ug8X7H3Z8H7yvjPmtxXSqKMvV2wrzXwCquoQs/fF4ME1cUBWEYZg8AfUPPV63RaEDTNGiaNtPwYJFe0UnNNM2VDC/WNA0nJyc4Pj4GcNbBPEmSkctwqKo61RIdsx53GU5PTwEA+/v7WafD/f19AOWmEGg2m+h0OtA0Dbdu3YJhGKWHlTYajazz+ywdzRVFQbPZhGmaWad6VVURxzGiKMLR0VE27FXX9Ww/kRd5Yn8xDQNw1iFz3k7Sg53xRafuwc7+u7u7fX9LkpTlp3g/RTOYi+NNyrcyxyjqcJsfii7LMtrtdt/0HIMd2ovIsgxFUQqnK4jjOHvfSZIgSZKhPPQ8r+/znGTwPQ9+BuJ93rhxYygts1xHg/kl0p3vtC2eE4NBpiXSnM+zBw8eTD2AQnxXxMADy7JmOj/QPwjn4OAAwNl3pOg6EcPd6/V6Nqv/NNfGIEmSskEM+bwPgiD7ro4a8FJ0X5mmE3mSJFNvV5TevFmvgWVeQxc+EBIfdv5CVRQlK3hnuWkUEYuVlnnYto0gCLC9vT1VQBZFUXbzUxQl+/KMml9mURzHQRAEOD4+hiRJcBwHpmnC9/2RX7RpvixljrsMURSh0Wj0FQqiwJlndE2j0UAYhn0FRZmgVVVVKIqSpSkfGIwLsCRJyt6T+FfM9eO6bt+j1+shTdPCAFSMBhETHALAycnJyJFkZYk0zhI0jpu7SOTTpOMt4hiDxOgk0zRRq9VmmmNpsGB3HCcbASYKrZOTk6F7l9h+3GixScTnP/h+B7/PZa6jZRHXocibKIqGgudJygQ/RYIggOd5WXA1KD/nTqfTmfs7ZNs2LMtCFEXZfF5xHGf3rXH5MHhfEffconuUeE6W5am32yQXPhASAYamaX3Piy+q+GVVlm3b6HQ6pR62bUNRFPi+P9WNQ7wXsa34xbHs2XfFjTmO4yzo6nQ6c9/slnXcWYjCY3CI+ePHjwHMHyjnzXOdybKc3dzyhVJRDU4RsW++JmUW4vOxLAue52F3d3fhw5BFGsvcRMf9Ap02mF3EMQRZlvH06dOsFq5er/fVqI0zWLCLwk5V1ewYtm0PFaIPHjyAqqpzfS75aTXGKXsdLYMkSWg0GgiCAHEc48GDBzNPGRHHcVbTmK/Fm5U47+C
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "both_eloss = np.append(up_energyloss_found,up_energyloss_lost)\n",
    "plt.hist(both_eloss, bins=100, density=True, histtype='bar', color=\"cornflowerblue\", label=\"Upstream\")\n",
    "plt.vlines(ak.mean(both_eloss),0,3,colors=\"red\", label=\"mean\")\n",
    "plt.xlabel(r\"Energyloss Ratio $E_\\gamma/E_0$\")\n",
    "plt.ylabel(\"counts (normed)\")\n",
    "plt.title(r'$B^0\\rightarrow K^{\\ast 0} e^+e^-$, $p>5$GeV, photons w/ brem_vtx_z$<9500$mm')\n",
    "plt.legend(title=\"LHCb Simulation\", title_fontsize=15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABjIAAAJPCAYAAADWqU3EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAACkIElEQVR4nOz9cYxbWX7Ye/6u1JZGdqebVZrkxQO8yTTLs4FfFnLjUuO38CLp8WsS6cUM9kVPVZIlLdbewCL9UL3OevJcjBJkNYLhLbOQbLDBFGyWsrvzx7aUFikIuxnv04bsSc++xIDRItMWvLH9Zord4933vDCmi5x+vVZLo9LdP6ovu1hF8vwuee7luazvByC6Vffw3N8959zLe3l4zvGCIAgEAAAAAAAAAADAQcdmHQAAAAAAAAAAAMAodGQAAAAAAAAAAABn0ZEBAAAAAAAAAACcRUcGAAAAAAAAAABwFh0ZAAAAAAAAAADAWXRkAAAAAAAAAAAAZ9GRAQAAAAAAAAAAnEVHBgAAAAAAAAAAcBYdGQAAAJg77XZ71iFgDvR6Pdna2pJOpzPrUOZOp9ORra0t6fV6sw4FM8b1GgAAaNCRAQAAUqXX68nGxoZ4niee58nCwoKsrKxIoVCQpaUl2djYmHWIfVtbW7K0tNSPtVAoSLPZFJG9L24KhUJ/28G4y+UyX54qrays9MsxfJXL5aFp2+22rKysyNLSkiwsLEgul5NCoSDlclm2trYkl8tF2nez2Ryo41wuJ/V6/VC6er0uuVyun0bzxV29Xh84toWFBSmXywNf/PZ6PSmXy0PbGKaztbUlL730kpRKJb5sn0CpVBp5Pm1sbMjS0pKUSiXZ2dlJODLY1ul0pFQqSalUkpWVlbHnjPZ6XS6XZWVlpf/5PurzUJsOAADMgQAAACCFfN8PRCSoVqv9v9VqtUBEguXl5RlGNqhSqQQiEmSz2UPb8vl84Pt+0O12+38Lj6dSqQTb29tBq9UKGo1GUuGmTrfbDbLZbLC8vDzwarVah9Kura3128f29vZAHsViMRCRYJLb426323/v/vZ4ULVaDTKZTOT8M5lMICJBrVYbmaZSqQxtY7Zsb28PtNOjImwzw9rTpI5KWWYymbHnw/LyciAiA+fivJvHum+1WofOkbW1tSCTyRyqW+312vf9gc/xRqMRZDKZidMBAID5wIgMAACQSuGvLvP5fP9vy8vLIiJDfxE/K9vb2yIih35xurKyItlsVlqtlmQymf7fz549K6VSSba3t6Ver8ubb74pZ8+eTTLkVFlfX5dKpSK1Wm3g5fv+QLqVlRXZ2Njop81ms/1tmUxGqtWqVKvViWLIZDKytrYmIiK1Wm1kularJdeuXYuc/+LiYn8/o/i+f+iYbVpZWTmSv5w/ffq09TyPQlk2m03p9Xpy4cKFkWnCdn2UzGPdX7169dD1p1Kp9EeL7ae5XpfLZWm323Lz5s3+3/L5vGSzWbl69WrkdAAAYH7QkQEAAFKn3W5Lr9eTbDY78IX0fjangtnY2Jh4uopwmp+ww6XT6fSnMxr2xbnv+5LL5aTT6cgHH3wgFy9eHPsFdlr0ej3J5XKytbVlNd9w/YJxUzXV63Wp1+vi+36/w2GYYrE4cWdA2EHRbDZHxnLnzh0pFosT5T9LKysrzGFvyVEpy2q1Kr7vz8W1y5a01X04Hd64z77w2jussz2fz0u9Xh94v+Z6vbW1Jdls9lDbuXjxorTb7f57tekAAMD8oCMDAACkzsHOgVA4EsP2F2jZbFYqlUrk9/V6Pel0OpLJZCSbzUq9XpdCoSA3b94c+YV2OFKjUCjItWvX5M0337T+5f8sZDIZeeutt2R7e9vaWiYbGxv9X/3mcjlZWFgYOhon/FVwqVQy5jlsREW4rka4Dsuw+dwzmUx/RND6+vqh7fV6XfL5fGxf7B78dXu4SHUul5Nms9n//3BNmYPrbJRKJSmXy1IqlWRpaanf5ur1ev8LwXD++/DfYXsO819YWBgoY025bWxs9Pd9sKMrPIZCodD/ArRQKMjCwoIUCoX+MYTrLYRriJhEKZvQzs5OP53nebKysjI03/BYCoXCofVKxpWl5v37467X69JsNtXxHKzXUZrNZn/dgjDPdrvdXwfm4Pou4XoHB7/srtfrqvNtf5wLCwuH2lCY16TtLK42tH/fCwsL/bLZXw7h2jUbGxsj6z6MJyzzer0uvV6vH2OU6+TS0pIUCoX+GhX716KI8qV+uLZTo9GQt956a+SPBcLjFxk+Wix8X5hGc70Of6QwrEM5/Nubb76pTmej/uNuQwAAIIJZz20FAAAQVT6fP7RmQLg+hu/7scy5Psn6A9VqNRCRIJ/PB8vLy0E+n1fPj762tqaa57tWqwXFYjFYW1sbeeyaNBq1Wi3I5/NBJpMJarVaUK1W+2uVRF2XJFzTYW1tbeI547vdblCr1YK1tbUgm83216koFov9NNvb2/2/T7LWSKvVCvL5fP/fYTvbv49h+zp4TL7vT7zWSXhs497farUG6mD/mh/5fD5YW1sLarVaf12C/e15eXk5WFtb6/+7Wq0GlUql/+9wnYj97aZWq/Xj2t+2fN/vx2MqtzDfUKPRGDjO7e3t/nGExxCuGRP+rVgsBo1GI9je3u4fm+m8iVI24Ro3+Xw+qFQqQavV6r9/fxm1Wq1Dc/OH57+pLLXv33+M+8tjWDxBYK7XUcIY97e3sMz35xfmefBvYV2bzuv9dVssFgeuJ2HbmbadxdWGhh3vwWvg9vb2QGyj6j4IPj3Hw235fD7yOg8Hr0lh2x12rRpmkmtyWI7Drv/h8YZtTnO9Dtv8sJjDtTiWl5fV6WzUfxJtCAAA6NCRAQAAUif8AiT8QktEgkwmE+ui2OGX+FGEX2iEizVrvkQMaTobisXiwBdIwxbX1aSJIvxCZ39HS/iF1cEvNDWq1WqQzWaDYrE49SK44ReK+7+EDb9sGvfFaqvVCiqVysArfL/v+4e+kArrc1h+YSfb/rLY3t6eaJHv0CQdGUHw6ReZB+v7YEdgJpM5VHeaL9/D/IfVu6bcwk6xULho+v78wi8lD35hGZ73+2MK02rOM23ZjEq3/8v2MJ5h14eDcY4qS+37wzZ+8BgPxhME5nodJeyQO9ieMpnMoQ7d8Mvig3/TXCvDa8nBRewPtvdp21lcbWhYzPtVKpWBYxvXkXHwS/pJrqX7rw/huZTJZMZeV7vdbr9jIeox79/PsN9Hhsc76vNm2PV6XF2H7dL3fXW6ILBT/0m0IQAAYMbUUgAAIFXCqVZ835dWq9V/icihqVhsWl5e7k/3ol1/I4zlrbfeEpG9qUa07x03nYfI3rHeuXOnv9BpOCXJ/um2NGmievDggYjszUMeTuFx8eJFEZGJ1hEpFouyvb0thUJBXn31VSmVShOvR7K8vCyNRkNEJNLC3b7vS7FYlHK53J8SJJ/P9+dyX19fl5WVlf4rFJbFfuH790/hU6lUJlrk25aDbSmclicsq2w2KxsbGwPT2IxbSyQUTifzpS99aeDv2nKr1Wr9c3f/34edIwenrgmPaf+UWuHftre3jbEffE/oYNmEDq4BkMlk+u00PN5h09yE+Y1rj5O8f9hUPgcXkZ60XrPZrPi+P3Satk6n0z/uXq/XX6tov3q9PnSqq1EOHvfBOpi2nYXiakP7Y95fZm+++WZ/ujmT8BpUr9elUqlMNJXh/ut6uNh1pVIZ2lY6nY6USiXJ5XKytLQk29vbqrZxUCaT6U+TuL/sw2nAREZ/lg27Xmum3uv1eup0B2Pdb5L6j7MNAQAAMzoyAABAqoRffOz/0sb3/f6X9aYv0MK54id5VatVaTabsrCwYOww2b8g+f5Fpm3Mnb21tSXNZlNu3rwpmUxGtra2pFwuS6PR6H+JokkziXa7LcvLywNfPoZfbB5cpyGK5eVlabVaA19ITrJgez6fF9/3+zHt/wJ6XAdJJpPpH1P433Bu91qtNvDqdrsSBMHQDqF8Pi/ZbLY/r7pIMot8H/wSe5zw+MLyqNVqkslkpFwuy9LSUuRFcg9+uactt/1rx7iyEPLBstEYF3fY/sb
      "text/plain": [
       "<Figure size 2000x600 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#energyloss in abh von der energie der elektronen\n",
    "#upstream\n",
    "fig, ((ax0, ax1)) = plt.subplots(nrows=1, ncols=2, figsize=(20,6))\n",
    "\n",
    "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",
    "ax0.set_ylim(0,5e4)\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(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",
    "ax1.set_ylim(0,5e4)\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\"$B\\rightarrow K^\\ast ee$, $p>5$GeV, Upstream photons w/ brem_vtx_z$<9500$mm\")\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABjIAAAJPCAYAAADWqU3EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAACX70lEQVR4nOz9YZAjZ37fef6yu6d7WqJJVPWc15qInR2iPKtQ2MeYSPTIt76zObMErFFIu3a7q5oieXvy+URAF8WVT+NVQW3vRk+vQltChbWKi1CFhGrfxry4IY8NdHR4LZ1pARzN2D5H7LEBjzp8K3tHAEe6W+vCHlZheLR7uqer814UEyxUAcgngQeJJ4HvJwJBduWDJ//55JOJJ/HgeR4vCIJAAAAAAAAAAAAADjoz7wAAAAAAAAAAAABGoSMDAAAAAAAAAAA4i44MAAAAAAAAAADgLDoyAAAAAAAAAACAs+jIAAAAAAAAAAAAzqIjAwAAAAAAAAAAOIuODAAAAAAAAAAA4Cw6MgAAAAAAAAAAgLPoyAAAAIDT2u32vEMAjPR6Pe3t7anb7c47lIXT7Xa1t7enXq8371AwZ3wmAACwnOjIAAAAc9Pr9bSzsyPP8+R5nlZWVrSxsaFCoaC1tTXt7OzMO8S+vb09ra2t9WMtFApqNpuSjr5UKRQK/W0n4y6Xy3yxaWhjY6NfjuGrXC4PTdtut7WxsaG1tTWtrKwol8upUCioXC5rb29PuVwu1r6bzebAOc7lcqrX66fS1et15XK5fhqTL9WazeapY1tbW+vHXigUtLOzw5e0Kba3t6dnn31WpVKJ8ziBUqk08prd2dnR2tqaSqWS9vf3E44MtnW7XZVKJZVKJW1sbIy9Zkw/E8rlsjY2NvptiFGfuabpAACAgwIAAIA5830/kBRUq9X+32q1WiApWF9fn2NkgyqVSiApyGazp7bl8/nA9/3g4OCg/7fweCqVStDpdIJWqxU0Go2kwk2dg4ODIJvNBuvr6wOvVqt1Ku3W1la/fnQ6nYE8isViICmYpKl7cHDQf+/x+nhStVoNMplM7PwzmUwgKajVaqfyy2azgaRga2srdr6u6XQ6A9fCsgjr5bA6O6llKctMJjP2mltfXw8kDVzvi24Rz32r1Tp1jWxtbQWZTObUuTX9TPB9f6Ct0Gg0gkwmM3E6AADgJi8IgmAuPSgAAAAfWllZUa/XU6fTUTab7f/d8zxJkivNlVKppL29PVWrVRWLxf7fNzY2tLq6qmq1OpC+3W73/7a2tqb33ntP169fVyaTSTLs1CiXy/rc5z6n9fX1sek2NjZUr9dVqVS0tbU1NM3e3p5KpdJEdadcLmtnZ0f5fF6NRmNomlKppLW1tZH7H2VtbU3dbleNRkP5fP7U9nCER7FYPFWf0iSXy6lWqw1cz8tgZ2dH5XJZrVZLvu9byXMZyrLZbKpQKOjg4GDk/TG8/578nFhki3juw1E3rVZr4O+e52l9fV21Wq3/N5PPhPB+fbLunNyPaToAAOAuppYCAABz1W631ev1lM1mR35ZY3Oalp2dnYmnkginkgq/gO52u/3pjIZ96ez7vnK5nLrdrt577z29+OKLC9GJ0ev1lMvltLe3ZzXfcG2BcVM11et11et1+b4/thOhWCxO/EXy9evXJR2d71Gx3L59e6Azy5a3335bmUxGe3t7qZ0HfmNjI7Wxu2ZZyrJarcr3/YW4P9qStnMfTrk37vM1vL9fvnz51LZ8Pq96vT7wfpPPhL29PWWz2VN158UXX1S73e6/1zQdAABwFx0ZAABgrk52DoTCtQlsf7mVzWZVqVRiv6/X66nb7SqTySibzaper6tQKOjWrVsjv9De2NhQNptVoVDQ9evX9eabb1r/8n8eMpmM3n77bXU6HWtrmYTrQ5TLZeVyOa2srAxdnyKcG71UKkXmefyXvaFwXY1wHZZhc61nMpn+L4C3t7dPba/X68rn8zP50jWTyfTr06uvvjqwrdfrqVQqqVwuq1AoDKzTIh1dS+Ec8hsbG5KOjjdc9+Pkeh7h3PPdbre/SHW4Lkiz2eyvAxLmNSyOcGRKWK/r9Xp/H+H89+G/w2um2Wxqb29PKysrA+fR5Nzs7Oz0932yMy08hkKh0P8CtFAo9NcgCTtEw/UWVlZWRq6/cvJ4w7IJYw/r6MbGxtCO1v39/X66YWV4shyHnc9xZWny/jjn9GR+J8/rKNPUuePq9brRNX08zpWVlVN1KMxr0no2qzp0fN8rKyv9sjleDuVyub/G0qhzH8YTlnm9Xlev1+vHGOdevLa2pkKh0F+j4vhaFHG+1A/Xj2o0Gnr77bfHjh4J8x127wzfF6Yx+UwIfwgxrNM6/Nubb75pnG7a8z/r+gMAwNKb78xWAABg2eXz+VNrBoTrY/i+P5P50IetcRGlWq0GkoJ8Ph+sr68H+XzeeO7yra0tozm4a7VaUCwWg62trZHHbpLGRK1WC/L5fJDJZIJarRZUq9X+WiVx1yWpVCpBNpsNtra2Jp7P/eDgIKjVasHW1lZ/rQhJQbFY7KfpdDr9v0+y1kir1Qry+Xz/32E9O76PYfs6eUy+70+81kl4bOPeH84hf7yp3mq1Ts3lHtbJSqXS/1u4RsPx/BuNxtC1N6rVav9vnU6nvwZBPp/v19lwvZHj+1hfXx/Iq1qtDo3heN2s1Wr9Yz9ef33f7x9f1LkJ8z15XOGxdjqdfrzHjyFMl8/ng2KxGDQajYHjjbo2j6+5EuZbq9X67z9+PwnX0cnn80GlUhlZhnHP58nr3OT9cc5pEESf11EmrXOh8FxH3TuOn9tisThwzwrrzrT1bFZ1aNjxnrzPdjqdgdhGnfsg+Og+Em7L5/Ox13k4ed8L6+6w++Ewk9z3w3Ic9hkTHm9Y50w+E8I6Pyzm8D66vr5unG7a859E/QEAYJnRkQEAAOYq/HIi/LJJUpDJZGa6KHb4JX4c4RcO4WLNJl/whUw6G4rF4sCXO8MWvjVJE0f4hcvxjpbwy6RJFpwOF6wuFotTL1Abftl3/AvS8MugcV96tlqtoFKpDLzC9/u+f+oLo/B8Dssv7GQ7XhadTmeiRb5DJh0ZxxccD+PyfX9onQ2vmbCOhR0wJ78ozGQypzrwTi6UHpb5ybp9/IviMK+T9cPky/fwi9Jhdcvk3IQdb6GwnI7nF34pefILy5PldDytybUcxn7yejvZETsq3ckyND2fo8rS9P2m5zQIos/rKNPUufBvJvfj8H51vNM7CE5fU9PWs1nVoWExH1epVAaObVxHxskv6Se5Xx+/B4XXUiaTGXvvPjg46HcsxD3m4/sZ9nvK8HhHfaYN+0wYd67Deun7vnG6IJj+/CdRfwAAWFZMLQUAAOYmnAbF9321Wq3+S9KpaVJsWl9f70/FYrr+RhjL22+/LeloGhDT90Yt1FooFHT79m3dunVLkvrThRyfbsskTVz37t2TdDRHeDi9xosvvihJE60jUiwW1el0VCgU9MILL6hUKk28Hsn6+np/oe04i177vq9isahyudyfsiOfz/fnWd/e3tbGxkb/FQrL4rjw/cen16lUKv01NJKwv7/fj33YtCjhlDlhGWWzWfm+P3Rarm632z8fvV6vvzbNScOmfdnf3+//fzab1c7OzsA0NiaLnof5fu5znzsVl8m5qdVqAwvyhn8fdh2ePIbwOFdXV0/9rdPpRMZ+8j2hsPxPLgp/cg2ATCbTL/s453OYSd4fdU6lyc/rtHWuXq8PnepqlJPHffIcTFvPQrOqQ8djPl5mb7755thFrY8L73P1el2VSmWi6RKPf3aE09hVKpWhdaXb7apUKimXy2ltbU2dTseobpx0fOq842UfTgMmjf68HPaZYDK9X6/XM053Mtbj4p7/WdYfAACWFR0ZAABgbsIvJY5/oeL7fv/L+qgvt8J53Cd5VatVNZtNraysRHaYHF+Q/Pgi0zbmtt7b21Oz2dStW7f6izyXy2U1Go3+lxwmaSbRbre1vr4+8MVg+KXj8S9b4lpfX1er1Rr4snCSBdvz+bx83+/HdPzL4XEdJJlMpn9M4X/DeddrtdrA6+DgQEEQDO0
      "text/plain": [
       "<Figure size 2000x600 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#downstream\n",
    "fig, ((ax0, ax1)) = plt.subplots(nrows=1, ncols=2, figsize=(20,6))\n",
    "\n",
    "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",
    "ax0.set_ylim(0,5e4)\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(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",
    "ax1.set_ylim(0,5e4)\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\"$B\\rightarrow K^\\ast ee$, $p>5$GeV, Downstream photons w/ brem_vtx_z$<9500$mm\")\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABkMAAAJPCAYAAADPIqNiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAACp7ElEQVR4nOz9f5Aj6X3feX6ye8hW743Z6KLsG/IubBVqGSHLskUmmjJHveYcOYB1Dt/KKzWq21qfbe4eByVvURZjJlRwa+/cbEWsS6jQcCmLdSbQf5gXsQ6rG2h5bYd2wwSGjKG3Nb5gAxpZ1o+zVChb55W4PqoKPR6r2cOuzvujmOj6gR/PF5UJJLLerwjETBe++eDJJ5988GQ+ePLxgiAIBAAAAAAAAAAAkFJnZp0BAAAAAAAAAACAODEYAgAAAAAAAAAAUo3BEAAAAAAAAAAAkGoMhgAAAAAAAAAAgFRjMAQAAAAAAAAAAKQagyEAAAAAAAAAACDVGAwBAAAAAAAAAACpxmAIAACIVK/Xm3UWAAAAAAAADmEwBAAARKrVaqnVas06GwAAAAAAAH0MhgAAAAAAAAAAgFRjMAQAAACIWafTmXUWkBK9Xk+1Wk3dbnfWWUmVbrerWq3Gox4hiTYbAIC0YjAEAAAkSqfTUblc1vLyspaWlrSxsTHrLI3U6/W0sbEhz/PkeZ4uXryo5eVlFQqFROa/VqtpaWmpn99CodB/rFmn01GhUOi/dzTv5XKZG7AOlpeX+2UYvsrl8sDYTqfTr+sXL15ULpdToVBQuVxWrVZTLpczfXar1Tp0fHO5nBqNxrG4RqOhXC7Xj3G98ddoNA7t38WLF1Uulw/dQO71eiqXywPrGE6mVqtpcXFRKysr3LQ3WllZGXo+bWxsaGlpSSsrK9rZ2ZlyzhCHbrerlZUVraysaHl5eeQ549pmh32T8Dt+2PehaxwAAJg+LwiCYNaZAAAA82tlZeXQv8OL/mw2e+jv1Wp1bFqdTkcvvviidnd3Jal/k9Vl21kLbyhXq1WVSiVJT28cF4tF1ev1GefwqY2NDZXLZWWzWW1tbR16r1AoaGdnR6+99poymYyk/RuwpVJJGxsbKhaL6vV62tnZUT6fn0Huk63X6ymXy8n3/UN/v379+rG/lcvlfplWKpX+ORMOJtRqNUmStbve6/V08eJFSTpUH4+q1Woql8v9883i4sWL6vV6qtfrKhaLA2M2NjZUrVaP1bGodLtdLSws9OvpaRHWm3a7faxOTeo0lOXFixdVqVSGng/Ly8tqNBra2to69v2VZmk89p1OR7lc7tA5Erap7Xb70PF1bbNzuZyy2Wz/u7zVaml5eVmvvfbaRHEAAGA2GAwBAACRajQaymQyE90oX15eVqfTie3maZzCm8NHb6R5nifJfkM7TisrK6rVasdulC8vL2thYeHY4FM4yCNJS0tL+oM/+ANdv349VTfPolIul/XhD3946ABBKLzxWqlUtLa2NjCmVqtpZWVloroT3jDP5/NqNpsDY1ZWVrS0tDT080dZWlpSt9tVs9kceq63Wi1Vq9XYBgJzuZzq9fqpunEtPR3MjHIwJO1l2Wq1VCgUtLu7O7TdCtvF0zYYksZjH84Aarfbh/7ued6xHye4tNlhe3q0/hz9HNc4AAAwOzwmCwAAJMa8PqO70+mo1+spm80OvaEU9SNtNjY2Jn70RvjIovAmdrfb7T+eadAsHN/3lcvl1O129Qd/8Ae6du1aKgZCwl8EhzMwohCu5TCqLjcaDTUaDfm+P3IgolQqTXyz+/r165L2j/WwvNy5c2for+STLhw4xcmdhrKsVqvyfT8V7VaU5u3Yh4/3G/XdF7a/ly5dOvZePp9Xo9E4tL1Lm12r1ZTNZo/Vn2vXrqnT6fS3dY0DAACzw2AIAACYuVqtpuXlZXW7XXW73f6ztsOb9o1Go7/uQK1W08WLFw89nqvX62llZUXlclmFQuHYGgXhgsOFQqF/46NQKOjixYsqFAr9gYrwufHhOgiujg4uhMK1GuK4CZfNZlWpVMzb9Xo9dbtdZTIZZbPZftneunVr5ONjstmsCoWCrl+/rtu3b0c6gDArmUxGr732mra2tiJZ32VjY6P/iKtcLqeLFy8OXK8jrFtHHzE3yKBZFeE6I+G6NIPqaiaT6f/SeX19/dj7jUZD+Xw+1pvDCwsLh/4dnoe5XK5/LofltLy8fGzdkfCcDmewhHWu0Wj0byqG6wGE/x7XVriU3cbGRv+zjw6WxdWWWMomtLOz04/zPE/Ly8sD0x3VNo4qS5ftD+a70Wio1WoNzc+oYzpMq9Xqr+EQptfpdPrr4hxd7yZc++HozfJGo+F0vh3M58WLF4/VnzCtSetY3N9F4edfvHixXz4HyyJcy2djY2PosQ/zFJZ7o9FQr9fr59PSTi4tLalQKPTX7Di4NodlYCBc66rZbOq1114bOYslTHdQ2xZuF8a4tNnhjx0GDUyHf7t9+7Zz3EnrwDTqEAAAqRYAAAA4qNfrQalUCtbW1gLf94Otra2hcc1mc6LPyGazQTabPZZeNpsNJB36fN/3gyAIgna7HWQymaDdbve3qVargaSgUqkEQRAEW1tbQalUCiQF+Xw+WFtbC9rtdtBsNvt/K5VKQbPZDLa2toJisRhIOpTmKPl8PpAU1Ov1Q/mWNLKsTupoWbkIyyafzwfFYjHI5/PB7u6u07ZhuY3iUk9c65KLer0e5PP5IJPJBPV6PahWq4Hv+4GkoFgsmtKqVCpBNpsN1tbWnMvkoN3d3aBerwdra2v9OhvW29DW1lb/75OcJ+12O8jn8/1/h/Xs4GcM+qyj++P7/sTnaRAE/f0blUa73T50DMJz8OB5WK/X++fbwfpcLBaDtbW1/r+r1Wr/fA6C/boo6VDdcWkrxpVdmG4obCPC/YyrLbGUTaVS6cdVKpWg3W73tz9YRi5t47CydN3+4D4eLI9B+Rl3TIcJ83ewroXlfTC9MM2jfwuP87hz+uBxLZVKh9qSsN6ctI7F/V109HOPtoFbW1uH8jfs2AfB03M8fC+fz5vzcbRdCuvuoPZqkEna5LAsB7X/4f6G9c6lzQ7r/KA8t9vt/me5xp20DkyrDgEAkFYMhgAAgLFKpdKhGwuZTCaoVqv99w6+8vl8/2L84MvFoMGQIHh6A+XoTa4g2L+pe/DmzsG/H7yRE96MOJqXo3EHY11u1AVB0L+BEt4UkxRkMpkT3Wx2EQ4EWIQ3RjKZjGkfgyAYO2gxqp5YYqzCG0MHB2vCm16D6sw41Wo1yGazQalUmmhQJBTekDx4Ize8YTXq5my73Q4qlcqhV7i97/vHbmqFx3JQeuFA3cFy2NraCjKZzMT7FQSTDYYEwdNz+egxPzqgmMlkjh07lxv449qKcWUXDqyFdnd3j6UXV1viWjbD4g7etA/z49I2DitL1+3Den50H4/mZ9wxHSYc1DtalzKZzLHvi/Bm89G/ubSTYTtycFA7CI7X9ZPWsTi/iwbl+6BKpXJo/0YNhhy90T9JW3qwfQjPpUwmM7Jd3d3d7Q9OTLLf4ecM+t1nuL/DvnMGtdmjjndYN33fd44LgpPXgWnVIQAA0ojHZAEAgJEKhYLu3LmjW7duSVL/kRnhI6Gq1eqhV7PZVLPZPPb3kwgfd/HhD3/40N/D53wPeixF+NiSo5999NEZ4WMzDj7SJ/yby0Lu4SNjfN9Xu93uvyQde6RM1IrFYv/RNa5rkoT5ee211yTtPzbFddtRjyYZV09cYyZx//59SfvPZQ/rwrVr1yRponVVSqWStra2VCgU9OKLL2plZWWidIrFYn/xcss54Pu+SqWSyuVy//Em+Xy+X9/X19f7j5I7+CiisBwOCrc/+DiiSqXSX1NkVo7WpfB8Dcsrm81qY2Pj0CN5XBZ6H9dWjCu7er1+aJHj8O+DzpGo25Kj24SOlk3o6JoImUymX08naRsPiqJtlfYf5RWa9Jhms1n5vj/wkXPhoxWl/WMUrt10UKPRGPgIsWGO7vPR8j9pHQvFVX+O5vtgud2+fXvkQuEHhe1Qo9FQpVKZ6LGMB9v2l156SdJ++zOornS
      "text/plain": [
       "<Figure size 2000x600 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#plot residual energy against energyloss and try to find a good split (eg energyloss before and after the magnet)\n",
    "#upstream\n",
    "nbins=60\n",
    "\n",
    "fig, ((ax0, ax1)) = plt.subplots(nrows=1, ncols=2, figsize=(20,6))\n",
    "\n",
    "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",
    "ax0.set_ylim(0,5e4)\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(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",
    "ax1.set_ylim(0,5e4)\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, Upstream photons w/ brem_vtx_z$<9500$mm\")\n",
    "\n",
    "\"\"\"\n",
    "\"\"\"\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABjQAAAJPCAYAAADbtz2DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAACVoUlEQVR4nOz9f4wjZ37neX6iSq3qsmWJVWp4t3v3NCXmNMZqW3A7WGrPooHWrkWO1/BhF1PFLLUX8GhutsWcubwbYPquky7jztXl2XGaeV4sMIfEmFmH29YM1tsqsmpvDvDu2mS1Tw0Ig3MXafkKbnlnTKq679be6WlnsTWyq0uqrOf+yA4qf/BHfMmIYJB8vwBCquSXT3zjiYcRT8TDJ8JzzjkBAAAAAAAAAACk2IlZJwAAAAAAAAAAADAOAxoAAAAAAAAAACD1GNAAAAAAAAAAAACpx4AGAAAAAAAAAABIPQY0AAAAAAAAAABA6jGgAQAAAAAAAAAAUo8BDQAAAAAAAAAAkHoMaAAAgJF6vd6sUwAAAAAAAGBAAwAAjNZsNtVsNmedBgAAAAAAWHIMaAAAAAAAAAAAgNRjQAMAAAAwarfbs04BCKXX62lnZ0fdbnfWqSycbrernZ0dbs0IjgkAACSIAQ0AABCrdrutcrms1dVVraysaGtra9YpjdTr9bS1tSXP8+R5ns6cOaPV1VUVCoVU5r+zs6OVlZV+voVCoX+LsHa7rUKh0H/vaO7lcpmLnCGsrq726zB4lcvlgbHtdrvf1s+cOaNcLqdCoaByuaydnR3lcjnTspvN5qHtm8vlVK/Xj8XV63Xlcrl+TNiLa81m89j6rays9PMvFAra2trigu2c2tnZ0bPPPqu1tTW24QTW1taGfme3tra0srKitbU17e7uJpwZotbtdrW2tqa1tTWtrq6O/M6EPSYEfZ+gDzHseBs2DgAA7POcc27WSQAAgPRYW1s79O/gxDqbzR76e7VaHVtWu93WSy+9pHv37knaP2nv9XqhPjtrwUXharWqUqkkaf+i8erqqorFomq12owz/NDW1pbK5bKy2aw6nc6h9wqFgnZ3d3Xr1i1lMhlJ+xc5S6WStra2VCwW1ev1tLu7q3w+P4Ps063X6ymXy8n3/UN/v3z58rG/lcvlfp1WKpX+d6bX6/UHNCTJ2v3u9Xo6c+aMJB1qj0ft7OyoXC73v28WZ86cUa/XU61WU7FYPFRmpVJRt9vVxsaGKpWKuew06Xa7Onv2bP+7sAyCdtlqtY612WksQ12eOXNGlUpl6HdudXVV9XpdnU7n2DFyUS3idm+328rlcoe+I8E+u9VqHdq2YY8JuVxO2Wy231cIBo9v3bo1URwAADjAAQAAjFCr1Vyj0Zjos8Vi0WWz2YgzSkYmk3GSXKfTOfR3SS5tXahSqeQkuWq1eujvxWLRlUqlY/GtVsuVSiVXKpVcpVJxGxsb7t69ewllO182NjZcrVYbG1csFp0kV6lUhsZUq9WJ287GxoaT5PL5/NCYYHtOIpvNOklDv+u+7ztJA9vTPPF9/9h3etFVKhUnybVarUjLXfS6bDQaTtLIfWOw713kejhqEbe77/vO9/1jf5fkisXiob+FOSYE++ujbefocsLGAQCAw7jlFAAAiM283lO63W6r1+spm80O/dVt1Ldv2dramvg2E8EtpoIZFt1ut3+ro0GzYXzfVy6XU7fb1Z//+Z/r5ZdfXohf2wa/nA1mQkQhePbAqLZcr9dVr9fl+742NjaGxpVKpYl/cXv58mVJ+9t6WC7Xr18f+kvyaQUzfHZ2dub2e726ujq3uafNMtRltVqV7/sLsW+Myrxt9+BWfKOOrcH+/fz588fey+fzqtfrhz4f5piws7OjbDZ7rO28/PLLarfb/c+GjQMAAIcxoAEAACK3s7Oj1dVVdbtddbvd/r2hgwvv9Xq9/6yHnZ0dnTlz5tCtrnq9ntbW1lQul1UoFA49FyJ4f2dnR4VCoX9xoVAo9O/5Hww2BPc4P3PmzNBnHgxydIAgEDy7II6LXNlsdqLb+fR6PXW7XWUyGWWz2X7dXrt2beRtUrLZrAqFgi5fvqzXX3890kGAWclkMrp165Y6nU4kzzsJnh1RLpeVy+V05syZgc+vCNrW0du1DTLoVmXBczeC57QMaquZTKZ/K6jNzc1j79frdeXz+dguvmYymX57evXVVw+9N+772mw2+/eZX11dlbS/zsGzQY4+8yO4P3232+1/14NnhzSbzf6zQoKyBuWxtramlZWVfruu1+v9ZQT3yA/+PW5/FGb7bG1t9Zd9dFAtjv3VwXoJ8g7a6Orq6sAB193d3X7coPoLsy3H1WWYMqLapsNM096OrmeY7/TBPM+cOXOs/QRlTdrG4j7etdttnTlzpl83B+uhXC73n780bLsH+QR1Xq/X1ev1+jla9sMrKysqFAr9Z1gcfFaF5eJ+8GypRqOhW7dujbwdWFDuoH1n8LkgJswxIfhBxKDB6+Bvr7/+eug4afo2EHcbAgAgcbOeIgIAAGajVqu5UqnkNjY2Rt5CYppbTmWz2WO3nKrVav3b2xxcfnB7hVar5TKZzKHbowS36glup9PpdPq3+sjn825jY8O1Wq3+LULy+bwrlUqu0Wi4TqfTvx1Q2Fuu5PN5J+nQbSVqtZqTFOvtNia5PVdQN/l83hWLRZfP50PfPiqot1HCtJOwbSmMWq3m8vm8y2QyrlaruWq12r/d0dFbf4xTqVRcNpud+JZa9+7dc7VazW1sbPTbrI7cdqnT6fT/Psn3pNVqHbqNVNDOBt3a6eCyBt2iZNLvqXPjbzkV5Kojt1wL83117sNbqxwsP/i+bmxsHFpOtVrt/+3g9/fgdz34/h9cRrFYPFRWtVodmMPB9hlmfzRu+wTlHl2vYF3j2F8F5R0ss1ar9T97cF8S3HIqn8+7SqUytP7CbsthdRm2jCi36TCTtrdAsJ3H7TcObtdSqXRofxW0m2nbWNzHu4PLPLqP7XQ6h3Ibtt2d+3AfEryXz+fNtzk7ut8L2m7YW91Nss8P6nHQ8SVY36DNhTkmBO192O0eg2WFjXNu+jaQRBsCACBJDGgAALCESqXSoZP3TCbTf/5C8GyF4JXP5/snuwdfYQwa0HDuw4sURy8iObd/YXbQcwKCi0TBxZLghP9oLkfjDsaGfb5AcJEiuOgkyWUymakuGIcRXMy3CC48BM/8sDxDYdzAw6h2YomxCi68HBxwCS4sDWoz41SrVZfNZl2pVJrqWSHBRb+DF0qDC0KjLn62Wi1XqVQOvYLP+75/7KJRsC0HlRcMth2sh06n4zKZzMTr5Vy4AY179+4dW9ew39dgMOboRcNMJnNsH1EsFo8NOgxq2wcvGgdlHW0f4wY0gphR+6Nx2ycYgAsE9XSwvDj2V0HeR79vRwdkh8Udrb+w29K54XUZtoyotukw07S34G9h9sXBvuroMxWOfp+mbWNxHu+O5nxQpVI5tG6jBjSOXqyfZF99cP8TfI8ymczI/fa9e/f6AwyTPEPo4H7tqGB9hx3TBh0TRm3roF36vh86LjBtG0iiDQEAkBRuOQUAwJIpFAq6fv26rl27Jkn920MEt1eqVquHXo1GQ41G49jfpxHc2uGFF1449PfgvtSDbsEQ3J7j6LKP3iYiuEXE2bNnj/2t0+mMzS24NYrv+2q1Wv2XpIG3X4lSsVjs36Il7DM6gnxu3bolaf8WIWE/O+o2HOPaSdiYSdy+fVvS/n3Eg7bw8ssvS9JEzxkplUrqdDoqFAp66aWXtLa2NlE5xWJRjUZD0vF2OIrv+yqVSiqXy/3beOTz+X5739zc7N+W7eAtd4J6OCj4/MHb7lQqlf4zNpKyu7tr+r5ms1n5vj/wll3Bremk/VujBM+vOWrQLWF2d3f7/5/NZrW1tXXoFjejnmlytNxh+6Nx26dWq/X3EQf/Puh7GPX+6mB8IKj7oK0Gjj4jIJPJ9Ot9kn3vUVHsv6Votum07a1erw+8JdcwR9f56DaYto0F4mg/R3M+WGevv/56/1Z34wT7uXq9rkqlMtEtFA8eO4Jb21UqlYHtpNvtam1tTbl
      "text/plain": [
       "<Figure size 2000x600 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#downstream\n",
    "fig, ((ax0, ax1)) = plt.subplots(nrows=1, ncols=2, figsize=(20,6))\n",
    "\n",
    "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",
    "ax0.set_ylim(0,5e4)\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(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",
    "ax1.set_ylim(0,5e4)\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(a0[3],ax=ax1)\n",
    "fig.suptitle(r\"$e^\\pm$ from $B\\rightarrow K^\\ast ee$, $p>5$GeV, Downstream photons w/ brem_vtx_z$<9500$mm\")\n",
    "\n",
    "\"\"\"\n",
    "\"\"\"\n",
    "plt.show()"
   ]
  },
  {
   "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}