diff --git a/BDT_select_and_fit.ipynb b/BDT_select_and_fit.ipynb index 7098684..c31424a 100644 --- a/BDT_select_and_fit.ipynb +++ b/BDT_select_and_fit.ipynb @@ -28,6 +28,15 @@ ] }, { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "#b = r.TBrowser()" + ] + }, + { "cell_type": "markdown", "metadata": {}, "source": [ @@ -36,15 +45,19 @@ }, { "cell_type": "code", - "execution_count": 2, + "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", - "test=4\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", @@ -52,7 +65,12 @@ "with open(BDT_PATH+'/variables_used.pickle', 'rb') as f: \n", " branches_needed_0 = pickle.load(f)\n", " \n", - "branches_needed=['cos_thetal']+branches_needed_0" + "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)" ] }, { @@ -64,21 +82,22 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 4, "metadata": { - "scrolled": true + "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "(34845, 2035, 1295)\n" + "14942\n", + "(13679, 1155, 1280)\n" ] }, { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEcdJREFUeJzt3W2MpWV9x/HvTxd5o21Ehpg4C9vExFgQFZZSQ0VpaQps06aIEsO2LqwlS2gt9o0PTUTSN74QjViTDcnyELYhm4pNU7ME8EVDSDAwLAj7YLVNkF0jMGJjaFIe4v774lwrZ8czM/eZp3POzPeTTPZ+Pv9z7Zn7d677OueeVBWSJL1p1AVIksaDgSBJAgwESVJjIEiSAANBktQYCJIkwECQJDUGgiQJMBAkSc2mURcwjNNPP722bNky6jIkaaI88cQTP6+qqcW2m6hA2LJlCzMzM6MuQ5ImSpKfdNnOS0aSJMBAkCQ1BoIkCTAQJEmNgSBJAgwESVJjIEiSAANBktQYCJIkYMK+qaz1Yeddj580v2fHBSOqRFI/ewiSJMBAkCQ1BoIkCTAQJEmNg8paE3MHkiWNH3sIkiTAQJAkNQaCJAkwECRJjYEgSQIMBElS0ykQklyW5GCSI0k+P2B9ktyW5HCSJ5Oc15ZvTvJw2/dHST7Xt89pSR5K8kySB5O8feWeliRpWIsGQpJTgd3A5cC5wFUnTvh9rgTOAs4GdgJ3tuWvA39TVecA5wOfTvKBtu4W4P6qeh9wf5uXJI1Ilx7ChcChqjpaVa8D+4Btc7bZBuytngPApiSbq+r5qnoaoKpeBp4G3tW3zz1teu+AY0qS1lCXQJgGjvbNH2vLhtomyRbgAuCRtmiqqmYB2r9nDHrwJNcnmUkyMzs726FcSdJSrMmgcpK3At8GbqqqXw6zb1XdXlVbq2rr1NTU6hQoSeoUCMeAzX3z021Zp22SnALcB9xbVd/p22Y2yVTbZgp4cbjSJUkrqUsgPAack2S6ndyvpjcI3G8/cA1AG3A+XlVHkwTYAxypqlsH7LO9TW8fcExJ0hpa9G6nVfVKkhuAB+gFyN6qmkmyq63fTa8HcEmSw8BrwLVt94uAvwSeSfJUW/bFqtoP3AzsS3Id8ALwiRV8XpKkIXW6/XU7ge+fs2x333QBNw7Y7xEg8xzzJeDSYYqVJK0ev6ksSQIMBElSYyBIkgADQZLUGAiSJMBAkCQ1BoIkCTAQJEmNgSBJAgwESVJjIEiSAANBktQYCJIkwECQJDUGgiQJMBAkSY2BIEkCDARJUmMgSJIAA0GS1BgIkiTAQJAkNQaCJAkwECRJjYEgSQIMBElSYyBIkgADQZLUGAiSJMBAkCQ1BoIkCTAQJEmNgSBJAgwESVJjIEiSAANBktQYCJIkwECQJDUGgiQJ6BgISS5LcjDJkSSfH7A+SW5LcjjJk0nO61t3R5IXkxycs8+Xk/w0yVPt54rlPx1J0lItGghJTgV2A5cD5wJX9Z/wmyuBs4CzgZ3AnX3r7gIum+fwX6+qD7Sf/UPWLklaQZs6bHMhcKiqjgIk2QdsAw70bbMN2FtVBRxIsinJ5qo6WlUPJ9mywnVrHdl51+Mnze/ZccGIKpE2ti6XjKaBo33zx9qyYbcZ5MYkP0zyz0neMWiDJNcnmUkyMzs72+GQkqSlGOWg8reAdwO/C/w3cNugjarq9qraWlVbp6am1rI+SdpQugTCMWBz3/x0WzbsNiepqtmq+lVVHac3RuF1AkkaoS6B8BhwTpLpJKcAVwP3z9lmP3ANQBtwPn5izGE+Sc7om/0YcLhz1ZKkFbfooHJVvZLkBuABegGyt6pmkuxq63cD9wGXJDkMvAZce2L/JPcCHwVOT3IMuLmq9gBfS3Iu8BbgOXqfTpIkjUiXTxnRPhK6f86y3X3TBdw4z76fnGf59u5lSpJWm99UliQBBoIkqTEQJElAxzEEaVhzv30safzZQ5AkAQaCJKkxECRJgIEgSWoMBEkSYCBIkhoDQZIEGAiSpMZAkCQBBoIkqTEQJEmAgSBJagwESRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSpMRAkSYB/U1maOIv9veo9Oy5Yo0q03thDkCQB9hC0QhZ71ypp/NlDkCQBBoIkqTEQJEmAgSBJagwESRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSpMRAkSYCBIElqOgVCksuSHExyJMnnB6xPktuSHE7yZJLz+tbdkeTFJAfn7HNakoeSPJPkwSRvX/7TkSQt1aKBkORUYDdwOXAucFX/Cb+5EjgLOBvYCdzZt+4u4LIBh74FuL+q3gfc3+YlSSPSpYdwIXCoqo5W1evAPmDbnG22AXur5wCwKclmgKp6GPjFgONuA+5p03sHHFOStIa6BMI0cLRv/lhbNuw2c01V1SxA+/eMDrVIklbJ2A8qJ7k+yUySmdnZ2VGXI0nrVpdAOAZs7pufbsuG3Wau2SRTAO3fFwdtVFW3V9XWqto6NTXVoVxJ0lJ0CYTHgHOSTCc5Bbia3iBwv/3ANQBtwPl4VR1lYfuB7W16+4BjSpLW0KKBUFWvADcADwBPA/9aVTNJdiXZ1Ta7D/hpksPAHcC1J/ZPci/wKPCeJMeS7Gyrbga2JXmG3oDyl1bqSUmShrepy0ZVtZ/eO/r+Zbv7pgu4cZ59PznP8peASztXKklaVWM/qCxJWhsGgiQJMBAkSY2BIEkCDARJUtPpU0aSRmfnXY+PugRtEPYQJEmAPQSNobnviPfsuGBElUgbiz0ESRJgIEiSGgNBkgQYCJKkxkFlaZ1xUF5LZQ9BkgQYCJKkxkCQJAEGgiSpMRAkSYCBIElqDARJEmAgSJIaA0GSBBgIkqTGW1dIY8i/kqZRsIcgSQIMBElS4yUjaUJd/PJ3f2PZw2/70xFUovXCHoIkCbCHIE2MQT0CaSXZQ5AkAQaCJKnxkpG0jsy9rOQgs4ZhIEjrXP+X3Pz7ylqIgaBO/MPt0vrnGIIkCTAQJEmNgSBJAhxDkMaWX0TTWrOHIEkCOvYQklwGfBV4M3B3VX1lzvoA3wAuBV4FdlbVgYX2TfJl4K+B2XaYL1bV/uU+Ia0N79cvrT+LBkKSU4HdwIeB54FHkzx44oTfXAmcBZwNfBC4E3h/h32/XlVfXbFnI0lasi6XjC4EDlXV0ap6HdgHbJuzzTZgb/UcADYl2dxxX0nSGOhyyWgaONo3fwz4aIdtpjvse2OSTwNPAJ+pqpc6Va0NxS/FSWtjlJ8y+hbwj0ABXwZuA66Zu1GS64HrAc4888w1LE+afL/5SSXDVPPrcsnoGLC5b366Leuyzbz7VtVsVf2qqo7TG2cY+EqtqturamtVbZ2amupQriRpKbr0EB4DzkkyDbwAXA3smrPNfmA78C9JzgOOV9XRJLPz7ZvkjKp6se3/MeDwsp+NNMlm7vz15MUvPzu6OrRhLRoIVfVKkhuAB+j1KPZW1UySXW39buA+4JIkh4HXgGsX2rcd+mtJzgXeAjwH7FzZpyZJGkanMYT2/YD9c5bt7psu4Mau+7bl24eqVBPPe/XPb+ddj9sr0Mh56wqtGm+9sAgvEWnMGAgamUGBYa9BGh0DQSvGHsEE6OuVALD12tHUobFkIEhrYe6JWBpD3u1UkgTYQ9CYWTefRLJHoAlkD0GSBNhDkJZmvfQAHGRWHwNB8xrHP4Jz8cvf5e5vvnFZ6VN/e8sIq5HWFy8ZSZIAewhaorX6zsHQjzPoUs5yL4Osl8tD0iIMBG08XjeXBjIQNNHu/ubNJ81/6kNbln9QewTaoAwE/do4DiKvCE/wUicGgmRgSICBoI68cd3y3P3osyfNr8ilrSVYtI4u4eiYy7plIGggA0BLthqf9NKaMBAkDcdLbOuWgaCemTv9q13SBmcgbFTr9F3eKK/V9z/2qMYIhjUuYxsaDwbCRrBOT/6aYH45cCwZCJJW32JvSgyIsWAgrAf2ALTeGBAjYSBMAk/4E2futXlpEhgI48aT/9hayQFYA2NIS/nCnL2MoRkIG5gnJUn9DIS15ruWNeXHKjcQB66XzUBYjpW474uXiKTR6HKLjQ0WIgbCavOErwmy4XtUG7yXYSAMw5P7urbhT4Ya3jq7O+zGDQTvyChpFMa4l7FxA2GQdd4D8FNFK9sGtqc6maDzysYJhAn6T9F48ISvjeZNoy5AkjQeDARJErCRLhmtUwt9MsZLHlppk/g3H9SdgSBpXgu9qVjsDYeBMXkMBEkj53dAxoNjCJIkoGMPIcllwFeBNwN3V9VX5qwP8A3gUuBVYGdVHVho3ySnAfuAdwI/A66uqv9ZiSe1kTluIGmpFg2EJKcCu4EPA88DjyZ58MQJv7kSOAs4G/ggcCfw/kX2vQW4v6q+luSzbf4zK/fU1idP+JoUy7kM5CWk0ejSQ7gQOFRVRwGS7AO2Af2BsA3YW1UFHEiyKclm4HcW2HdbOzbAXuD7GAjShjEub24Mnzd0CYRp4Gjf/DHgox22mV5k36mqmgWoqtkkZ3Suep1Z6KN84/JLIy3Xcl7Lw+y72O/QYif8lfx9XE647Lzr8ZPm9+y4YMnH6mrsP2WU5Hrg+jb7v0n+c4mHOh34+cpUtXp2TEidTE6dMDm1WucK2PHG5MA6d8xd0O1Yy61lAdd1as87lncPvLO6bNQlEI4Bm/vmp9uyQdt8f842pyyw72ySqdY7mAJeHPTgVXU7cHuHOheUZKaqti73OKvNOlfepNRqnSvLOofX5WOnjwHnJJlOcgpwNXD/nG32A9cAJDkPON7GDRbadz+wvU1vH3BMSdIaWrSHUFWvJLkBeIBegOytqpkku9r63cB9wCVJDgOvAdcutG879M3AviTXAS8An1jZpyZJGkanMYSq2k/vHX3/st190wXc2HXftvwlet9bWCvLvuy0Rqxz5U1Krda5sqxzSOmdyyVJG523rpAkARMeCEnuSPJikoN9y/Ylear9PJvkqbZ8S5L/61u3u2+f85M8meRwktvarThWu86LkvwgyaEkTye5qG/dF5IcSXIwyZ+MY51j2J5bkxxodf57kt/qWzdO7TmwzhG35+YkD7f2+VGSz7XlpyV5KMkzSR5M8va+fda8TYetc1RtukCdH2//78eTbJ2zz0heo7+hqib2B7gYOA84OM/6W4EvtektC2z3NHB+m/434MrVrhN4BLi8TV8BPNKmzwdm6H1kdxp4Fjh1DOsct/Z8BvhIm74OuHVM23O+OkfZnu8Ezm3TbwN+DHwA+Cbw9235Z4HbRtmmS6hzJG26QJ3vBd4D/AewtW/7kb1G5/5MdA+hqh4GfjFoXUvSTwD3LnSMJGcCb66qJ9qivfRuq7HadR4DTryL/W3guTa9DdhXVa9X1THgEPB7Y1jnQCOs893Aw236IeDP2vS4ted8dQ60RnU+X1VPt+mX6Z2E3tUe554BjzuSNl1CnQONqs6qOlJVg75YO7LX6FwTHQiL+DDwQlX9uG/Zlnb549Ekf9SWzXfbjdX2OeDWJEfp3Q32C4vUM251wni15xHgz9v0x4EzF6ln3OqEMWjPJFuAC+j1DE+6vQxw4vYyI2/TjnXCiNt0Tp3zGXl7nrCeA+GTnNw7+BkwXVXvp/cR2Xv6r4mOwB7g76pqM71u7p4R1rKQ+eoct/b8K+Cmdr3+dHq3YR9H89U58vZM8lbg28BNVfXLtXzsYQxR50jbdFLas9/Y38toKZJsondL7vNPLKuqV2m/fFV1oP1Cvpdut+ZYDR8C/rhN/wu9W4azQD1jVee4tWdVHaTdOLG9K7uirRqr9pyvzlG3Z3p3ErgPuLeqvtMWz3d7mZG16TB1jrJN56lzPmPzGl2vPYRLgR+263EAJHlHkje16S3AOcB/VdVzwPH0brkBvVtwrMVtNH4CfKRN/yG9gSTofYnv6iSnJJludT42bnWOW3smOb39G+CLvNGTGav2nK/OUbZnq2UPcKSqbu1bNd/tZUbSpsPWOao2XaDO+YzPa3Q1R6xX+4feJaGfAa/TS86dbfldwK45215Fb7DmGeAg8PG+dVuBp4DDwD/RvrC3mnUCFwE/aI/5FPD7fdv/A71rzYdon/AZtzrHsD1vAn7YavlK/2OOWXsOrHPE7fkHQNEb/Hyq/VwBvAP4Xqvpe8Bpo2zTYescVZsuUOdftNfBq/Ru1/PAqF+jc3/8prIkCVi/l4wkSUMyECRJgIEgSWoMBEkSYCBIkhoDQZIEGAiSpMZAkCQB8P9XpuL7eLeHHAAAAABJRU5ErkJggg==\n", + "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": [ "
" ] @@ -91,15 +110,68 @@ ], "source": [ "#MC_Dplus_sig_dict, MC_Ds_sig_dict, _ = load_datasets(l_index)\n", + "\"\"\"\n", + "Data\n", "\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]=='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", + "\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", @@ -121,10 +193,10 @@ " upper_cut=2010\n", " if mother_index_fit==1:\n", " #fitting only the Dplus\n", - " lower_cut=1810\n", - " upper_cut=1930\n", + " lower_cut=1860\n", + " upper_cut=1880\n", " if mother_index_fit==2:\n", - " lower_cut=1810\n", + " lower_cut=1820\n", " upper_cut=2040\n", "\n", "if l_index==0:\n", @@ -134,15 +206,15 @@ " upper_cut=2040\n", " if mother_index_fit==1:\n", " #fitting only the Dplus\n", - " lower_cut=1780\n", - " upper_cut=1930\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, \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", @@ -153,11 +225,13 @@ "\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, mother_index_fit=mother_index_fit)\n", - " \n", - " \n", - "plt.hist(np.concatenate((mc_Dplus_mass, mc_Ds_mass)),bins=70,density=True,alpha=0.7);\n", - "plt.hist(data_mass,bins=70,density=True,alpha=0.4);\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" ] }, { @@ -169,18 +243,20 @@ }, { "cell_type": "code", - "execution_count": 4, - "metadata": {}, + "execution_count": 5, + "metadata": { + "scrolled": true + }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "BDT cut at 0.1 fit\n", + "BDT cut at 0.0615789473684 fit\n", "BDT selection..\n", - "(340, 476)\n", - "BDT Signal selection efficiency: 0.245045045045\n", - "chi2 1.10479563477\n", + "(499, 474)\n", + "BDT Signal selection efficiency: 0.399589322382\n", + "chi2 0.747467122077\n", "Saving fit results...\n" ] } @@ -194,91 +270,161 @@ " # 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.10,0.90,num=10)\n", + "x_cut_values = np.linspace(0.01,0.99,num=20)\n", + "data={}\n", + "MC_Dplus={}\n", + "MC_Ds={}\n", "fit_results={}\n", - "splots={}\n", "#Just fit it!\n", - "for i, x_cut in enumerate(x_cut_values[0:1]):\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], splots[i] = select_and_fit(BDT_PATH, mother_index_fit, l_index, lower_cut, upper_cut,\n", - " data_mass, mc_Dplus_mass, mc_Ds_mass, data_to_select, mc_Dplus_to_select, mc_Ds_to_select,\n", - " data_dict, MC_Dplus_tuple_dict, MC_Ds_tuple_dict, branches_needed,\n", - " x_cut, i, test)\n", - "\n", - " with open(FIT_PATH+'/'+mother_ID[mother_index_fit]+'/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))\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": 5, - "metadata": {}, - "outputs": [], - "source": [ - "sw_idx=0\n", - "with open('/disk/lhcb_data/davide/Rphipi/BDT_selected_data/'+l_flv[l_index]+l_flv[l_index]+'/BDT_sel_data_'+l_flv[l_index]+l_flv[l_index]+'_geom_var.pickle', 'rb') as f:\n", - " data_cut=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[1]+'_'+l_flv[l_index]+l_flv[l_index]+'_geom_var.pickle', 'rb') as f:\n", - " MC_Dplus_cut=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]+'_geom_var.pickle', 'rb') as f:\n", - " MC_Ds_cut=pickle.load(f)\n", - "\n", - "key=data_cut.keys()[0]\n", - "nentries=data_cut[key].shape[0]" - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "metadata": {}, - "outputs": [], - "source": [ - "data_weights={}\n", - "if mother_ID[mother_index_fit]=='Dplus':\n", - " data_weights[0]=np.array([splots[sw_idx].GetSWeight(i, \"nDplus_sw\") for i in range(nentries)])\n", - "\n", - "if mother_ID[mother_index_fit]=='Ds':\n", - " data_weights[1]=np.array([splots[sw_idx].GetSWeight(i, \"nDs_sw\") for i in range(nentries)])\n", - "\n", - "if mother_ID[mother_index_fit]=='both':\n", - " data_weights[0]=np.array([splots[sw_idx].GetSWeight(i, \"nDs_sw\") for i in range(nentries)])\n", - " data_weights[1]=np.array([splots[sw_idx].GetSWeight(i, \"nDplus_sw\") for i in range(nentries)])\n" - ] - }, - { - "cell_type": "code", - "execution_count": 7, + "execution_count": 23, "metadata": {}, "outputs": [ { - "name": "stdout", - "output_type": "stream", - "text": [ - "Saved: e/fits/both/0/Ds_sweighted_cos_thetal_0.png\n", - "Saved: e/fits/both/0/Ds_sweighted_Ds_ENDVERTEX_CHI2_0.png\n", - "Saved: e/fits/both/0/Ds_sweighted_Ds_IPCHI2_OWNPV_0.png\n", - "Saved: e/fits/both/0/Ds_sweighted_Ds_OWNPV_CHI2_0.png\n", - "Saved: e/fits/both/0/Ds_sweighted_Ds_DIRA_OWNPV_0.png\n", - "Saved: e/fits/both/0/Ds_sweighted_Ds_PT_0.png\n", - "Saved: e/fits/both/0/Ds_sweighted_Ds_FD_OWNPV_0.png\n", - "Saved: e/fits/both/0/Ds_sweighted_Ds_FDCHI2_OWNPV_0.png\n", - "Saved: e/fits/both/0/Ds_sweighted_pi_PT_0.png\n", - "Saved: e/fits/both/0/Ds_sweighted_e_plus_PT_0.png\n", - "Saved: e/fits/both/0/Ds_sweighted_e_minus_PT_0.png\n" - ] + "data": { + "text/plain": [ + "" + ] + }, + "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": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" } ], "source": [ - "for variable in branches_needed:\n", + "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": [ + "
" + ] + }, + "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", - " 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)" + "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()" ] }, { @@ -287,6 +433,11 @@ "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", @@ -298,6 +449,29 @@ "#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": {