{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/hep/davide/miniconda3/envs/root_env/lib/ROOT.py:301: RuntimeWarning: numpy.dtype size changed, may indicate binary incompatibility\n",
      "  return _orig_ihook( name, *args, **kwds )\n"
     ]
    }
   ],
   "source": [
    "import ROOT as r\n",
    "import root_numpy as rn\n",
    "import pickle\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import array as array\n",
    "from xgboost import XGBClassifier\n",
    "from tools.data_processing import *\n",
    "from tools.select_and_fit_funct import select_and_fit, plot_MC_vs_data\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#b = r.TBrowser()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PATHs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "test=1\n",
    "l_index=0\n",
    "trig_index=5\n",
    "mother_index_fit=2\n",
    "\n",
    "\n",
    "l_flv=['e','mu']\n",
    "mother_ID=[\"Ds\",\"Dplus\",\"both\"]\n",
    "trig_cat=[\"Any\",\"L0l_TOS\",\"TISnotTOS\",\"L0G_TIS\",\"TISorTOS\",\"TOSnotTIS\",\"TOSandTIS\"]\n",
    "\n",
    "\n",
    "BDT_PATH=l_flv[l_index]+'/BDTs/test_'+str(test)\n",
    "FIT_PATH=l_flv[l_index]+'/fits'\n",
    "with open(BDT_PATH+'/variables_used.pickle', 'rb') as f:  \n",
    "        branches_needed_0 = pickle.load(f)\n",
    "        \n",
    "branches_needed_BDT=['cos_thetal']+branches_needed_0\n",
    "branches_old = ['cos_thetal'] + return_branches(data_index=1,mother_index=0,l_index=l_index,meson_index=0)\n",
    "branches=[]\n",
    "for label in branches_old:\n",
    "    if 'NDOF' not in label:\n",
    "        branches.append(label)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data loading, applying mass cuts, preprocessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "14942\n",
      "(13679, 1155, 1280)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAHUNJREFUeJzt3Xtw1eW97/H3R0BAtF6xegyK7daOWy4Ro+yDtUVr5yg4tloUFbtHZNeCtD3q9GIv1ku1tXugnk21k+GM1XOKXA7V1tktHGQXrTLV0giRq5eeKTapCilubRWQS77nj/UQF2El+SWsZK2VfF4zmfx+z+95fuu7niS/73p+lyeKCMzMzA4pdQBmZlYenBDMzAxwQjAzs8QJwczMACcEMzNLnBDMzAxwQjAzs8QJwczMACcEMzNL+pc6gM447rjjYvjw4aUOw8ysorzwwgt/jYihHdWrqIQwfPhw6urqSh2GmVlFkfRalno+ZWRmZoATgpmZJZkSgqSLJa2XtEnSbQW2S9IcSRslrZE0JpUPklQnqV7Sq5L+hySlbcdIWi5pnaQnJR1d3LdmZmad0eE1BEkDgVrgfOBN4DlJT0bE6rxqVwCnAGcCZwEPA6OB94FPRsR7kgYAK4ELgBXAXcDSiPiRpFvS+leK9s7MrOh2795NY2MjO3fuLHUoVsCgQYOoqqpiwIABXWqf5aLyWGBDRDQASFoETATyE8JEYF7k/rnCakn9JQ1Lbd5LdQYA/YCteW3GpuV5wPM4IZiVtcbGRo444giGDx9OGuxbmYgItm3bRmNjI6eeemqX9pHllFEV0JC33pjKMtWR1E9SPblE8HRErE91hkZEE0D6fnyhF5d0YzrtVNfU1JQhXDPrLjt37uTYY491MihDkjj22GMPavTW7ReVI2JvRFSTSxDnS7qgk+3nRkRNRNQMHdrhbbRm1s2cDMrXwf5ssiSERmBY3npVKutUnYh4G/g18E+pqEnSUID0fStmZlYyWa4hrAJGSKoCtgCTgemt6iwBrgMWpzuMmiOiQdJxwPsR8XdJg4FPA//aqs396fvSg343Ztajpj3yh6Lu76Hrz+mwjiSmTJnCvHnzANizZw8nnngiY8eO5Ve/+hUAS5cu5Y477uD9999n9+7dXHjhhTzwwAP77eeRRx5h6tSpLF++nIsuugiAX/7yl1x++eUsXryYSZMmMX78eGbNmkVNTQ0PPvggc+bMYeDAgezdu5d7772Xz372s1x//fVceumlLfXfeOMNhgwZwp49e/jEJz7BPffcw1FHHXXA+zj88MN59913DyjPf82e1mFCiIidkmYAy8iNKOZFRJ2k6Wl7LfAYcIGkjcAuYGpq/l+A/51uNR0ELIiIf0/b7gAWSbqBXKK5qojvy8pY/kEkywHALN+QIUNYv349O3bsYPDgwSxfvpyTTjqpZfvzzz/PzTffzLJlyxg+fDh79+5l7ty5Bfc1cuRIFi5c2JIQFixYwOjRow+o96c//YnZs2ezZs0ajjzySLZv305b1zQfffRRampq2LNnD9/73vf4zGc+w29/+9sivPPul+kaQkQsiYgzI+KMiLg3ldWmZEDkzIyIf4yI6oioS+Vr0/roiPhYRNyZt89tEXFRRIxM39/qhvdnZr3QhAkT+PWvfw3kDuLXXHNNy7b777+f7373u+yb96xfv37MmDGj4H7OP/98Vq1axe7du3n33Xf54x//SHV19QH1mpqa+NCHPsQRRxwBwGGHHcYpp5zSboz9+/fnrrvu4o033uDFF18sWOfWW29l9OjRjBs3ji1btuy3rbm5meuvv57vfOc7APzkJz/hox/9KOPGjeMLX/gCX/rSl9p9/a7wk8pmVnGuvvpqFi5cyM6dO1m7di1jx45t2bZ27VrOPvvsTPuRxEUXXcSyZct44oknuOyyywrWGzNmDEcffTQf+chHmDp1Ko8//ji5u+w7NmbMGF566aUDyt977z3OPfdcXnzxRSZOnMjtt9/esm3Pnj1MmTKF0047jXvuuYfXXnuNH/zgB6xevZpnn32Wl19+OdNrd5YTgplVnFGjRrF582YWLFjAhAkTDmpf+5LLwoUL9xtp5Ovfvz8rVqxg/vz5nH766Xz9619v+eTekbYSxyGHHMKkSZMAuOaaa1i5cmXLti9+8YuMGDGCb3/72wD8/ve/58ILL+TII4+kX79+Le2KzQnBSmraI3/Y78ssq8suu4yvfvWrBxzER44cyerVq9todaBzzz2XdevW8de//pXTTz+9zXqSGDduHN/85jdZuHAhjz32WKb919fXc8YZZ3RYL/+W0XHjxvHUU0/1+BPhTghmVpFuuOEG7rjjDkaOHLlf+c0338zdd9/Na6/lZnxubm6mtra23X3dd999fP/7329z++uvv8769etb1uvr66mqav187v727NnD3XffzQknnMCoUaMO2N7c3MwvfvELABYtWsR5553Xsm3atGlMmDCBq666ij179jB27Fieeuop/va3v7F3714ef/zxdl+7qyrq/yGYWXkp5V1iVVVVfOUrB852M27cOGbNmsWkSZPYtWsXe/fu5VOf+lS7+7rkkkva3b57925uuukmtm3bRkRw3HHHtXnn0pQpUzjssMNabjt94oknCtYbMmQIzz33HPfeey+DBw9uSQ773Hrrrbzzzjt8/vOf59FHH+VrX/sa1dXVnHDCCZx22mkMHjy43Zi7QlkvjJSDmpqa8D/IqXztnRrybaid17o/u7MPN23alOn0hxXf9u3bWxLN5ZdfzpQpU7j66qsPqFfoZyTphYjo8MEGnzIyM6sAt99+O2eddRannXYaJ554IldeeWXRX8OnjMzMKsDs2bO7/TU8QjAzM8AjBLNex1ODWFd5hGBmZoATgpmZJT5lZGZdN39ycfd37aIOqxRr+ms7kBOCmVWUYk5/bfvzKSMzqzjFmv7a9ueEYGYVp1jTX9v+nBDMrOIUc/pr+4ATgplVpGJNf20fcEIws4pUzOmvLcd3GZlZ12W4TbS7FHP6a8txQjCzivLuu+8eUDZ+/HjGjx/fsn7ppZdy6aWX9mBUvYNPGZmZGeCEYGZmiROCmZkBTghmZpY4IZiZGZAxIUi6WNJ6SZsk3VZguyTNkbRR0hpJY1L5MEnPpLavSPpGXps7Jf1FUn368uOGZmYl1OFtp5IGArXA+cCbwHOSnoyI/EcBrwBOAc4EzgIeBkYDu4EvRcRaSUcAqyUti4j61O7+iJhVvLdjZj1p8SuLi7q/K0/v+B/H9+vXj5EjRyKpZSrsW265hUMOafvz7ebNm/nd737HtddeW8xwe50sI4SxwIaIaIiI3cAiYGKrOhOBeZGzGugvaVhEvBkRawEi4u/AWuAkzMy6aPDgwdTX17NmzRp+85vfsGLFCu66665222zevJn58+f3UISVK0tCqAIa8tYbU1mn6kgaDpwDrMwrninpJUmPSjq20ItLulFSnaS6pqamDOGaWV9x9NFHM3fuXB544AEigs2bN/Pxj3+c6upqzjzzTJ5++mkAbrvtNp599lmqq6u5//7726zX1/XIk8qSDgd+DtwcEe+k4geB7wEB3AnMAaa0bhsRc4G5ADU1NdET8ZpZ5TjppJMYMGAAW7du5cMf/jArVqzg0EMP5dVXX+WKK65g3bp13HfffcyaNavlP6rt2LGjYL2+LktCaASG5a1XpbJCdZ5vXUfSAOAxYEFEPL6vQUS0fNyXVAs83cnYzcwAiMh9Vty+fTs33XQT69at49BDD+Xll18uWD9rvb4myymjVcAISVXp4D4ZWNqqzhLSp/t0h1FzRDRIEvAQsCkiZuc3kHR83urngI1dfA9m1oe9/vrr7N27l+OPP57Zs2dz8skns2HDBurq6mhubi7YJmu9vqbDEUJE7JQ0A1hGLoHMi4g6SdPT9lpyI4ALJG0EdgFTU/PzgM8D6yTtu7PoWxGxBPiRpFHAocCfgWlFfF9m1ge8/fbbTJ8+nZkzZyKJHTt2UFVVhSQWLFjA3r17gdyF6O3bt7e0a6teX5fpGkI6gC9pVVabtxzAzALtVgJqY5/XdSpSMys7WW4TLbYdO3ZQXV3dctvptddey6233grAjBkzuOyyy3j00Uf59Kc/zZAhQwCorq5m165djBgxgmnTprVZr6/z9NdmVlHa+zR/+umn89JLL7Ws//CHPwRg4MCBrFy5cr+6her1dZ66wszMACcEMzNLnBDMrFP23eJp5edgfzZOCGaW2aBBg9i2bZuTQhmKCLZt28agQYO6vA9fVDazzKqqqmhsbMTTyJSnQYMGUVXVemah7JwQzCyzAQMGcOqpp5Y6DOsmPmVkZmaAE4KZmSVOCGZmBjghmJlZ4oRgZmaAE4KZmSVOCGZmBjghmJlZ4oRgZmaAE4KZmSVOCGZmBjghmJlZ4oRgZmaAE4KZmSWe/tp6xLRH/lDqEHoV96d1B48QzMwMcEIwM7PECcHMzAAnBDMzSzIlBEkXS1ovaZOk2wpsl6Q5kjZKWiNpTCofJumZ1PYVSd/Ia3OMpOWS1kl6UtLRxXtbZmbWWR0mBEkDgVrgEmAUMGnfAT/PFcApwJnANODhVL4b+FJEjADOBv5FUnXadhewNCJGAkvTupmZlUiWEcJYYENENETEbmARMLFVnYnAvMhZDfSXNCwi3oyItQAR8XdgLXBSXpufpeV5BfZpZmY9KEtCqAIa8tYbU1mn6kgaDpwDrExFQyOiCSB9P77Qi0u6UVKdpLqmpqYM4ZqZWVf0yEVlSYcDPwdujoh3OtM2IuZGRE1E1AwdOrR7AjQzs0xPKjcCw/LWq1JZoTrPt64jaQDwGLAgIh7Pa9MkaWhENEkaCmztQvzWy+Q/gfvQ9eeUMBKzvifLCGEVMEJSVTq4TyZ3ETjfEmAKQLrg3BwRDZIEPARsiojZBdpcl5avK7BPMzPrQR2OECJip6QZwDJyCWReRNRJmp6215IbAVwgaSOwC5iamp8HfB5YJ6k+lX0rIpYAdwCLJN0AbAGuKuL7MjOzTso0uV06gC9pVVabtxzAzALtVgJqY5/bgIs6E6yZmXUfP6lsZmaAp78269V8kd46wyMEMzMDnBDMzCxxQjAzM8AJwczMEicEMzMDnBDMzCxxQjAzM8AJwczMEicEMzMDnBDMzCxxQjAzM8AJwczMEicEMzMDnBDMzCzx9NdmvcSXt3xnv/Uff/ieEkVilcojBDMzA5wQzMws8SkjKymf5jArHx4hmJkZ4BGC9bDWIwIzKx8eIZiZGeCEYGZmiU8ZmVUgn3qz7uARgpmZARkTgqSLJa2XtEnSbQW2S9IcSRslrZE0Jm/bTyVtlbS+VZs7Jf1FUn36mnDwb8fMzLqqw4QgaSBQC1wCjAIm5R/wkyuAU4AzgWnAw3nbHgEubmP390dEdfpa0snYzcysiLKMEMYCGyKiISJ2A4uAia3qTATmRc5qoL+kYQAR8QzwVjGDNjOz4suSEKqAhrz1xlTW2TqFzJT0kqRHJR1bqIKkGyXVSapramrKsEszM+uKUl5UfhD4B+Afgf8HzClUKSLmRkRNRNQMHTq0J+MzM+tTsiSERmBY3npVKutsnf1ERFNE7I2IZnLXKM7JEIuZmXWTLAlhFTBCUpWkAcBkYGmrOkuAKQDpgnNzRDTQDknH561+DtiYOWozMyu6Dh9Mi4idkmYAy8glkHkRUSdpetpeCzwGXCBpI7ALmLqvvaQFwHjgOEmNwB0R8RDwI0mjgEOBP5O7O8nMzEok05PK6ZbQJa3KavOWA5jZRttr2ii/LnuYZmbW3fykspmZAU4IZmaWOCGYmRnghGBmZokTgpmZAU4IZmaWOCGYmRnghGBmZokTgpmZAU4IZmaWOCGYmRnghGBmZokTgpmZAU4IZmaWOCGYmRnghGBmZokTgpmZARn/Y5qZVZ4vb/lOq5JlJYnDKodHCGZmBniEYFYxpj3yh5blL5cwDuu9nBCse82fDMCXt7xd4kDMrCM+ZWRmZoATgpmZJT5lZN1m2iN/8KkiswrihGBWCeZPdnK1budTRmZmBmRMCJIulrRe0iZJtxXYLklzJG2UtEbSmLxtP5W0VdL6Vm2OkbRc0jpJT0o6+uDfjpmZdVWHCUHSQKAWuAQYBUzKP+AnVwCnAGcC04CH87Y9AlxcYNd3AUsjYiSwNK2bmVmJZBkhjAU2RERDROwGFgETW9WZCMyLnNVAf0nDACLiGeCtAvudCPwsLc8rsE8zM+tBWRJCFdCQt96Yyjpbp7WhEdEEkL4fnyEWMzPrJmV/UVnSjZLqJNU1NTWVOhwzs14ry22njcCwvPWqVFaozvPt1GmtSdLQiGiSNBTYWqhSRMwF5gLU1NREhnjNrID8uZAeuv6cEkZi5SpLQlgFjJBUBWwBJgPTW9VZAlwHLE4XnJsjooH27Wtzf/q+tDOBW++UP2Vz/Q/hxx++B/ABzKwndHjKKCJ2AjPITaa+FvhFRNRJmi5pX2J4DPiLpI3AT4Gp+9pLWgA8B3xMUqOkaWnTHcBESevIXVD+brHelJmZdV6mJ5UjYgm5T/T5ZbV5ywHMbKPtNW2UbwMuyhyplb80s+k+frLWrLKU/UVlMzPrGU4IZmYGOCGYmVnihGBmZoATgpmZJU4IZmYGOCGYmVni/5hmVo5aPdNRDPlPgTP/KLh2UdFfwyqbE4JZhVsx+L0Dyi7cMaQEkVil8ykjMzMDnBDMzCxxQjAzM8DXEMwqTqFrBmbF4BGCmZkBTghmZpY4IZiZGeBrCGa9UuvrDH4uwbLwCMHMzACPEKzMtUy3MP+o3HdPt2DWbTxCMDMzwCMEs7K3uPk/AWgYvL3EkVhv5xGCmZkBTghmZpY4IZiZGeCEYGZmiROCmZkBGROCpIslrZe0SdJtBbZL0hxJGyWtkTSmo7aS7pT0F0n16WtCcd6SmZl1RYe3nUoaCNQC5wNvAs9JejIiVudVuwI4BTgTOAt4GBidoe39ETGraO/GzMy6LMtzCGOBDRHRACBpETARyE8IE4F5ERHAakn9JQ0DTs3Q1voQz7HTscWvLIb07IFZT8qSEKqAhrz1RmB8hjpVGdrOlPQvwAvAVyJiW6aordfwP4g3Kx+lfFL5QeB7QAB3AnOAKa0rSboRuBHg5JNP7sHwrEPzJ5c6AjMroiwJoREYlrdelcoK1Xm+VZ0BbbWNiKZ9hZJqgacLvXhEzAXmAtTU1ESGeM16jYa3ijNdReuRWDVHFWW/1rtkSQirgBGSqoAtwGRgeqs6S4DrgMXpDqPmiGiQ1NRWW0nHR8TW1P5zwMaDfjdWdvz/f80qR4cJISJ2SpoBLCN3m+q8iKiTND1trwUeAy6QtBHYBUxtr23a9Y8kjQIOBf4MTCvuW7NSWJx3MdSTsXVs8SuLSx2CWYtM1xAiYgm5UUB+WW3ecgAzs7ZN5dd1KlKz3mzV/yx1BGae/tqKp77hbY8KzCqYE4Jl57uKeo3Fzf8JrU5XXXn6lSWKxsqFE4J12eJWD08Va3TQmx5ea32NoFwOug1vbeep321uWf/nccNLFouVD09uZ2ZmgEcIVgFWDH6PV5vTIyivLC6bT9mdki4aL9538fjcL5QwmJwL/vbEByurDqP1/U4V2c92UJwQrCK0PKD1f/8N6n7+wYZrF5UmIArfMuqDqFUyJwSrOPtdu0gH5e44EJfr+X+z7uKEYL1CRw949dTBPPODZn7uwMqQE4K1rdVtpq3vKjKz3sUJwfqkrkwZ4WkmrLdzQrAWBxzwetGIwAdzs445IZh1hz5wjcB3WfU+Tgh9hO+YsQ61SmIeU/U9Tgh9VK85hVLok3gZPPTVV5XL3V7WNU4IFaArn+57zQG/HLVOQk5AReXRbOk4IZShjg7mPthbX+IE0XOcEOwDveVCaE9/gu8t/dZaD/SjP9yUFyeEg5DlLossv/D+xGN9RTESgEcM3ccJoQyU7FNSb/1kW8FaJvGzzHz7a/E4IVjv54vAfU5nRxFOKjl9NiF01y9AXzsn6k+0vUP+z3HYMYeVMJLu0VN/l5V+222fTQiF9LWDuSUdnTrziCKnjz/z0ReuXfSZhOCDfe/RelTSGz/RWun1xO3f5Xaqqs8kBLMWnb2Y7ovvViTl/sHUCaG36uPDe+shnU2W/h0sa4eUOgAzMysPmUYIki4GZgH9gP8VEfe12i7g34CLgPeBaRGxur22ko4BFgEnAG8AkyOi90zA39N8WqNi+U6tPL7AX1IdJgRJA4Fa4HzgTeA5SU/uO+AnVwCnAGcCZwEPA6M7aHsXsDQifiTplrT+leK9NTtAkZJGuR3Aevstkz3N/dl3ZRkhjAU2REQDgKRFwEQgPyFMBOZFRACrJfWXNAw4tZ22E9O+AeYBz+OEkJ1HBBWt3JJqj6mEhwQrIcZukiUhVAENeeuNwPgMdao6aDs0IpoAIqJJ0vGZo+5pB/sLkuXg3XqfJf6lrNQDVjl/uq3EPm0v5qL0b3fc8dXZv6Xu+HB1sH+/8ycfWHbtoq7Hk1HZ32Uk6UbgxrT6rqSXu7ir44C/Fieq5cXZTdv7LBBrd7zmQStin3Y7x1p8ZRpnwb+VvFg7+7dUjL+9Tu2jcL9O+T8HE8ApWSplSQiNwLC89apUVqjO863qDGinbZOkoWl0MBTYWujFI2IuMDdDnO2SVBcRNQe7n55QKbFWSpzgWLtDpcQJjjWrLLedrgJGSKqSNACYDCxtVWcJMAVA0higOV03aK/tEuC6tHxdgX2amVkP6nCEEBE7Jc0AlpFLIPMiok7S9LS9FngMuEDSRmAXMLW9tmnXdwCLJN0AbAGuKu5bMzOzzsh0DSEilpD7RJ9fVpu3HMDMrG1T+TZyzy30lIM+7dSDKiXWSokTHGt3qJQ4wbFmotyx3MzM+jpPXWFmZkCFJwRJP5W0VdL6vLJFkurT12ZJ9al8uKQdedtq89qcLWmNpI2S5qSpOLo7zvMkvShpg6S1ks7L2/ZNSZskrZf033oqzs7GWoZ9WiNpdYrz3yV9KG9bufVpwVhL3KfDJD2T+ugVSd9I5cdIWi5pnaQnJR2d16Yk/drZWMu0X69MP/9mSTWt2pTm9zUiKvYL+AQwBljfxvbZwHfT8vB26q0Fzk7LTwBXdHecwErgkrQ8AViZls8G6sjdslsFbAYG9kScXYi13Pp0HfDJtHwDMLuM+7StWEvZpycAo9LyEcCrQDXwY+DWVH4LMKfU/dqFWMuxX88APgY8DdTk1S9Zv1b0CCEingHeKrQtZc6rgAXt7UPSyUC/iHghFc0jN61Gd8fZCOz7BHsk8Oe0PBFYFBG7I6IR2ACc2xNxdiHWgkrYp/8APJOWlwOXpeVy7NO2Yi2oh/r0zYhYm5b/Tu7gc1J6nZ8VeN2S9WsXYi2olLFGxKaIKPSgbcn6taITQgfOB7ZExKt5ZcPTqY/nJH0qlbU17UZ3+wYwW1IDudlgv9lBPKWKE9qOFcqrTzcBn0nLVwIndxBPKfu0rVihDPpU0nDgHHKjw/2mmQH2TTNTFv2aMVYov35tS8n6tTcnhGvYf3TwBlAVEaPJ3SL7s/xzoSXwEPDfI2IYuaHtQyWMpSNtxVpuffrPwM3pXP1x5KZiL1dtxVryPpV0OPBz4OaIeKcnX7uzOhGr+zWDsp/LqCsk9Sc3JffZ+8oi4n3SH11ErE5/iGeQbWqO7vBfgU+n5cXkpgynnXhKFSe0EWu59WlErCdNnpg+iU1Im8quT9uKtdR9qtyMAo8BCyLi8VTc1jQzJe3XzsRapv3alpL1a28dIVwEvJTOvwEg6VhJh6Tl4cAI4I8R8WegWbkpNyA3BUdPTKPxGvDJtHwhuQtHkHuIb7KkAZKqUpyrShhnm7GWW59KOi59F/AtPhjJlF2fthVrKfs0xfIQsCkiZudtamuamZL1a2djLdN+bUvpfl+LeYW6p7/InRJ6A9hNLlNOS+WPANNb1Z1E7uLMOmA9cGXethqgHtgIPEB6YK874wTOA15Mr1kP/FNe/W+TO8e8gXR3T0/E2dlYy7BPbwZeSrHcl/+aZdinBWMtcZ9+HAhyFz3r09cE4FjgP1JM/wEcU+p+7WysZdqvl6ffh/fJTd+zrNT96ieVzcwM6L2njMzMrJOcEMzMDHBCMDOzxAnBzMwAJwQzM0ucEMzMDHBCMDOzxAnBzMwA+P9ZqbNp5N+XVAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#MC_Dplus_sig_dict, MC_Ds_sig_dict, _ = load_datasets(l_index)\n",
    "\"\"\"\n",
    "Data\n",
    "\n",
    "\"\"\"\n",
    "if l_flv[l_index]=='e' and trig_cat[trig_index]!='Any':\n",
    "    with open('/disk/lhcb_data/davide/Rphipi/data_for_BDT_selection/'+l_flv[l_index]+l_flv[l_index]+'/TrigCats/data_for_BDT_selection_'+l_flv[l_index]+l_flv[l_index]+'_'+trig_cat[trig_index]+'.pickle', 'rb') as f:\n",
    "        data_dict=pickle.load(f)\n",
    "        \n",
    "#    with open('/disk/lhcb_data/davide/Rphipi/MC/for_fit/TrigCats/'+l_flv[l_index]+l_flv[l_index]+'/'+'MC_Dplus_'+l_flv[l_index]+l_flv[l_index]+'_'+trig_cat[trig_index]+'.pickle', 'rb') as f:\n",
    "#        MC_Dplus_tuple_dict=pickle.load(f)\n",
    "#        \n",
    "#    with open('/disk/lhcb_data/davide/Rphipi/MC/for_fit/TrigCats/'+l_flv[l_index]+l_flv[l_index]+'/'+'MC_Ds_'+l_flv[l_index]+l_flv[l_index]+'_'+trig_cat[trig_index]+'.pickle', 'rb') as f:\n",
    "#        MC_Ds_tuple_dict=pickle.load(f)\n",
    "    \n",
    "elif l_flv[l_index]=='e' and trig_cat[trig_index]=='Any':\n",
    "    with open('/disk/lhcb_data/davide/Rphipi/data_for_BDT_selection/'+l_flv[l_index]+l_flv[l_index]+'/'+'data_for_BDT_selection_'+l_flv[l_index]+l_flv[l_index]+'.pickle', 'rb') as f:\n",
    "        data_dict=pickle.load(f)\n",
    "        \n",
    "if l_flv[l_index]=='mu' and trig_cat[trig_index]!='Any':\n",
    "    with open('/disk/lhcb_data/davide/Rphipi/data_for_BDT_selection/'+l_flv[l_index]+l_flv[l_index]+'/TrigCats/data_for_BDT_selection_'+l_flv[l_index]+l_flv[l_index]+'_'+trig_cat[trig_index]+'.pickle', 'rb') as f:\n",
    "        data_dict=pickle.load(f)\n",
    "        \n",
    "#    with open('/disk/lhcb_data/davide/Rphipi/MC/for_fit/TrigCats/'+l_flv[l_index]+l_flv[l_index]+'/'+'MC_Dplus_'+l_flv[l_index]+l_flv[l_index]+'_'+trig_cat[trig_index]+'.pickle', 'rb') as f:\n",
    "#        MC_Dplus_tuple_dict=pickle.load(f)\n",
    "#        \n",
    "#    with open('/disk/lhcb_data/davide/Rphipi/MC/for_fit/TrigCats/'+l_flv[l_index]+l_flv[l_index]+'/'+'MC_Ds_'+l_flv[l_index]+l_flv[l_index]+'_'+trig_cat[trig_index]+'.pickle', 'rb') as f:\n",
    "#        MC_Ds_tuple_dict=pickle.load(f)\n",
    "    \n",
    "elif l_flv[l_index]=='mu' and trig_cat[trig_index]=='Any':\n",
    "    with open('/disk/lhcb_data/davide/Rphipi/data_for_BDT_selection/'+l_flv[l_index]+l_flv[l_index]+'/'+'data_for_BDT_selection_'+l_flv[l_index]+l_flv[l_index]+'.pickle', 'rb') as f:\n",
    "        data_dict=pickle.load(f)\n",
    "\n",
    "print(data_dict[\"Ds_ConsD_M\"].shape[0])\n",
    "\"\"\"\n",
    "MC signal\n",
    "\n",
    "\"\"\"    \n",
    "\n",
    "with open('/disk/lhcb_data/davide/Rphipi/MC/for_fit/'+l_flv[l_index]+l_flv[l_index]+'/MC_'+mother_ID[1]+'_'+l_flv[l_index]+l_flv[l_index]+'.pickle', 'rb') as f:\n",
    "    MC_Dplus_tuple_dict=pickle.load(f)\n",
    "with open('/disk/lhcb_data/davide/Rphipi/MC/for_fit/'+l_flv[l_index]+l_flv[l_index]+'/MC_'+mother_ID[0]+'_'+l_flv[l_index]+l_flv[l_index]+'.pickle', 'rb') as f:\n",
    "    MC_Ds_tuple_dict=pickle.load(f)\n",
    "\n",
    "\"\"\"\n",
    "MC MISID bkg\n",
    "\n",
    "\"\"\" \n",
    "\n",
    "mc_Ds_MISID_mass=rn.root2array(\n",
    "        \n",
    "        filenames=\"/disk/lhcb_data/davide/Rphipi/MC/mc_MISID/\"+l_flv[l_index]+l_flv[l_index]+\"/\"+mother_ID[0]+\"23pi_\"+l_flv[l_index]+\"MISID.root\",\n",
    "        treename = 'DecayTree',\n",
    "        branches = 'Dsp_0_M',\n",
    "    )*1000\n",
    "\n",
    "mc_Dplus_MISID_mass=rn.root2array(\n",
    "        \n",
    "        filenames=\"/disk/lhcb_data/davide/Rphipi/MC/mc_MISID/\"+l_flv[l_index]+l_flv[l_index]+\"/\"+mother_ID[1]+\"23pi_\"+l_flv[l_index]+\"MISID.root\",\n",
    "        treename = 'DecayTree',\n",
    "        branches = 'Dp_0_M',\n",
    "    )*1000\n",
    "\n",
    "\"\"\"\n",
    "Normalising the chi2s in data and MC\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "MC_Dplus_tuple_dict, MC_Ds_tuple_dict, data_dict = norm_chi2(MC_Dplus_tuple_dict, MC_Ds_tuple_dict, data_dict)\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Applying mass cuts on data, \n",
    "retrieving mass distributions for sig MC and data\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "if l_index==1:\n",
    "    if mother_index_fit==0:\n",
    "        #fitting only the Ds\n",
    "        lower_cut=1930\n",
    "        upper_cut=2010\n",
    "    if mother_index_fit==1:\n",
    "        #fitting only the Dplus\n",
    "        lower_cut=1860\n",
    "        upper_cut=1880\n",
    "    if mother_index_fit==2:\n",
    "        lower_cut=1820\n",
    "        upper_cut=2040\n",
    "\n",
    "if l_index==0:\n",
    "    if mother_index_fit==0:\n",
    "        #fitting only the Ds\n",
    "        lower_cut=1900\n",
    "        upper_cut=2040\n",
    "    if mother_index_fit==1:\n",
    "        #fitting only the Dplus\n",
    "        lower_cut=1820\n",
    "        upper_cut=1920\n",
    "    if mother_index_fit==2:\n",
    "        lower_cut=1750\n",
    "        upper_cut=2100\n",
    "        \n",
    "        \n",
    "data_mass, mc_Dplus_mass, mc_Ds_mass, data_dict = mass_cut_for_fit(lower_cut, upper_cut, mother_index_fit=mother_index_fit, l_index=l_index,\n",
    "                                                 branches_needed=branches_needed_BDT, \n",
    "                                                 data_dict=data_dict, MC_Dplus_dict=MC_Dplus_tuple_dict, MC_Ds_dict=MC_Ds_tuple_dict)\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Preprocessing\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "\n",
    "\n",
    "data_to_select, mc_Dplus_to_select, mc_Ds_to_select = preprocess_for_XGBoost(MC_Dplus_tuple_dict, MC_Ds_tuple_dict, data_dict, branches_needed_BDT, mother_index_fit=mother_index_fit)\n",
    "\n",
    "\n",
    "plt.hist(np.concatenate((mc_Dplus_MISID_mass, mc_Ds_MISID_mass)),bins=70,density=True,alpha=0.7,label=\"MC MISID bkg\");\n",
    "plt.hist(np.concatenate((mc_Dplus_mass, mc_Ds_mass)),bins=70,density=True,alpha=0.7,label=\"MC\");\n",
    "plt.hist(data_mass,bins=70,density=True,alpha=0.4,label=\"Data\");\n",
    "plt.legend();\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Just RooFit it!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BDT cut at 0.0615789473684 fit\n",
      "BDT selection..\n",
      "(499, 474)\n",
      "BDT Signal selection efficiency: 0.399589322382\n",
      "chi2 0.747467122077\n",
      "Saving fit results...\n"
     ]
    }
   ],
   "source": [
    "#Just fit it!\n",
    "\n",
    "#if __name__=='__main__':\n",
    "\n",
    "  # parser = argparse.ArgumentParser(description='Fit at different BDT cuts')\n",
    "  # parser.add_argument('-x_cut', metavar='BDT cut', type =float, help='BDT cut')\n",
    "  # args = parser.parse_args()\n",
    "    \n",
    "x_cut_values = np.linspace(0.01,0.99,num=20)\n",
    "data={}\n",
    "MC_Dplus={}\n",
    "MC_Ds={}\n",
    "fit_results={}\n",
    "#Just fit it!\n",
    "for i, x_cut in enumerate(x_cut_values[1:2]):\n",
    "    \n",
    "    if not os.path.exists(FIT_PATH+'/'+mother_ID[mother_index_fit]+'/{0}/'.format(i)):\n",
    "        os.mkdir(FIT_PATH+'/'+mother_ID[mother_index_fit]+'/{0}/'.format(i))\n",
    "    \n",
    "    fit_results[i], data[i], MC_Dplus[i], MC_Ds[i] = select_and_fit(BDT_PATH, mother_index_fit, l_index, trig_index, lower_cut, upper_cut,\n",
    "                                               data_mass, mc_Dplus_mass, mc_Ds_mass, mc_Dplus_MISID_mass, mc_Ds_MISID_mass,\n",
    "                                               data_to_select, mc_Dplus_to_select, mc_Ds_to_select,\n",
    "                                               data_dict, MC_Dplus_tuple_dict, MC_Ds_tuple_dict, branches_needed_BDT, branches, \n",
    "                                               x_cut, i, test, plots=True)\n",
    "    if l_flv[l_index]=='mu':\n",
    "        with open(FIT_PATH+'/'+mother_ID[mother_index_fit]+'/{0}/fit_results_{1}.txt'.format(i, i), 'w') as f:\n",
    "            print('Saving fit results...')\n",
    "            for key, value in fit_results[i].items():\n",
    "                f.write('%s:%s\\n' % (key, value))\n",
    "                \n",
    "    elif l_flv[l_index]=='e' and trig_cat[trig_index]=='Any':\n",
    "        with open(FIT_PATH+'/'+mother_ID[mother_index_fit]+'/{0}/fit_results_{1}.txt'.format(i, i), 'w') as f:\n",
    "            print('Saving fit results...')\n",
    "            for key, value in fit_results[i].items():\n",
    "                f.write('%s:%s\\n' % (key, value))\n",
    "    elif l_flv[l_index]=='e' and trig_cat[trig_index]!='Any':\n",
    "        with open(FIT_PATH+'/'+mother_ID[mother_index_fit]+'/{0}'.format(i)+'/TrigCats/'+trig_cat[trig_index]+'/fit_results_{0}.txt'.format(i), 'w') as f:\n",
    "            print('Saving fit results...')\n",
    "            for key, value in fit_results[i].items():\n",
    "                f.write('%s:%s\\n' % (key, value))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<BarContainer object of 9 artists>"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAD/RJREFUeJzt3X2MnWlZx/HvD1obFzUqHTQ6s50NKGpLgTpYFUmsJNh1RMK6Wk1XtC423ayvMYECGiTRZDErmArJpL7sH8xGmqVgDJ3NLgl/IBFsZrvL0k5BQ1KZIS7MhoRgdNOFXv5xbpOT6XTOMy+dM8t8P8lmz3M/93XOdZ/snt885znnOakqJEl63rAbkCRtDQaCJAkwECRJjYEgSQIMBElSYyBIkgADQZLUGAiSJMBAkCQ1O4bdwGrs3r27xsfHh92GJD2nPPbYY09X1cigec+pQBgfH2d2dnbYbUjSc0qS/+wyz7eMJElAx0BIcjjJxSSXk5xcZn+SnEoyl+TxJAe61Cb5vSRPtv33r385kqS1GviWUZJdwBTwGuAp4FNJHq2qC33T7gD2AHuBVwIPAC9fqTbJJHAYmKiqq0l2b+TCJEmr0+UI4SBwqarmq+pZ4AwwuWTOJDBdPReAHUnGBtT+DvCXVXUVoKqe3oD1SJLWqEsgjALzfdsLbazLnJVqfwR4XZLPJPl0klcv9+BJjieZTTK7uLjYoV1J0loM86Ty84DvAl4B/D7wwSTPXzqpqk5X1URVTYyMDPzUlCRpjboEwgIw1rc92sa6zFmpdh74cHub6TxwFfi+7q1LkjZSl0A4D+xLMppkJ3AEeHjJnBngKED7hNG1qpofUHsOONRqfhi4BfjKOtcjSVqjgZ8yqqpnktwDPEIvQKarajbJibZ/CjgLHEoyR+8v/WMr1ba7fh/wD0kute1jVfWNDVybJGkVUlXD7qGziYmJWus3lcdPntvgblZ25b6lH8SSpOFI8lhVTQya5zeVJUmAgSBJagwESRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSpMRAkSYCBIElqDARJEmAgSJIaA0GSBBgIkqTGQJAkAQaCJKkxECRJgIEgSWoMBEkSYCBIkhoDQZIEGAiSpMZAkCQBBoIkqTEQJEmAgSBJajoFQpLDSS4muZzk5DL7k+RUkrkkjyc5MKg2yZ8l+VKSJ9o/v7AxS5IkrcWOQROS7AKmgNcATwGfSvJoVV3om3YHsAfYC7wSeAB4eYfa91bV/Ru2GknSmnU5QjgIXKqq+ap6FjgDTC6ZMwlMV88FYEeSsY61kqQtoEsgjALzfdsLbazLnEG19yb5XJIHk7ywc9eSpA03zJPK7wdeAvwY8AXg1HKTkhxPMptkdnFxcTP7k6RtpUsgLABjfdujbazLnBvWVtViVX2zqq7RO8/wquUevKpOV9VEVU2MjIx0aFeStBZdAuE8sC/JaJKdwBHg4SVzZoCjAO0TRteqan6l2iQv6qv/ZWBuXSuRJK3LwE8ZVdUzSe4BHqEXINNVNZvkRNs/BZwFDiWZA64Cx1aqbXf9niT7gW8DvgjcvbFLkyStxsBAAKiqGXpHAf1jU323C7i3a20bv2tVnUqSbiq/qSxJAgwESVJjIEiSAANBktQYCJIkwECQJDUGgiQJMBAkSY2BIEkCDARJUmMgSJKAjtcy0sYaP3lu0x7ryn3+QJ2kbjxCkCQBBoIkqTEQJEmAgSBJagwESRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSpMRAkSYCBIElqDARJEmAgSJKaToGQ5HCSi0kuJzm5zP4kOZVkLsnjSQ6sovaPk1SS3etbiiRpPQYGQpJdwBRwO7AfuLP/Bb+5A9gD7AXuBh7oUptkDHgd8MV1r0SStC5djhAOApeqar6qngXOAEt/dWUSmK6eC8CO9mI/qPa9wFuAWu9CJEnr0yUQRoH5vu2FNtZlzg1rk7wB+FJVfWaVPUuSboKh/IRmkluAt9N7u2jQ3OPAcYBbb731JncmSdtXlyOEBWCsb3u0jXWZc6PxFwO3AZ9JcqWNX0jy/UsfvKpOV9VEVU2MjIx0aFeStBZdAuE8sC/JaJKdwBHg4SVzZoCjAO2k8bWqmr9RbVV9tqpeVFXjVTVOLyQOVNVTG7MsSdJqDXzLqKqeSXIP8Ai9AJmuqtkkJ9r+KeAscCjJHHAVOLZS7c1ZiiRpPTqdQ6iqGXpHAf1jU323C7i3a+0yc8a79CFJunmGclJZW8P4yXOb+nhX7lv6aWVJW4mXrpAkAQaCJKkxECRJgIEgSWoMBEkSYCBIkhoDQZIEGAiSpMZAkCQBBoIkqTEQJEmAgSBJagwESRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSpMRAkSYCBIElqDARJEmAgSJIaA0GSBBgIkqTGQJAkAR0DIcnhJBeTXE5ycpn9SXIqyVySx5McGFSb5M+TPJnkUpJ/SfLijVmSJGktBgZCkl3AFHA7sB+4s/8Fv7kD2APsBe4GHuhQ++6q2l9Ve4GHgHeufzmSpLXqcoRwELhUVfNV9SxwBphcMmcSmK6eC8COJGMr1VbV1/vqXwA8tc61SJLWYUeHOaPAfN/2AvCzHeaMDqpN8hfAm4D/pRcekqQh6RIIN01VvQN4R5K3Ae8FfmvpnCTHgeMAt95666b2p80zfvLcpj7elfuWHuRK6vKW0QIw1rc92sa6zOlSC/Ag8FPLPXhVna6qiaqaGBkZ6dCuJGktugTCeWBfktEkO4EjwMNL5swARwHaSeNrVTW/Um2S2/rq3wBcXNdKJEnrMvAto6p6Jsk9wCP0AmS6qmaTnGj7p4CzwKEkc8BV4NhKte2u39M+aroTuAK8eUNXJklalU7nEKpqht5RQP/YVN/tAu7tWtvG37iqTiVJN5XfVJYkAQaCJKkxECRJgIEgSWoMBEkSYCBIkhoDQZIEGAiSpMZAkCQBBoIkqTEQJEmAgSBJagwESRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSpMRAkSUDH31SWtpPxk+c27bGu3De5aY8lDeIRgiQJMBAkSY2BIEkCDARJUmMgSJIAA0GS1BgIkiSgYyAkOZzkYpLLSU4usz9JTiWZS/J4kgODapO8p41dTnIuye6NWZIkaS0GBkKSXcAUcDuwH7iz/wW/uQPYA+wF7gYe6FD7UeBlVfWjwEXgT9a9GknSmnX5pvJB4FJVzQMkOQNMAhf65kwC01VVwIUkO5KMAbfdqLaqPt5X/0ngTetejfQtZDO/MQ1+a1rd3jIaBeb7thfaWJc5XWoBjgP/vNyDJzmeZDbJ7OLiYod2JUlrMfSTykneAXwDmF5uf1WdrqqJqpoYGRnZ3OYkaRvpEggLwFjf9mgb6zJnxdokvwm8Hjja3m6SJA1Jl0A4D+xLMppkJ3AEeHjJnBngKEA7aXytnTe4YW2Sw8BbgddX1f9syGokSWs28KRyVT2T5B7gEXoBMl1Vs0lOtP1TwFngUJI54CpwbKXadtfvA3YBH0sC8OmqOrGhq5Mkddbp9xCqaobeUUD/2FTf7QLu7Vrbxl+yqk4lSTfV0E8qS5K2BgNBkgQYCJKkxt9UljSQvzO9PXiEIEkCDARJUmMgSJIAA0GS1BgIkiTAQJAkNQaCJAkwECRJjYEgSQIMBElS46UrJD1nbOYlNGD7XUbDIwRJEmAgSJIaA0GSBBgIkqTGQJAkAQaCJKkxECRJgIEgSWoMBEkSYCBIkhoDQZIEdAyEJIeTXExyOcnJZfYnyakkc0keT3JgUG2SX0lyKcm1JBMbsxxJ0loNDIQku4Ap4HZgP3Bn/wt+cwewB9gL3A080KH2Yqv7xPqXIUlary5HCAeBS1U1X1XPAmeApZcAnASmq+cCsCPJ2Eq1VXW5qj6/YSuRJK1Ll0AYBeb7thfaWJc5XWolSVvAlj+pnOR4ktkks4uLi8NuR5K+ZXUJhAVgrG97tI11mdOldkVVdbqqJqpqYmRkZDWlkqRV6PKLaeeBfUlGgS8DR4ATS+bMAHcBD7WTxteqaj7JYodaSXrO+Vb89baBgVBVzyS5B3iE3hHFdFXNJjnR9k8BZ4FDSeaAq8CxlWoBkrwR+BtgBDiX5Imq+vkNX6EkqZNOv6lcVTP0jgL6x6b6bhdwb9faNv4R4COraVaSdPNs+ZPKkqTNYSBIkgADQZLUGAiSJMBAkCQ1BoIkCTAQJEmNgSBJAgwESVJjIEiSAANBktQYCJIkwECQJDUGgiQJMBAkSY2BIEkCDARJUmMgSJIAA0GS1BgIkiTAQJAkNQaCJAkwECRJjYEgSQIMBElSYyBIkoCOgZDkcJKLSS4nObnM/iQ5lWQuyeNJDgyqTfK9ST6W5LNJHk3yPRuzJEnSWgwMhCS7gCngdmA/cGf/C35zB7AH2AvcDTzQofZdwMNV9TLg4bYtSRqSLkcIB4FLVTVfVc8CZ4DJJXMmgenquQDsSDI2oHYS+EC7Pb3MfUqSNlGXQBgF5vu2F9pYlzkr1Y5U1SJA+/eLurctSdpoO4bdwCBJjgPH2+Z/J/n8JrewG3h6tUV5903oZA1uUh/P6ecEbkovPifLW/Xz4nOyvHX2sqfLpC6BsACM9W2PtrHl5nx6yZydK9QuJhmpqsUkI8BXlnvwqjoNnO7Q502RZLaqJob1+FuRz8n1fE6W5/Nyva38nHR5y+g8sC/JaJKdwBF6J4H7zQBHAdpJ42tVNT+gdga4q92+a5n7lCRtooFHCFX1TJJ7gEfoBch0Vc0mOdH2TwFngUNJ5oCrwLGVattdvxM4k+S3gS8Dv7qxS5MkrUaqatg9bGlJjre3rdT4nFzP52R5Pi/X28rPiYEgSQK8dIUkqTEQbmDQ5Tq2oyRjST7Rnpd/T/LWYfe0VSR5frtsy0eH3ctWkOS7kzyU5Mkkn0vy08PuadiSvCvJfyT5fJKzSV4w7J6WMhCW0fFyHdvRs8DvVtU+4MeBNyd5xZB72ir+ALg87Ca2kL8F/qmq9gP7gEtD7meokrwEeBOwv6peCnwT+PXhdnU9A2F5XS7Xse1U1VNV9WS7/XXgSeAHh9vV8CUZpfffx98Nu5etIMkLgVdW1YMAVfWNqvrakNsatq/S+4Pq25PsAG4Bvjjclq5nICyvy+U6trUk48CrgE8Ot5Mt4a+BtwDXht3IFvFD9L54+lCSS0k+kOQ7h93UMFXVV4H76YXAfwFfq6pHh9vV9QwErVqS7wA+BPzhdv/LL8kvAl+pqseG3csW8jx6fyzcX1V76f11/KfDbWm4krwY+CPgNuAHgBckuWvlqs1nICyvy+U6tqX2jfOzwD9W1YeH3c8W8Grgl5JcAT4I/FyS6eG2NHTzwJeq6t/a9oeA7X6u6SeAf62qxfY29IeBnxlyT9cxEJbX5XId206SAH8PXK6qvxp2P1tBVb2tqkarahz4NeDjVbXl/vLbTO2yNU8neWkbei3wuSG2tBV8AfjJJLe0/49e28a2lC1/tdNhGHDJje3s1cBvAJ9N8kQbe3tVzQyxJ21NdwMPJvn/k6dHh9zPUFXV+SQfovdBjGvAE8D7h9vV9fymsiQJ8C0jSVJjIEiSAANBktQYCJIkwECQJDUGgiQJMBAkSY2BIEkC4P8Axqw6HEUrNl4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.bar(range(9),[fit_results[i][\"punzi_fom\"] for i in range(9)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAD0FJREFUeJzt3X+MXWldx/H3R1tJZAkROkik7Tb7DyEQxKWEbPxVIgm7LJgImHV3IZGsqUtWDRIJqAmQ6B8kyB9WTZpm+aFW1jUsAaPFAiFmwahkWsrSbRFjAmmNSwuYdc0GdnG//nHP4HU6c3/MnPtjnnm/ksnce85zzvn23NPPfe5zzzmTqkKS1JYfWHQBkqT+Ge6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBu1Z1Ib37dtXhw4dWtTmJWlHOnPmzDeramVcu4WF+6FDh1hdXV3U5iVpR0ry9UnaOSwjSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNWtgVqpLUrI/cNnr+HffPvAR77pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktSgseGe5ECSB5OcT/LVJO/YoE2SHEtyIckXk9w4m3IlSZOY5M/sPQn8WlU9lOQZwNkkp6vq3FCb1wHXAy8EfgL4EPDjvVcrSZrI2J57VT1SVQ91jx8DHgKet67ZrcDJGjgL7ElyoPdqJUkTmWrMPckh4GXA59fN2g9cGnp+uZu2fvmjSVaTrF69enW6SiVJE5s43JNcB3wUeGtVPbqVjVXViao6XFWHV1ZWtrIKSdIEJgr3JHuBB4D7qupjGzS5DAwPw+zvpkmSFmCSs2UCfAC4WFXv36TZKeDOrv2NwFNVdWmTtpKkGZvkbJmfBN4EfDnJ2hkyvwMcBKiq4wx69a9IcgF4AnjzDGqVJE1obLhX1eeBjGlTwD19FSVJ2h6vUJWkBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhq0Z9EFbMlHbhs9/47751OHJC0pe+6S1KCx4Z7kg0muJDm/yfwjSR5Ncq77eVf/ZUqSpjHJsMyHgT8G/mxEm89V1Wt6qUiStG1je+5V9SDw7TnUIknqSV9j7jclOZ/ks0le0tM6JUlb1MfZMmeAA1X1eJJXAR9PckNVPbW+YZKjwFGAgwcP9rBpSdJGtt1zr6rHqurx7vFp4AnguZu0PVFVh6vq8MrKynY3LUnaxLbDPcnK0OOXAtcBV7a7XknS1o0dlklyH3AE2JfkMvBuYC9AVR0Hbu+GW2DQa7+jqr43m3IlSZMYG+5VdfuY+ceAY71VJEnaNq9QlaQGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDVoz6ILkKQd6SO3LbqCkey5S1KDDHdJapDhLkkNMtwlqUFjwz3JB5NcSXJ+k/lJcizJhSRfTHJj/2VKkqYxSc/9w8DNI+a/DrgeeCFwF/Ch7ZclSdqOseFeVQ8C3x7R5FbgZA2cBfYkOdBXgZKk6fUx5r4fuDT0/HI3TZK0IHP9QjXJ0SSrSVavXr06z01L0q7SR7hfBoaHYfZ3065RVSeq6nBVHV5ZWelh05KkjfQR7qeAOwG6M2WeqqpLoxeRJM3S2HvLJLkPOALsS3IZeDewF6CqjgMPAK9IcgF4AnjzzKqVJE1kbLhX1e1j5hdwT28VSZK2zStUJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkH8gW9qCQ+/826mX+dp7b51BJZqZJf8D2OPYc5ekBhnuktQgw12SGuSYu6Yyaqx5XmPKy1CDtOwMd83FZoHcdxjPazvSsjPcpTnxE4fmqc1wH3cK0x33z6cOSVqQNsNdTdjKueTLvB1pnjxbRpIaZM9dG7I3u7v4fUB7DHctlG8i0mwY7tIu4pvp7mG4a1cw1HSNHX5jsHEMd/XGAJWWh+EuLQGvrFXfDHdpiW3lLBY/QQk8z12SmmTPfRezhye1a3eG+6hvyb3vjPT/eIHTzrQ7w32XsYfeJl9XjWK4S2pX4+eyj2K4S9oyT+FcXoa7pN45Tr94ngopSQ0y3CWpQYa7JDXIcJekBvmFqqS58gyb+TDcJe1cu/g89nEM9x3GqxIlTWKiMfckNyc5n+RiknduMP9IkkeTnOt+3tV/qZKkSY3tuSd5GnAc+GngEeAfk3yqqs6ua/q5qnrNDGqcr3Ef87yxmKQdYJKe+8uBh6vqUlU9CdwP+M2HJC2xScbc9wOXhp5fBo5s0O6mJOeBK8Dbqurc+gZJjgJHAQ4ePDh1sUvB2wVL2gH6+kL1DHCgqh5P8irg40luqKqnhhtV1QngBMDhw4erp21LasBmJwvcu/d9vPIFPzrnana+SYZlLgMHhp7v76Z9X1U9VlWPd49PA08Az+2rSEnSdCYJ9y8AL0qyP8le4Dbgk8MNkqwMPX4pcB2D4RlJ0gKMHZapqu8keQtwmsGbwcmqWk1ydzf/OHB7N54Og177HVX1vVkVLUkabaIx96o6BZxaN+340ONjwLF+S5MkbZVXqEpaqHv3vm/RJTTJcO+TF0BJWhKG+xLy/jGStsv7uUtSg+y5z9O6YZvPXPzGhs3u3XvttF958u2zqEhSowz3HWLcl06Gv6RhDstIUoPsuUuaKU91XAzDvREO22iWDOidx3DfJUb95xwX/NtZVtJiGO6yVyaPgQYZ7tqWRQ4H+YlC2pzhLvVsWd907J3vLoa7Zmo7PfvthNEsP1Essq7dGtCbXfDnX2janOE+Y5sdlBpY1rBa1rqkSRnu2pUMb7XOcO+JPXT1wTed6Yz6f7fbh2y8/YAkNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQV7EJKlJu/1+NIb7FLwKVdJO4bCMJDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yCtU1/EqVEktmCjck9wM/AHwg8CfVtV7180P8IfAK4HvAndV1dmea5Wkbdstf1R77LBMkqcBx4FbgBcDb0hy47pmrwOuB14I3AV8qOc6JUlTmKTn/nLg4aq6BJDkfuBWYLhnfitwsqoKOJtkT5IDa8tI0k7Q0p0kJ/lCdT8wHNKXu2nTtpEkzclcv1BNchQ42j397yT/ssVV7QO+2U9VvVrWumB5a7Ou6VjXdJazrjv/ajt1XT9Jo0nC/TJwYOj5/m7aRm3+aUQbquoEcGKSwkZJslpVh7e7nr4ta12wvLVZ13Ssazq7ua5JhmW+ALwoyf4ke4HbgE+ua3MKuBOg+7L1KcfbJWlxxvbcq+o7Sd4CnGbwZnCyqlaT3N3NPw48ALwiyQXgCeDNM6xZkjTGRGPuVXWKQe98eNrxoccF3NNvaSNte2hnRpa1Llje2qxrOtY1nV1bVwa5LElqifeWkaQGLW24J/nFJA8neSrJpt8qJ7k5yfkkF5O8c2j6s5J8OsmXk3wqyY/0VNfY9SZ5fpJzQz//leSt3bz3JPn3oXmvnlddXbuvdW3OJVmddvlZ1JXkQJIHu9fxq0neMTSv1/212fEyND9JjiW5kOSLw1djj1t2xnW9qduH55OcGf4/sdlrOqe6jiR5dOj1edeky864rrcP1XQ+yf8keVY3b5b764NJriQ5v8n8+R1fVbWUP8ALgOcDfw8c3qTN04CvMTgNcy+wCtzYzfsj4G3d498EjvVU11TrZXA/nkeA67vn7wF+awb7a6K6uv21b7v/rj7rAp4LvLh7/AzgX4GX9L2/Rh0vQ21eD3wCCHAj8KVJl51xXS8Hntk9vgU4N+41nVNdR4C/2cqys6xrXfvXAp+d9f7q1v0z3XFzfpP5czu+lrbnXlUXq2rcRU7fvzVCVT0JrN0age73n3ePTw5N365p1/tzwL9V1dd72v5mtvvvXdj+qqpHquqh7vFjwEPA83ra/rBRx8twvSdr4CywJ8mBCZedWV1V9c9V9Wj39PPMZv9MXdeMlu173bcD9/W07ZGq6kHg2yOazO34Wtpwn9Co2x6sVNVVgO73c3ra5rTr/SWuPbDuSfKVJH+R5NlzrquAtWGSX9/C8rOqC4Akh4CXMQiwNX3tr+3cSmOWt9iYdt2/Cvz10PPNXtN51XVTN5zw2SQvmXLZWdZFkh8GbmZwuvaaWe2vSczt+Fro/dyTfIbBR/L1freqPjHvetaMqmvK9fwQ8PPAbw9N/hPg9xgcYO8BjtFdADanum6qqkeSPAf4uyRfqapPT7H8rOoiyXXAR4G3DvVSt7y/WpTkCIM7r/7U0OTeX9MpnAEOVNXjSV4FfDzJDXPa9iReC/xDVQ33phe5v+ZmoeFeVa/c5ipG3RrhapKVqrqaZAW40kddSaZZ7y3A2ar6/q3m1nqx3bqOM/hOYW51VdUj3e8rST7KoJf8aRa8vzK4+vkB4L6q+tjQure8vzawnVtp7J1g2VnWRZIXAx8Abqmqb61NH/Gazryubhht7fHpJE8weKOf6N80q7qGXPPJeYb7axJzO752+rDMqFsjnALe2D1+I9feMmGrplnvNWN9XW9hzeuBC/OqK8nTu4+pJHk6g4+rFyZdfoZ1hUFoXayq96+b1+f+2s6tNCZZdmZ1JTkIfAx4U1V9dWj6qNd0HnWtDD1+KXAdgzfwhe6vrp5nAj/L4AvMtWmz3F+TmN/xNYtvjPv4AX6BwTvXd4FvAKe76T8GnBpq92rgYeAig+GctenPBj4DfLn7/aye6tpwvRvU9XTgW3RnOAxNP8ngC8OvAJ9i8JF2LnUBN3Tb/hKDM1J+n/+7kG1h+4vBEEN1tZ3rfl49i/210fEC3A3c3T0Og6GgC10dh0ct2+PxPq6ue4H/HNo/q+Ne0znV9RvA+e7nLHBkGfZX9/yXgb9ct9ys99d9wH8ATzLIr7sWdXx5haokNWinD8tIkjZguEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1KD/BQAh8Rrk7RZYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(data[0][\"cos_thetal\"],density=True,bins=40);\n",
    "plt.hist(data_dict[\"cos_thetal\"],alpha=0.7,density=True,bins=40);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "iteration=0\n",
    "variable=\"Ds_ConsD_M\"\n",
    "\n",
    "variable_MC_Ds=variable\n",
    "variable_MC_Dplus=variable_MC_Ds.replace(\"Ds\",\"Dplus\")\n",
    "\n",
    "inf=data[iteration][variable].min()\n",
    "sup=data[iteration][variable].max()\n",
    "\n",
    "data_entries=data[iteration][variable].shape[0]\n",
    "\n",
    "mc_Ds_entries=MC_Ds[iteration][variable_MC_Ds].shape[0]\n",
    "mc_Dplus_entries=MC_Dplus[iteration][variable_MC_Dplus].shape[0]\n",
    "\n",
    "f=r.TFile(\"test.root\",\"RECREATE\")\n",
    "t=r.TTree(\"histo\",\"histo\")\n",
    "\n",
    "data_hist = r.TH1F(\"sWeighted data \"+ variable, \n",
    "                   \"sWeighted data \"+ variable, \n",
    "                   70,inf,sup)\n",
    "data_hist.Sumw2()\n",
    "\n",
    "mc_hist = r.TH1F(\"MC \"+ variable, \n",
    "                 \"MC \"+ variable, \n",
    "                   70,inf,sup)\n",
    "mc_hist.Sumw2()\n",
    "\n",
    "for i in range(data_entries):\n",
    "    data_hist.Fill(data[iteration][variable][i],data[iteration][\"nDs_sw\"][i])\n",
    "        \n",
    "for i in range(mc_Ds_entries):    \n",
    "    mc_hist.Fill(MC_Ds[iteration][variable_MC_Ds][i])\n",
    "    \n",
    "for i in range(mc_Dplus_entries):    \n",
    "    mc_hist.Fill(MC_Dplus[iteration][variable_MC_Dplus][i])\n",
    "\n",
    "\n",
    "n1 = data_hist.Integral(\"width\") \n",
    "data_hist.Scale(1/n1)\n",
    "n2 = mc_hist.Integral(\"width\")\n",
    "mc_hist.Scale(1/n2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_hist.Draw(\"E\")\n",
    "#mc_hist.Draw(\"E Same\")\n",
    "data_hist.GetXaxis().SetTitle(variable)\n",
    "data_hist.GetYaxis().SetTitleOffset(1.5)\n",
    "data_hist.GetYaxis().SetTitle(\"distr dN/d\"+variable)\n",
    "data_hist.SetLineColor(r.kRed)\n",
    "r.gStyle.SetOptStat(\"11\")\n",
    "legend = r.TLegend(0.89, 0.89, 0.75, 0.8)\n",
    "legend.AddEntry(data_hist,\"data\",\"f\")\n",
    "legend.AddEntry(mc_hist,\"mc\",\"f\")\n",
    "legend.Draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#for variable in branches_needed:\n",
    "#\n",
    "#    plot_MC_vs_data(data_cut, data_weights, MC_Ds_cut, \n",
    "#                variable, \n",
    "#                sw_idx=sw_idx, l_index=l_index, mother_index_fit=mother_index_fit)\n",
    "#plt.hist(data_cut[\"cos_thetal\"], weights=data_weights_Ds, density=True,\n",
    "#         histtype='step', hatch='/', fill='True', linewidth='2', \n",
    "#         alpha=0.65,  bins=70, label='Weighted data Cos(theta_l)');\n",
    "#\n",
    "#plt.hist(MC_Ds_cut[\"cos_thetal\"],density=True, \n",
    "#         histtype='step', hatch='\\\\', fill='True', linewidth='2', \n",
    "#         alpha=0.45, bins=70, label='MC Ds Cos(theta_l)');\n",
    "#plt.legend(fontsize=12)\n",
    "#fig=plt.gcf()\n",
    "#fig.set_size_inches(10,6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#sw_data_array= np.array( \n",
    "#                           [\n",
    "#                           tuple( data[label][k] for label in branches+[\"nDplus_sw\",\"nDs_sw\"] )\n",
    "#                               \n",
    "#                           for k in range(s_weights[0][\"nDplus_sw\"].shape[0])\n",
    "#                           ],\n",
    "#    \n",
    "#                           dtype=[(label, np.float32) for label in branches+[\"nDplus_sw\",\"nDs_sw\"]]\n",
    "#    \n",
    "#                         )\n",
    "# rn.array2root(sw_data_array,\n",
    "#                      filename='/disk/lhcb_data/davide/Rphipi/BDT_selected_data/'+l_flv[l_index]+l_flv[l_index]+'/sw_data_'+l_flv[l_index]+l_flv[l_index]+'.root',\n",
    "#                      treename='decay_tree',\n",
    "#                      mode='recreate',\n",
    "#                     )"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}