diff --git a/.ipynb_checkpoints/bkg_reduction-checkpoint.ipynb b/.ipynb_checkpoints/bkg_reduction-checkpoint.ipynb new file mode 100644 index 0000000..e908b54 --- /dev/null +++ b/.ipynb_checkpoints/bkg_reduction-checkpoint.ipynb @@ -0,0 +1,264 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "import ROOT as r\n", + "import ctypes\n", + "import numpy as np\n", + "from array import array" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "r.TBrowser()" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "data_mumu = r.TFile(\"/disk/lhcb_data/davide/Rphipi/small/data/Ds_phipi_mumu/Ds_phipi_mumu_MagDown.root\")\n", + "MC_mumu = r.TFile(\"/disk/lhcb_data/davide/Rphipi/MC/Ds_phipi_mumu/Ds_phipi_mumu_MagDown.root\")" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "t_data_mumu = data_mumu.Get(\"Ds_OfflineTree/DecayTree\")\n", + "t_data_mumu" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "t_MC_mumu = MC_mumu.Get(\"Ds_OfflineTree/DecayTree\")\n", + "t_MC_mumu" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [], + "source": [ + "t_data_mumu.SetBranchStatus(\"*\",0)\n", + "t_data_mumu.SetBranchStatus(\"Ds_ENDVERTEX_CHI2\",1)\n", + "t_data_mumu.SetBranchStatus(\"Ds_ENDVERTEX_NDOF\",1)\n", + "t_data_mumu.SetBranchStatus(\"Ds_OWNPV_CHI2\",1)\n", + "t_data_mumu.SetBranchStatus(\"Ds_OWNPV_NDOF\",1)\n", + "t_data_mumu.SetBranchStatus(\"Ds_IP_OWNPV\",1)\n", + "t_data_mumu.SetBranchStatus(\"Ds_IPCHI2_OWNPV\",1)\n", + "t_data_mumu.SetBranchStatus(\"Ds_DIRA_OWNPV\",1)\n", + "t_data_mumu.SetBranchStatus(\"Ds_ConsD_M\",1)\n", + "\n", + "t_data_mumu.SetBranchStatus(\"Ds_Hlt1Phys_TOS\",1)\n", + "t_data_mumu.SetBranchStatus(\"Ds_Hlt1TrackMVADecision_TOS\",1)\n", + "#t_data_mumu.SetBranchStatus(\"Ds_Hlt2RareCharmD2PiMuMuOSDecision_TOS\",1)\n", + "\n", + "t_MC_mumu.SetBranchStatus(\"*\",0)\n", + "t_MC_mumu.SetBranchStatus(\"Ds_ENDVERTEX_CHI2\",1)\n", + "t_MC_mumu.SetBranchStatus(\"Ds_ENDVERTEX_NDOF\",1)\n", + "t_MC_mumu.SetBranchStatus(\"Ds_OWNPV_CHI2\",1)\n", + "t_MC_mumu.SetBranchStatus(\"Ds_OWNPV_NDOF\",1)\n", + "t_MC_mumu.SetBranchStatus(\"Ds_IP_OWNPV\",1)\n", + "t_MC_mumu.SetBranchStatus(\"Ds_IPCHI2_OWNPV\",1)\n", + "t_MC_mumu.SetBranchStatus(\"Ds_DIRA_OWNPV\",1)\n", + "t_MC_mumu.SetBranchStatus(\"Ds_ConsD_M\",1)\n", + "\n", + "t_MC_mumu.SetBranchStatus(\"Ds_Hlt1Phys_TOS\",1)\n", + "t_MC_mumu.SetBranchStatus(\"Ds_Hlt1TrackMVADecision_TOS\",1)\n", + "#t_MC_mumu.SetBranchStatus(\"Ds_Hlt2RareCharmD2PiMuMuOSDecision_TOS\",1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#end, start = fit region" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "def fcn( npar, gin, f, par, iflag):\n", + " \n", + " global ncount\n", + " global chisq\n", + " global pulls\n", + " global fitfunavg\n", + " global delta\n", + " \n", + " pull = np.array([0. for i in range(nbins)])\n", + " residuals = np.array([0. for i in range(nbins)])\n", + " \n", + " fitfunavg=np.array([func((xup[i]+xdown[i])/2) for i in range(nbins)])\n", + " \n", + " chisq, delta = 0., 0.\n", + " \n", + " delta = np.array([data[i]-fitfun[i] for i in range(0,nbins)])\n", + " \n", + " chisq = np.dot(np.transpose(delta),np.dot(Vxinv,delta))\n", + " \n", + " for i in range(nbins):\n", + " \n", + " pulls[i] = (y[i]-fitfun[i])/(np.sqrt(Vx[i,i]))\n", + " \n", + " print('X^2', chisq)\n", + " f[0] = chisq\n", + " ncount+=1" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "###DEFINE FIT FUNCTIONS HERE\n", + "###AS FUNC" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "##___MINUIT______\n", + "\n", + "def fit():\n", + " \n", + " gMinuit = r.TMinuit(5)\n", + " gMinuit.SetFCN(fcn)\n", + " \n", + " arglist = np.array([0. for i in range(10)], dtype='float32')\n", + " ierflg = ctypes.c_int(1982)\n", + " \n", + " arglist[0]=1\n", + " gMinuit.mnexcm( \"SET ERR\", arglist, 1, ierflg )\n", + " \n", + " # Set starting values and step sizes for parameters\n", + " vstart = np.array([], dtype='float32')\n", + " step = np.array([], dtype='float32')\n", + " \n", + " gMinuit.mnparm( 0, \"\", vstart[0], step[0], lower_bound, upper_bound, ierflg)#par0\n", + " \n", + " #if needed\n", + " #gMinuit.FixParameter(0)\n", + " \n", + " #Minimization Step\n", + " \n", + " arglist[0] = 500\n", + " arglist[1] = 1e-20\n", + " gMinuit.mnexcm( \"MIGRAD\", arglist, 3, ierflg)\n", + " gMinuit.mnexcm(\"HESSE\", arglist, 4, ierflg)\n", + " \n", + " amin, edm, errdef, par0, err_par0 = map(ctypes.c_double, (0.18, 0.19, 0.20, 0.21, 0.22, 0.23))\n", + " nvpar, nparx, icstat = map(ctypes.c_int, (1983, 1984, 1985))\n", + " gMinuit.mnstat( amin, edm, errdef, nvpar, nparx, icstat )\n", + " gMinuit.mnprin( 4, amin.value )\n", + " gMinuit.GetParameter(0,par0,err_par0)\n", + " \n", + " fitpars = np.array([par0.value])\n", + " err_fitpars = np.array([err_par0.value])\n", + " ratio = chisq/(fit_bins-2)\n", + " \n", + " return firpars, err_fitpars, ratio, pulls" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "if __name__ == '__maun__':\n", + " \n", + " binwidth = (end-start)/nbins\n", + " \n", + " xup = np.array([start +1.0*(i+1.)*binwidth] for i in range(0,nbins))\n", + " xdown = np.array([start +1.0*(i+0.)*binwidth] for i in range(0,nbins))\n", + " \n", + " fitpars, err_fitpars, ratio, pulls = fit()\n", + " \n", + " #save data and plot here\n", + " " + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.6.6" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/bkg_reduction.ipynb b/bkg_reduction.ipynb new file mode 100644 index 0000000..e908b54 --- /dev/null +++ b/bkg_reduction.ipynb @@ -0,0 +1,264 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "import ROOT as r\n", + "import ctypes\n", + "import numpy as np\n", + "from array import array" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "r.TBrowser()" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "data_mumu = r.TFile(\"/disk/lhcb_data/davide/Rphipi/small/data/Ds_phipi_mumu/Ds_phipi_mumu_MagDown.root\")\n", + "MC_mumu = r.TFile(\"/disk/lhcb_data/davide/Rphipi/MC/Ds_phipi_mumu/Ds_phipi_mumu_MagDown.root\")" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "t_data_mumu = data_mumu.Get(\"Ds_OfflineTree/DecayTree\")\n", + "t_data_mumu" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "t_MC_mumu = MC_mumu.Get(\"Ds_OfflineTree/DecayTree\")\n", + "t_MC_mumu" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [], + "source": [ + "t_data_mumu.SetBranchStatus(\"*\",0)\n", + "t_data_mumu.SetBranchStatus(\"Ds_ENDVERTEX_CHI2\",1)\n", + "t_data_mumu.SetBranchStatus(\"Ds_ENDVERTEX_NDOF\",1)\n", + "t_data_mumu.SetBranchStatus(\"Ds_OWNPV_CHI2\",1)\n", + "t_data_mumu.SetBranchStatus(\"Ds_OWNPV_NDOF\",1)\n", + "t_data_mumu.SetBranchStatus(\"Ds_IP_OWNPV\",1)\n", + "t_data_mumu.SetBranchStatus(\"Ds_IPCHI2_OWNPV\",1)\n", + "t_data_mumu.SetBranchStatus(\"Ds_DIRA_OWNPV\",1)\n", + "t_data_mumu.SetBranchStatus(\"Ds_ConsD_M\",1)\n", + "\n", + "t_data_mumu.SetBranchStatus(\"Ds_Hlt1Phys_TOS\",1)\n", + "t_data_mumu.SetBranchStatus(\"Ds_Hlt1TrackMVADecision_TOS\",1)\n", + "#t_data_mumu.SetBranchStatus(\"Ds_Hlt2RareCharmD2PiMuMuOSDecision_TOS\",1)\n", + "\n", + "t_MC_mumu.SetBranchStatus(\"*\",0)\n", + "t_MC_mumu.SetBranchStatus(\"Ds_ENDVERTEX_CHI2\",1)\n", + "t_MC_mumu.SetBranchStatus(\"Ds_ENDVERTEX_NDOF\",1)\n", + "t_MC_mumu.SetBranchStatus(\"Ds_OWNPV_CHI2\",1)\n", + "t_MC_mumu.SetBranchStatus(\"Ds_OWNPV_NDOF\",1)\n", + "t_MC_mumu.SetBranchStatus(\"Ds_IP_OWNPV\",1)\n", + "t_MC_mumu.SetBranchStatus(\"Ds_IPCHI2_OWNPV\",1)\n", + "t_MC_mumu.SetBranchStatus(\"Ds_DIRA_OWNPV\",1)\n", + "t_MC_mumu.SetBranchStatus(\"Ds_ConsD_M\",1)\n", + "\n", + "t_MC_mumu.SetBranchStatus(\"Ds_Hlt1Phys_TOS\",1)\n", + "t_MC_mumu.SetBranchStatus(\"Ds_Hlt1TrackMVADecision_TOS\",1)\n", + "#t_MC_mumu.SetBranchStatus(\"Ds_Hlt2RareCharmD2PiMuMuOSDecision_TOS\",1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#end, start = fit region" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "def fcn( npar, gin, f, par, iflag):\n", + " \n", + " global ncount\n", + " global chisq\n", + " global pulls\n", + " global fitfunavg\n", + " global delta\n", + " \n", + " pull = np.array([0. for i in range(nbins)])\n", + " residuals = np.array([0. for i in range(nbins)])\n", + " \n", + " fitfunavg=np.array([func((xup[i]+xdown[i])/2) for i in range(nbins)])\n", + " \n", + " chisq, delta = 0., 0.\n", + " \n", + " delta = np.array([data[i]-fitfun[i] for i in range(0,nbins)])\n", + " \n", + " chisq = np.dot(np.transpose(delta),np.dot(Vxinv,delta))\n", + " \n", + " for i in range(nbins):\n", + " \n", + " pulls[i] = (y[i]-fitfun[i])/(np.sqrt(Vx[i,i]))\n", + " \n", + " print('X^2', chisq)\n", + " f[0] = chisq\n", + " ncount+=1" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "###DEFINE FIT FUNCTIONS HERE\n", + "###AS FUNC" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "##___MINUIT______\n", + "\n", + "def fit():\n", + " \n", + " gMinuit = r.TMinuit(5)\n", + " gMinuit.SetFCN(fcn)\n", + " \n", + " arglist = np.array([0. for i in range(10)], dtype='float32')\n", + " ierflg = ctypes.c_int(1982)\n", + " \n", + " arglist[0]=1\n", + " gMinuit.mnexcm( \"SET ERR\", arglist, 1, ierflg )\n", + " \n", + " # Set starting values and step sizes for parameters\n", + " vstart = np.array([], dtype='float32')\n", + " step = np.array([], dtype='float32')\n", + " \n", + " gMinuit.mnparm( 0, \"\", vstart[0], step[0], lower_bound, upper_bound, ierflg)#par0\n", + " \n", + " #if needed\n", + " #gMinuit.FixParameter(0)\n", + " \n", + " #Minimization Step\n", + " \n", + " arglist[0] = 500\n", + " arglist[1] = 1e-20\n", + " gMinuit.mnexcm( \"MIGRAD\", arglist, 3, ierflg)\n", + " gMinuit.mnexcm(\"HESSE\", arglist, 4, ierflg)\n", + " \n", + " amin, edm, errdef, par0, err_par0 = map(ctypes.c_double, (0.18, 0.19, 0.20, 0.21, 0.22, 0.23))\n", + " nvpar, nparx, icstat = map(ctypes.c_int, (1983, 1984, 1985))\n", + " gMinuit.mnstat( amin, edm, errdef, nvpar, nparx, icstat )\n", + " gMinuit.mnprin( 4, amin.value )\n", + " gMinuit.GetParameter(0,par0,err_par0)\n", + " \n", + " fitpars = np.array([par0.value])\n", + " err_fitpars = np.array([err_par0.value])\n", + " ratio = chisq/(fit_bins-2)\n", + " \n", + " return firpars, err_fitpars, ratio, pulls" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "if __name__ == '__maun__':\n", + " \n", + " binwidth = (end-start)/nbins\n", + " \n", + " xup = np.array([start +1.0*(i+1.)*binwidth] for i in range(0,nbins))\n", + " xdown = np.array([start +1.0*(i+0.)*binwidth] for i in range(0,nbins))\n", + " \n", + " fitpars, err_fitpars, ratio, pulls = fit()\n", + " \n", + " #save data and plot here\n", + " " + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.6.6" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}