diff --git a/fitter.ipynb b/fitter.ipynb new file mode 100644 index 0000000..ae8d8ee --- /dev/null +++ b/fitter.ipynb @@ -0,0 +1,212 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import ROOT as r\n", + "import pickle\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "l_index=0\n", + "l_flv=['e','mu']\n", + "\n", + "data = np.load('/disk/lhcb_data/davide/Rphipi/selected_data/'+l_flv[l_index]+l_flv[l_index]+'/sel_data_NN_'+l_flv[l_index]+l_flv[l_index]+'.npy')" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "#cut on mass\n", + "cut_indices=[]\n", + "\n", + "for i in range(len(data)):\n", + " if 1790" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "data_cut=data[cut_indices]/1000\n", + "plt.hist(data_cut,range=(1.79,2.03),bins=100);\n", + "fig=plt.gcf();\n", + "fig.set_size_inches(16,10)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "h1 = r.TH1D(\"h1\",\"data histogram\",100,1.79,2.03)\n", + "for i in range(len(data_cut)): h1.Fill(data_cut[i])" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "m = r.RooRealVar(\"m\",\"m\",1.8,2.5)\n", + "l = r.RooArgList(m)\n", + "data_hist = r.RooDataHist(\"data histogram\", \"data\", l, h1)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "#Creating the signal PDF\n", + "\n", + "\n", + "mean_dplus= r.RooRealVar(\"mean_D+\",\"mean_D+\",1.87,1.77,1.97)\n", + "sigma_dplus = r.RooRealVar(\"sigma_D+\",\"sigma_D+\",0.020,0.,0.2)\n", + "sig_dplus = r.RooGaussian(\"signal_D+\",\"signal_D+\", m, mean_dplus, sigma_dplus)\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "#Creating the signal PDF\n", + "\n", + "mean_ds= r.RooRealVar(\"mean_Ds\",\"mean_Ds\",1.97,1.87,2.07)\n", + "sigma_ds = r.RooRealVar(\"sigma_Ds\",\"sigma_Ds\",0.020,0.,0.2)\n", + "sig_ds = r.RooGaussian(\"signal_Ds\",\"signal_Ds\", m, mean_ds, sigma_ds)\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "coef0 = r.RooRealVar(\"c0\",\"coefficient #0\",1.0,-1.,1)\n", + "coef1 = r.RooRealVar(\"c1\",\"coefficient #1\",0.1,-1.,1)\n", + "coef2 = r.RooRealVar(\"c2\",\"coefficient #2\",-0.1,-1.,1)\n", + "bkg = r.RooChebychev(\"bkg\",\"background p.d.f.\",m,r.RooArgList(coef0,coef1,coef2))" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "#Add 2 sources of signal\n", + "\n", + "# first add the sig and the peak together with fraction fpeak\n", + "fpeak = r.RooRealVar(\"fpeak\",\"peaking signals fraction\",0.5,0.,1.)\n", + "# sigpeak(x) = fpeak*sig(x) + (1-fpeak)*bkg(x)\n", + "sigpeaks = r.RooAddPdf(\"sigpeak\",\"two mass peaks\",r.RooArgList(sig_ds,sig_dplus),r.RooArgList(fpeak))\n", + "\n", + "\n", + "# now we can add (sig+peak) to bkg with the fraction fpeak\n", + "fbkg = r.RooRealVar(\"fbkg\",\"background fraction\",0.5,0.,1.)\n", + "# model2(x) = fbkg*bkg(x) + (1-fbkg)*sigpeak(x)\n", + "model = r.RooAddPdf(\"model\",\"bkg+(sig_ds+sig_dplus)\",r.RooArgList(bkg,sigpeaks),r.RooArgList(fbkg))\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "##Plotting composite models\n", + "#\n", + "#xframe = m.frame()\n", + "#sigpeaks.plotOn(xframe)\n", + "#bkg_component = r.RooArgSet(bkg)\n", + "#model.plotOn(xframe,r.RooFit.Components(bkg_component),r.RooFit.LineStyle(2))\n", + "#\n", + "#xframe.Draw()" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [], + "source": [ + "model.fitTo(data_hist,r.RooFit.Extended(),r.RooFit.PrintLevel(-1))\n", + "\n", + "\n", + "xframe = m.frame(r.RooFit.Title(\"Fit to muon data\"))\n", + "sigpeaks.plotOn(xframe)\n", + "bkg_component = r.RooArgSet(bkg)\n", + "model.plotOn(xframe,r.RooFit.Components(bkg_component),r.RooFit.LineStyle(2))\n", + "data_hist.plotOn(xframe)\n", + "model.plotOn(xframe)\n", + "xframe.Draw()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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 +}