Navigation Menu

Skip to content

Instantly share code, notes, and snippets.

@wiso

wiso/gist:a7bdc5768ed6f8fce170

Created May 18, 2015 21:13
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save wiso/a7bdc5768ed6f8fce170 to your computer and use it in GitHub Desktop.
Save wiso/a7bdc5768ed6f8fce170 to your computer and use it in GitHub Desktop.
Download ZIP
discrete profiling
Raw
gistfile1.txt
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import ROOT\n",
"import numpy as np\n",
"from collections import OrderedDict\n",
"\n",
"import matplotlib as mpl\n",
"import matplotlib.pyplot as plt\n",
"\n",
"def safe_factory(func):\n",
" def wrapper(self, *args):\n",
" result = func(self, *args)\n",
" if not result:\n",
" raise ValueError('invalid factory input \"%s\"' % args)\n",
" return result\n",
" return wrapper\n",
"\n",
"ROOT.RooWorkspace.factory = safe_factory(ROOT.RooWorkspace.factory)\n",
"\n",
"import urllib\n",
"urllib.urlretrieve(\"https://gist.githubusercontent.com/mazurov/6194738/raw/67e851fdac969e670a11296642478f1801324b8d/rootnotes.py\",\n",
" \"rootnotes.py\")\n",
"#from rootnotes import default_canvas as TCanvas\n",
"from ROOT import TCanvas\n",
"\n",
"try:\n",
" ROOT.gROOT.ProcessLine(\".x ~/atlasstyle/AtlasStyle.C\")\n",
" ROOT.SetAtlasStyle()\n",
"except AttributeError:\n",
" print \"no ATLAS atyle for you\"\n",
"\n",
"#ROOT.gROOT.ProcessLine(\".x /home/turra/cms_classes/loadClasses.C\")\n",
"%matplotlib inline\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"colors = ROOT.kBlue, ROOT.kRed, ROOT.kOrange, ROOT.kViolet, ROOT.kGreen, ROOT.kYellow, ROOT.kSpring\n",
"from itertools import cycle\n",
"colors = cycle(colors)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"nsignal_expected = 1000\n",
"nbackground = 100000"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Define observable and POI\n",
"-------------------------"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"False"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ws = ROOT.RooWorkspace(\"workspace\")\n",
"ws_import = getattr(ws, \"import\")\n",
"\n",
"mass = ws.factory(\"mass[100,160]\")\n",
"mass.setBins(480)\n",
"mH = ws.factory(\"mH[125,110,150]\")\n",
"mH_shiftted = ws.factory(\"sum::mH_shifted(mH, mH_shift[-125])\")\n",
"mass_shifted = ws.factory(\"sum::mass_shifted(mass, mass_shift[-125])\")\n",
"bkg_index = ROOT.RooCategory(\"bkg_index\", \"bkg_index\")\n",
"ws_import(bkg_index)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Define signal\n",
"--------------"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"signal = ws.factory(\"Gaussian:signal(mass, mH, res[1.3])\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Define background\n",
"-----------------"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"bkgs = OrderedDict()\n",
"\n",
"bkgs['exp'] = 1, ws.factory(\"RooExponential::bkg_exp(mass, bkg_exp_par0[-0.05, -1, 0.01])\")\n",
"\n",
"bkgs['poly1'] = 1, ws.factory(\"RooBernstein::bkg_poly1(mass, {bkg_poly1_par0[1], bkg_poly1_par1[0.1, 0, 1000]})\")\n",
"bkgs['poly2'] = 2, ws.factory(\"RooBernstein::bkg_poly2(mass, {bkg_poly2_par0[1], bkg_poly2_par1[0.1,0,1000], bkg_poly2_par2[0,0,1000]})\")\n",
"bkgs['poly3'] = 3, ws.factory(\"RooBernstein::bkg_poly3(mass, {bkg_poly3_par0[1], bkg_poly3_par1[0.1,0,1000], bkg_poly3_par2[0,0,1000], bkg_poly3_par3[0,0,1000]})\")\n",
"bkgs['poly4'] = 4, ws.factory(\"RooBernstein::bkg_poly4(mass, {bkg_poly4_par0[1], bkg_poly4_par1[0.1,0,1000], bkg_poly4_par2[0,0,1000], bkg_poly4_par3[0,0,1000], bkg_poly4_par4[0,0,1000]})\")\n",
"bkgs['poly5'] = 5, ws.factory(\"RooBernstein::bkg_poly5(mass, {bkg_poly5_par0[1], bkg_poly5_par1[0.1,0,1000], bkg_poly5_par2[0,0,1000], bkg_poly5_par3[0,0,1000], bkg_poly5_par4[0,0,1000], bkg_poly5_par5[0,0,1000]})\")\n",
"\n",
"bkgs['landau'] = 2, ws.factory(\"Landau::bkg_landau(mass, bkg_landau_par0[80, 70, 130], bkg_landau_par1[10, 0, 1000])\")\n",
"\n",
"bkgs['exppol2'] = 2, ws.factory(\"EXPR::bkg_exppol2('exp(@1*(@0-100)/100.0 + @2*(@0-100)*(@0-100)/10000.0)', mass, bkg_exppol2_par0[-0.0025,-100.,100.],bkg_exppol2_par1[-0.02,-100,100.])\")\n",
"bkgs['exppol3'] = 3, ws.factory(\"EXPR::bkg_exppol3('exp(@1*(@0-100)/100.0 + @2*(@0-100)*(@0-100)/10000.0 + @3*(@0-100)*(@0-100)*(@0-100)/1000000.0)', mass, bkg_exppol3_par0[-0.0025,-100.,100.],bkg_exppol3_par1[-0.02,-100,100.],bkg_exppol3_par2[-0.02,-100,100.])\")\n",
"bkgs['exppol4'] = 4, ws.factory(\"EXPR::bkg_exppol4('exp(@1*(@0-100)/100.0 + @2*(@0-100)*(@0-100)/10000.0 + @3*(@0-100)*(@0-100)*(@0-100)/1000000.0 + @4*(@0-100)*(@0-100)*(@0-100)*(@0-100)/100000000.0)', mass, bkg_exppol4_par0[-0.0025,-100.,100.],bkg_exppol4_par1[-0.02,-100,100.],bkg_exppol4_par2[-0.02,-100,100.], bkg_exppol4_par3[-0.02,-100,100.])\")\n",
"\n",
"bkgs['laurent2'] = 2, ws.factory(\"EXPR::bkg_laurent2('@1 + @2/@0 + @3/(@0 * @0)', mass, bkg_laurent2_par0[0.002, -100, 100], bkg_laurent2_par1[-0.3, -100, 100], bkg_laurent2_par3[60, -100, 100])\")\n",
"bkgs['laurent3'] = 3, ws.factory(\"EXPR::bkg_laurent3('@1 + @2/@0 + @3/(@0 * @0) + @4/(@0 * @0 * @0)', mass, bkg_laurent3_par0[0.002, -100, 100], bkg_laurent3_par1[-0.3, -100, 100], bkg_laurent3_par3[60, -100, 100], bkg_laurent3_par4[1, -1000, 1000])\")\n",
"\n",
"bkgs['laurent_half'] = 4, ws.factory(\"EXPR::bkg_laurent_half('@1 + @2/@0 + @3/(@0 * @0) + @4/(sqrt(@0)) + @5 * @0', mass, bkg_laurent_half_par0[0.002, -100, 100], bkg_laurent_half_par1[-0.3, -100, 100], bkg_laurent_half_par3[60, -100, 100], bkg_laurent_half_par4[1, -1000, 1000], bkg_laurent_half_par5[1, -1000, 1000])\")\n",
"\n",
"ws.factory(\"RooBernstein::bkg_poly5_clone(mass, {bkg_poly5_par0_clone[1], bkg_poly5_par1_clone[0.5,0,1000], bkg_poly5_par2_clone[0.4,0,1000], bkg_poly5_par3_clone[0.3,0,1000], bkg_poly5_par4_clone[0.2,0,1000], bkg_poly5_par5_clone[0.2,0,1000]})\")\n",
"ws.factory(\"EXPR::bkg_exppol3_clone('exp(@1*(@0-100)/100.0 + @2*(@0-100)*(@0-100)/10000.0 + @3*(@0-100)*(@0-100)*(@0-100)/1000000.0)', mass, bkg_exppol3_par0_clone[-0.0025,-100.,100.],bkg_exppol3_par1_clone[-0.02,-100,100.],bkg_exppol3_par2_clone[-0.02,-100,100.])\")\n",
"ws.factory(\"Landau::bkg_landau_clone(mass, bkg_landau_clone_par0[80, 70, 130], bkg_landau_clone_par1[10, 0, 1000])\")\n",
"ws.factory(\"SUM::bkg_complicated(nc1[0.4] * bkg_poly5_clone, nc2[0.4] * bkg_landau_clone, bkg_exppol3_clone)\")\n",
"\n",
"bkgs['complicated'] = 999, ws.pdf(\"bkg_complicated\")\n",
"\n",
"#bkg_multi = ROOT.RooMultiPdf(\"bkg_multi\", \"bkg_multi\", bkg_index, bkg_list)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"mu = ws.factory(\"mu[1, 0, 3]\")\n",
"nsignal = ws.factory(\"prod::nsignal(mu, nsignal_expected[%d])\" % nsignal_expected)\n",
"\n",
"models = OrderedDict()\n",
"models['exp'] = ws.factory(\"SUM::model_exp(nbkg[%d, 0, 1000000] * bkg_exp, nsignal * signal)\" % nbackground)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def create_model_from_template(template, new_name, old_bkg, new_bkg):\n",
" command = \"EDIT::%s(%s, %s=%s)\" % (new_name, template, old_bkg, new_bkg)\n",
" return ws.factory(command)\n",
"\n",
"def create_model_from_exp(new_name, new_bkg):\n",
" return create_model_from_template('model_exp', new_name, 'bkg_exp', new_bkg)\n",
"\n",
"cmfe = create_model_from_exp\n",
"\n",
"for bkg_name, bkg in bkgs.iteritems():\n",
" models[bkg_name] = cmfe('model_%s' % bkg_name, bkg[1].GetName())"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1 exp 1\n",
"2 poly1 1\n",
"3 poly2 2\n",
"4 poly3 3\n",
"5 poly4 4\n",
"6 poly5 5\n",
"7 landau 2\n",
"8 exppol2 2\n",
"9 exppol3 3\n",
"10 exppol4 4\n",
"11 laurent2 2\n",
"12 laurent3 3\n",
"13 laurent_half 4\n",
"14 complicated 999\n"
]
}
],
"source": [
"for i, model in enumerate(models, 1):\n",
" print \"{:<5d}{:<30}{:d}\".format(i, model, bkgs[model][0])"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# fit the constant to landau\n",
"obs_set = ROOT.RooArgSet(mass)\n",
"data_landau = models['landau'].generateBinned(obs_set, ROOT.RooFit.ExpectedData())\n",
"for k, pdf in models.iteritems():\n",
" pdf.fitTo(data_landau)\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"False"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mu.setVal(1)\n",
"mH.setVal(125)\n",
"true_model = models['complicated']\n",
"ws.var('bkg_poly5_par5_clone').setVal(ws.var('bkg_poly5_par5_clone').getVal() * 1.5)\n",
"ws.var('bkg_poly5_par4_clone').setVal(ws.var('bkg_poly5_par4_clone').getVal() * 1.5)\n",
"ws.var('bkg_poly5_par3_clone').setVal(ws.var('bkg_poly5_par3_clone').getVal() * 1.5)\n",
"data = true_model.generateBinned(obs_set,\n",
" ROOT.RooFit.ExpectedData()\n",
")\n",
"ws_import(data)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"bkg only fit\n",
"------------"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"<ROOT.TCanvas object (\"c1\") at 0x4d26f30>"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ws.var('mu').setVal(0)\n",
"ws.var('mu').setConstant()\n",
"\n",
"canvas = TCanvas()\n",
"frame = mass.frame()\n",
"\n",
"legend = ROOT.TLegend(0.4, 0.7, 0.9, 0.9)\n",
"legend.SetNColumns(2)\n",
"\n",
"data_bkg_only = true_model.generateBinned(obs_set,\n",
" ROOT.RooFit.ExpectedData()\n",
")\n",
"data_bkg_only.plotOn(frame)\n",
"\n",
"\n",
"for color, pdfname in zip(colors, models):\n",
" pdf = ws.pdf('model_' + pdfname)\n",
" pdf.fitTo(data_bkg_only)\n",
" pdf.plotOn(frame, ROOT.RooFit.LineColor(color), ROOT.RooFit.Name(\"curve_sb_%s\" % pdfname))\n",
" curve = frame.findObject(\"curve_sb_%s\" % pdfname)\n",
" curve.SetFillColor(0)\n",
" curve.SetMarkerSize(0)\n",
" legend.AddEntry(curve, pdfname)\n",
" \n",
"frame.addObject(legend)\n",
"frame.SetTitle(\"bkg-only fit\")\n",
"frame.Draw()\n",
"canvas"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"signal + background fit\n",
"------------------------"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false,
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" ndof bkg min nll chi/ndof mH mu \n",
"================================================================================\n",
"exp 1 -659715 0.295886 125.02 ± 0.22 0.577 ± 0.10 \n",
"poly1 1 -659005 3.156457 124.66 ± 26.66 0.000 ± 1.88 \n",
"poly2 2 -659736 0.209648 125.08 ± 0.20 0.731 ± 0.10 \n",
"poly3 3 -659785 0.004059 125.02 ± 0.17 1.011 ± 0.10 \n",
"poly4 4 -659786 0.000141 125.00 ± 0.17 0.992 ± 0.10 \n",
"poly5 5 -659786 0.000012 125.00 ± 0.17 1.000 ± 0.11 \n",
"exppol2 2 -659784 0.006975 125.02 ± 0.17 0.951 ± 0.10 \n",
"exppol3 3 -659786 0.002118 125.00 ± 0.17 0.961 ± 0.10 \n",
"exppol4 4 -659786 0.000216 125.00 ± 0.17 0.991 ± 0.10 \n",
"laurent2 2 -659782 0.016514 125.03 ± 0.17 0.927 ± 0.10 \n",
"laurent3 3 -659783 0.013812 125.02 ± 0.17 0.928 ± 0.10 \n",
"laurent_half 4 -659766 0.082723 125.06 ± 0.18 0.844 ± 0.10 \n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAUsAAAILCAYAAABhB75vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH8BJREFUeJzt3XuUZWV55/HvQzdBkZukIyqCaMRBjCgqF6+UYrRxDKxo\nFG11xhthZqKOE2fGWzJ0NNGYrKCjTkxjAGOiqDFqIEG8REqRhVyUqzQEAiwbMSgK3iBKh2f+2Lug\nOH3q1Ft99t7n7FPfz1q16lx2P89bVV2/2pd37x2ZiSRptB0mPQBJ6gPDUpIKGJaSVMCwlKQChqUk\nFTAsJanA2kkPYHtEhPOdJLUiM2OpN3r3UQ27zfontFyfjR18j1rtYf3J9+h7/Wn8GkZli5vhklTA\nsJSkAoblUHNtN5hvu0EHPaw/+R59r99Fj8bqR72d3isRkbnUTthG6pOZtFZf0nQalS2uWUpSAcNS\nkgoYlpJUwLCUpAKGpSQVMCwlqYBhKUkFDEtJKmBYSlIBw1KSChiWklTAsJSkAoalJBUwLCWpgGEp\nSQUMS0kqYFhKUgHDUpIKGJaSVMCwlKQChqUkFTAsJamAYSlJBQxLSSpgWEpSAcNSkgoYlpJUwLCU\npAKGpSQVMCwlqYBhKUkFDEtJKmBYSlKBTsMyIk6JiJsj4vIRy8xFxMURcUVEzHc4PElaUmRmd80i\nngb8FPhIZj5myPt7AOcCz8nMGyNiXWbeMmS5zMxob5xkJq3VlzSdRmVLp2uWmXkOcOuIRTYAf5eZ\nN9bLbxOUkjQJ07bPcn9gz4g4OyIuioiXT3pAkgSwdtIDGLAj8HjgSGBn4LyI+HpmXjPZYUla7aYt\nLLcAt2TmHcAdEfFV4LHANmEZERsXPZ3PzPlORihpZkTEHDBXtGyXB3gAImI/4IwlDvAcAHwAeA6w\nE3A+cGxmXjmwnAd4JDVuVLZ0umYZEacBRwDrImILcALVpjeZuSkzr4qIs4DLgLuADw0GpSRNQudr\nlk1wzVJSG6Zm6pAk9ZVhKUkFDEtJKmBYSlIBw1KSChiWklTAsJSkAoalJBUwLCWpgGEpSQUMS0kq\nYFhKUgHDUpIKGJaSVMCwlKQChqUkFTAsJamAYSlJBQxLSSpgWEpSAcNSkgoYlpJUwLCUpAKGpSQV\nMCwlqYBhKUkFDEtJKmBYSlIBw1KSChiWklTAsJSkAoalJBUwLCWpgGEpSQUMS0kqYFhKUgHDUpIK\nGJaSVMCwlKQChqUkFTAsJamAYSlJBQxLSSpgWEpSAcNSkgoYlpJUwLCUpAKGpSQVMCwlqYBhKUkF\nDEtJKtBpWEbEKRFxc0Rcvsxyh0TE1oh4fldjk6RRul6zPBVYP2qBiFgDvBs4C4guBiVJy+k0LDPz\nHODWZRZ7HfAp4Pvtj0iSykzVPsuI2Bs4Bvhg/VJOcDiSdLe1kx7AgPcCb87MjIhgxGZ4RGxc9HQ+\nM+dbHpukGRMRc8Bc0bKZ3a68RcR+wBmZ+Zgh713HPQG5DrgdOC4zTx9YLjOztf2ZEWSm+0ul1WZU\ntkzVmmVmPnzhcUScShWqp4/4J5LUiU7DMiJOA44A1kXEFuAEYEeAzNzU5VgkaSU63wxvgpvhktow\nKlum6mi4JE0rw1KSChiWklTAsJSkAoalJBUwLCWpgGEpSQUMS0kqYFhKUgHDUpIKGJaSVMCwlKQC\nhqUkFTAsJamAYSlJBQxLSSpgWEpSAcNSkgoYlpJUwLCUpAKGpSQVMCwlqYBhKUkFDEtJKmBYSlIB\nw1KSChiWklTAsJSkAoalJBUwLCWpgGEpSQUMS0kqYFhKUgHDUpIKGJaSVMCwlKQChqUkFTAsJamA\nYSlJBQxLSSpgWEpSAcNSkgoYlpJUwLCUpAKGpSQVMCwlqYBhKUkFDEtJKmBYSlIBw1KSChiWklSg\n07CMiFMi4uaIuHyJ918aEZdGxGURcW5EHNTl+CRpKV2vWZ4KrB/x/nXA0zPzIOAdwEmdjEqSltFp\nWGbmOcCtI94/LzN/VD89H3hIJwOTpGVM8z7LVwNnTnoQkgSwdtIDGCYingG8CnjKiGU2Lno6n5nz\nLQ9L0oyJiDlgrmjZzGx1MNs0jNgPOCMzH7PE+wcBnwbWZ+a1SyyTmRntjZHMpLX6kqbTqGyZqs3w\niNiXKihftlRQStIkdLpmGRGnAUcA64CbgROAHQEyc1NE/CXwm8C3639yZ2YeOqSOa5aSGjcqWzrf\nDG+CYSmpDb3ZDJekaWVYSlIBw1KSChiWklTAsJSkAoalJBUwLCWpgGEpSQUMS0kqYFhKUgHDUpIK\nGJaSVMCwlKQChqUkFTAsJamAYSlJBQxLSSpgWEpSAcNSkgoYlpJUwLCUpAJrJz0ANSeCOWCufjoH\nzNeP5zPvfixpO3gr3KH1+38r3D5+DYa9Js37hq+4fv+CZlDfv4a+j1/95H3DJWlMhqUkFfAAj4q5\nT1Grmfssh9bv//6ytr+GvtfvK/9gtcsDPCuu3/9f1L6H2Sz8DNrm96h5HuCRpDEZlpJUwLCUpAKG\npSQVMCwlqYBhKUkFDEtJKuAZPFKDnDQ+u5yUPrR+/yf79n3SuD+DyddfjZyULkljcjNcq4abyBqH\nm+FD6/d/86bvm4B9r99Fj1n4fzpt3AyXpDEZlpJUwLCUpAKGpSQVMCwlqYBhKUkFnGcp6W5dzEXt\n63xX51kOrd//+Wt9n+PX9/pd9Oh7/a56rITzLCVpTJ2GZUScEhE3R8TlI5Z5X0RcExGXRsTBXY5v\nVkRwUv35zAj2mPR4pFnQ9ZrlqcD6pd6MiOcCj8jM/YHfBj7Y1cBmzCPrz0cBmyY5EGlWdBqWmXkO\ncOuIRY4G/qpe9nxgj4jYq4uxzZjb688XAsdPciDSrCgOy4jYd9Hjp0fEk1oYz97AlkXPbwQe0kKf\nWbeh/vzsTG6b6EikGbGSqUPHR8TjgV8AlwA7Aee1MKbBI1H9O1w/YZncFlF9nvRYtsfifa7Ahr5+\nHZotI8MyInYAXpyZH8vMt9Wv7QQcBjy0hfF8B9hn0fOH1K8NG9s898zP2g+4ITM31u9tBNje57CR\niD/Y2FS9yTw/AdhIc/W6qw/XPBv2h3qfa0Rs7tP4F54v/J23/vQ+j4g5Fv4jsMwcz8xc8gN4A9UB\nF4CXA7+26L3/OOrfjqi5H3D5Eu89Fzizfnw48PUllsvt6V0+xmy1fhcfbX4NkCdBJuSZkHu0UP/M\nuv4FbdTv6mfc55/BLHyPtm88LDme5TbD3w+8CLgW+AHwyoh4DLAzsHtE/BQ4LzN/sUwdACLiNOAI\nYF1EbKH6875jPcJNmXlmRDw3Iq4Ffga8sqSuOjd4tP3YhutvoDoQ6D7XpbX9M9CA7TqDJyJ2ptoU\nfwrw8Mx8VdMDW6Z/pmfwjNTm11DvSzyK6mh7K4Hm2SnL1u79z6CrHisxKls83XFo/en6AW6Pln9R\n96Ba87t/W2t+huWytXv/M+iqx0p4uqMatfDL6Sby5Pgz6J5hKUkFDEtJKmBYSlIBw1KSChiWUgu8\nTN7sMSyldniZvBljWErt8DJ5M8ZJ6UPrT9dE2e3R90ndM1C/95PGnZR+b65ZatXpYn+ik8Znj2Gp\n1cj9iVoxw1KrkfsTtWLusxxaf7r2o2yPvu/P6vtFKOo+vf0edVG/rR4RzAFz9dM57rmo73zm6Av8\netWhFdc3LK0//T36Xr+LHiut7wEeSRqTYSlpIvp2lpNhKWkbHQVZr2YlGJaShukiyHo1K8GwlDRM\nF0G2of7cixvTeTR8aH2Phlt/+ns4var5+h4Nl7Qinq65LcNSkgqsnfQA1JyBMxe+EsHG+vGyZy5I\nGs19lkPr93+fZdumbV/TtNXvokff63fRw32WktQxw1KSChiWklTAsJSkAh7gGVrfAzzDjHOdwGmo\nP9Crlwcv/B61W9/rWa64vmE562YhCNo2C98jj4ZLUscMS0kqYFhKUgHDUpIKGJaSVMCwlKQChqUk\nFXCe5dD6/Z4fp+W19TPuctJ425xnOfCeYTmsvmE56/wZL8+wvDc3wyWpgFdK16rhleRXl8X3Pgc2\njHs/ITfDh9Z3E03q+2Z4BPPAEfXTT2Zy7PL/xs1wSatPo/c+d81yaH3XLKUZWLNc8b3PXbOUtOo0\nfe9zD/B0ZJbm30mrkZvhQ+tP19wvaRJmYeK+k9INS6kVs7YFZFgalpIKeAaPJHWs07CMiPURcVVE\nXBMRbxry/rqIOCsiLomIKyLiFV2OT5KW0tlmeESsAa4GngV8h2qi6Esyc/OiZTYCO2XmWyJiXb38\nXpm5daCWm+GSltXXzfBDgWsz84bMvBP4OHDMwDLfBXarH+8G/GAwKCVpErqcZ7k3sGXR8xuBwwaW\n+RDw5Yi4CdgVeFFHY5OkkbpcsyzZ3n8rcElmPhh4HPD/ImLXdoclScvrcs3yO8A+i57vQ7V2udiT\ngT8CyMx/iYjrgf8AXDRYrN6/uWA+M+ebHKyk2RcRc9wzr3T0sh0e4FlLdcDmSOAm4AK2PcBzIvCj\nzPyDiNgL+AZwUGb+cKCWB3gkLavJAzydrVlm5taIeC3weWANcHJmbo6I4+v3NwHvBE6NiEupdhH8\n78GglKRJ8AyeofVds5RmQV+nDklSbxmWklTAsJSkAoalJBUwLCWpgGEpSQUMS0kqYFhKUgHDUpIK\nGJaSVMCwlKQCnhs+tL7nhkt9Nc7tfL0V7orrG5bSauSFNCRpTIalJBUwLDsWwUn15zMj2GPS45FU\nxrDs3iPrz0cBmyY5EEnlDMvu3V5/vhA4fpIDkVTOo+FD67d3tLre9L4VuH8mt7XRQ9L2cerQius7\ndUhajZw6JEljMiwlqYBhKUkFDEtJKmBYSlIBw1KSChiWklTAsJSkAoalJBUwLCWpgGEpSQUMS0kq\nYFhKUoG1kx7A9opgY/1w2Tu2SdK4vETb0Ppeok1ajbxEmySNybCUpAKGpSQVMCwlqYBhKUkFDEtJ\nKmBYSlIBw1KSChiWklTAsJSkAoalJBXw3PCh9Zs/dzuCOWCufjoHd1/8wwuBSFNiVLYYlkPre6EL\naTXyQhqSNCbDUpIKGJaSVKDTsIyI9RFxVURcExFvWmKZuYi4OCKuiIj5LscnSUvp7ABPRKwBrgae\nBXwHuBB4SWZuXrTMHsC5wHMy88aIWJeZtwyp5QEeSY2blgM8hwLXZuYNmXkn8HHgmIFlNgB/l5k3\nAgwLSkmahC7Dcm9gy6LnN9avLbY/sGdEnB0RF0XEyzsbnSSN0OXdHUu293cEHg8cCewMnBcRX8/M\nawYXjIiNi57OZ+Z8E4OUtHpExBz3nCwyUpdh+R1gn0XP96Fau1xsC3BLZt4B3BERXwUeC2wTlpm5\nsaVxSlol6pWs+YXnEXHCUst2uRl+EbB/ROwXEb8EHAucPrDM3wNPjYg1EbEzcBhwZYdjlKShOluz\nzMytEfFa4PPAGuDkzNwcEcfX72/KzKsi4izgMuAu4EOZaVhKmjjPDR9a36lD0mo0LVOHJKm3ehuW\nEZwZwR6THoek1aG3YQkcBWya9CAkrQ59DssLgeMnPQhJq0NvD/BA3j+T29qp7wEeaTWayQM8bQWl\nJA3T27CUpC4ZlpJUwLCUpAKGpSQVMCwlqYBhKUkFDEtJKmBYSlIBw1KSChiWklTAsJSkAoalJBUw\nLCWpgGEpSQUMS0kqYFhKUgHDUpIKGJaSVMCwlKQChqUkFTAsJamAYSlJBQxLSSpgWEpSAcNSkgqs\nnfQAtlcEG+uH85nMT3AoklaByMxJj2HFIiIzM9qrT2bSWn1J02lUtrgZLkkFDEtJKmBYSlIBw1KS\nChiWklTAsJSkAoalJBUwLCWpgGEpSQUMS0kqYFhKUgHDUpIKGJaSVMCwlKQChqUkFTAsJamAYSlJ\nBQxLSSrQaVhGxPqIuCoiromIN41Y7pCI2BoRz+9yfJK0lM7CMiLWAB8A1gMHAi+JiEctsdy7gbPA\n++BImg5drlkeClybmTdk5p3Ax4Fjhiz3OuBTwPc7HJskjdRlWO4NbFn0/Mb6tbtFxN5UAfrB+qX+\n3XpS0kzq8r7hJcH3XuDNmZkREYzYDI+IjYuezmfm/HjDk7TaRMQcMFe0bFf3DY+Iw4GNmbm+fv4W\n4K7MfPeiZa7jnoBcB9wOHJeZpw/U8r7hkho3Klu6DMu1wNXAkcBNwAXASzJz8xLLnwqckZmfHvKe\nYSmpcaOypbPN8MzcGhGvBT4PrAFOzszNEXF8/f6mrsYiSSvV2Zplk1yzlNSGUdniGTySVMCwlKQC\nhqUkFTAsJamAYSlJBQxLSSpgWEpSAcNSkgoYlpJUwLCUpAKGpSQVMCwlqYBhKUkFDEtJKmBYSlIB\nw1KSChiWklTAsJSkAoalJBUwLCWpgGEpSQUMS0kqYFhKUgHDUpIKGJaSVMCwlKQChqUkFTAsJamA\nYSlJBQxLSSpgWEpSAcNSkgoYlpJUwLCUpAKGpSQVMCyHmm+1ekTMtdqggx7Wn3yPvtfvokeT9Q3L\noebbbjDXdoMOelh/8j36Xr+LHo3VNywlqYBhKUkFIjMnPYYVi4j+DVpSL2RmDHu9l2EpSV1zM1yS\nChiWklTAsJSkAoblEiLilQ3VeVREHBkRuwy8vr6J+nWtp0bEgfXjuYj4nxFxZFP1h/T7SIu1nxYR\nb4yIZ7fVQ9MhIn65wVp/XX9+Q1M1t+nRhwM8ERHA84GnAgmcA3w2Wxx8RGzJzH3GrPF64HeAzcDB\nwH/PzM/W712cmQc3MM53Ac8A1gBnA08H/hH4deCMzPzTMeufQfU9X3yE8JnAl4HMzKPHrH9BZh5a\nPz6O6vv1GeDZwD9k5rvGrL8WeA3wEOBzmXnuovd+LzP/cJz6dZ3dgdcBtwB/CbwNOAS4GHhnZt4x\nbo8hPf85Mx/ZUK2DgJOovkdnAm/KzFvr9+7++YzZ45l1j1uA1wN/Dayt335xZl44Zv0rgWcBZzFk\nInpm/nCc+tCfsPwg8KvAaVS/tC8CrsvM/zZm3ctHvP3IzNxpzPpXAIdn5k8jYj/gU8DfZOZ7GwzL\nK4GDgF8CbgYekpk/ioj7Audn5kFj1r8YuJIqBO6i+v6fBrwYIDO/Mm79he9DRFwEHJWZ34+I+9Xj\n/7Ux658M3Be4EHgZ8JXM/N3B3mP2+AxwXd3nIOAy4OPA0cCemfmqMev/hG3/YO0M3E71B2u3Meuf\nC7wDOB94NfAq4OjMvLbB79E3gFcAuwCfA34jM8+JiMcD/zcznzZm/dcD/xV4OHDTwNuZmQ8fp/5C\nlan/AK4Cdlj0fAfgqgbq3ky1xrffkI+bGqj/rYHnuwCfB94DXNLQ9+aSYY+HPd/O+muA3wW+BBxc\nv3Z9gz/by4A9gV8GLm5h/Jcverwj8CHg08B9BvuN8zXUnwP414X/q/Xzyxqo/z7gI8ADF9Vt9Gcw\n8PwZwLXA4Q1+jy5e9HjzUu810Ocvmqo1+LGwGjztrgX2BW6on+9bvzaufwR2ycyLB9+IiLHWmGrf\ni4jHZeYlAFmtYT4POJlqDaQJP4+InTPzduAJCy9GxB5Ua4Jjycx/B06MiE8C74mI70Gj/292A76x\n0C4iHpSZ342IXRuqv+PCg8y8EzguIk4A/onqj1cT7qrrZ0R8LjMXPx+7eGa+PiKeCHwsIv4e+MDY\nRQdaRMTumfmjut/ZEfF8qj8q92+ox+LjI29ZeFDvYttx28VXJiL2rB++bdHju+Uq2gz/KtU+oAuo\nNkcOpdqs+jEN7DdrS0TsA9yZmf868HoAT8nMrzXQ4z6Z+W9DXl8HPCgzR+1q2J5+zwOenJlvbbLu\nkD47A3tl5vVj1vko1a6Pzw28/hrgg5nZxC/qycAbMvMnA68/AvhwZj513B51vTXAa4HfAh6RmQ9q\nqO5LqXZrnTfw+r7A72fmcQ30OAb4Umb+bOD1XwVekJl/Mmb9G6iyYajMfNg49aE/YTk35OWFfTiZ\n4+83OxE4OTO/NU6dSdXvokff609KREQ2/EsWEQ8GHpeZZzZc94WZ+bcDr70oMz/Zpx5t6UtYHpiZ\nVw68NpeZ8w3VP45q5/OOwCnAaQubJH2o30WPvtaPiBew7cGRBZmZn572HgP1F35hF3o18jXUfbY5\nmNPUAZ6Oe9wf2J9qvzQAmfnVsev2JCyvoJpq8CdURxzfDRySmYc33OcAql/YDcDXgA9l5tl9qd9F\nj77Vj4gPM3rzbOz5tG336KD+UcBzgWOpjuIvBPGuwIHZzNSh1nvUfY6jmpq0D9XUrcOB8zLzmWPX\n7klY3o8qIJ9ItVP+Y8AfL+xIb6jHGuA3gFdSzTf7JNW8ztsz89hpr99Fj77X13AR8ViqWSFvB36f\ne4Lsx8DZWc+5nPYedZ8rqI5vnJeZj6v/+L4rM39z7OJtHWZv8gPYCfhT4FKqo+Avbrj+e+q6JwGH\nDrx39bTXn4WvoYP6e9Q9vlF//Bmwe8P/j1rt0UH9HZv8fkyiB3BR/fkS4D714yubqN2X0x0vAP6N\nas3yacCGiPjb0f9kRS4DHpuZv52ZFwy8d1gP6nfRo+/1T6Fai3kh1UkNPwFObaBulz3arn9YRHwx\nIq6JiOvrj+sarN9Fjy31PsvPAl+MiNO5Z8rhWPqyGX4Y8EjgYZn59oh4KPCfMvMdY9Z9AtvuOGfh\neWZ+c5rrd9Gj7/UX9bk0Mx+73GvT3KOD+lcDbwC+Cfz7wuuZeUsT9bvqsajXHNU83rMy8xfj1uvL\npPRXUn1jj6Ta5/ET4BiqU7TG8WeM2HFOdSbDNNfvokff6y+4IyKelpnnAETEU6lOF2xS2z3arn9b\nDsxHbUHrPep933tRnYIawAOBb49dtydrlhdn5sFx7/OIG10r0GyLiMdRnTK4e/3SrcB/zsxL+9Kj\ng/p/THV666eBny+83tTafRc9IuJ1wAnA97j3mutjxq7dk7A8H3gy1c7bgyPiV4AvZENzsyLil6hO\nwn96/dI81Tmmd/ahfhc9+l5/UZ/dADLzx03W7bJHW/UjYp4ha/mZ2dTafes9IuJfqA4Q/qCJeveq\n3ZOwfBnVDu0nAH9FdbrX72VDs/7r09XW1rUDeDmwNTNf04f6XfSYgfrrqNY4Fl/m7+1N/lK13aOL\nr6HvIuJs4NlN/5GFnoQlQEQ8imqfJcA/ZebmBmtflgOXMhv22rTW76LHDNT/EvAV4G+owngDMJeZ\nz2qifhc9Oqj/QOCPgL0zc31UF5R+Umae3ET9NntExBvrhwcCBwD/ACwc1MnMPHGc+tCjK6Vn5ubM\n/ED90VhQ1rZGddEDgIWT+7f2qH4XPfpe/4GZ+Y7MvD4zr8vqor97NVi/ix5t1/8w8AXgwfXza4D/\n0WD9NnvsSnXCyreBL1Jd33WX+qORK1j15Wh42/4X8OV6vldQXc+ykdtKdFS/ix59r/+FiHgJ8In6\n+Qupfmmb1HaPtuuvy8xPRMSbobqkXUQ0/Ue9lR6ZubFkuYh4f2a+bnt69GYzvE1RXVX8jVS3S7gN\nuAg4MYdc+mwa63fRYwbq/5Tq6uILp8juACxcLixzzKuNd9Gjg/rzwAuoLqV2cEQcDrw7M48Yp27X\nPZbpv90X7TAsgfpsoB9z731Bu2fmC/tQv4sefa+v5dUnCLwfeDTwLeBXgN9qeHpV6z2W6b/dYelm\neOXRmXngoudfjureNn2p30WPXtePuNdN7+4CvpaZn2mqfhc92qxfT+R+ev1xANUfrKubOPOlyx5t\n6s0BnpZ9MyKetPCk3jT4xojlp61+Fz36Xv/PgeOpzkH/FvBfIuLPG6zfRY/W6md1+5ANmbk1M6/I\nzMubDrEuerTJzXAgIq6iOvd8C9X8tX2Bq6mOxua401fart9Fjxmpf2DWl/WLiB2orkZzwDh1u+zR\nQf33UF18+RNU+0IbPT+/ix4x/Ersd78WEa/IzA9vT203wyvre16/ix59r9/WTe+67NF2/YOp/lC9\nfeD1xs7g6aDHW4HBK5Ld/dr2BiW4ZqlVIjq46V3bPbr4GvoqOrgSu2uWWi3+z5DXhl0abpp7tFo/\nqlsEb1MvMwfXAqexx01U+7iPqT8vvhJ7IxPrDUutFt/LFm9611GPtuv/jHsC7L7A84CmZ2200qOe\nenRpRHw0WzgvHNwM1yoRHdz0ru0eXXwNA/12orq6V2sTxpvuEdU1Pk+gOgNsYWUwM/Ph49Z26pBW\ni8Oo7vh3HtU+v+9SXfavTz26+BoWux+wd4v12+hxMnAi1VzUQ+qPRu4c6Wa4VoutwB1Ua2T3Aa7L\nBu8O2lGPVutHxOWLnu4APIBtj1pPe4/WrsTuZrhWhYi4FDid6hdzHbAJ+HnDp5y22qOD+vsteroV\nuLnp/X9t94gWr8RuWGpViJZuetdljy6+hrrPA6jWXAHIzLHvX9NVj2jxSuyGpVaFiPgL6pveZeYB\nEbEn1YGFJ/alRwf1j6a6gdyDqe5h81Bgc2Y+uon6XfVoiwd4tFoclpm/Q7XPj8z8IdVpd33q0Xb9\nPwSeBPxzZj6M6s4E5zdYv/UeEfHAiDg5Is6qnx8YEa9uorZhqdXiF/VVbwCI6qZ3TR/gabtH2/Xv\nzOr+3TtExJrMPBtobM27ox4fpqWrvRuWWi3eD3wGeEBEvBM4F3hXz3q0Xf/WiNiV6kZoH42I9wE/\nbbB+Fz3WZeYnqG+DWx88auRq7+6z1KoRLd70rqsebdaPiF2oNvF3AF4K7AZ8NJu9A2arPdq8Erth\nKWlmtHkldielS6tcfW+fpdaaMpu7P1HbPVq9ErtrlpJmRkRcmJmHtFLbsJQ0K9q8ErthKWlmeAaP\nJE2YB3gkzYw2r/ZuWEqaJa1d7d3NcEkzq8krsXu6o6RZ1tiV2N0MlzQz2rwSu5vhkmZGm1diNywl\nzZw2rsTuPktJMyMijo6Ia4Drga8ANwCN3MDMsJQ0S1q7ErthKWmWtHYldo+GS5olg1di/x4NXYnd\nAzySZkabV2I3LCWpgJvhknqvkyuxu2YpScvzaLgkFTAsJamAYSlJBQxLSSpgWEpSgf8PIXHFDDRG\nNfQAAAAASUVORK5CYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x1789b90>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"<ROOT.TCanvas object (\"c1\") at 0x4d94270>"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mu.setConstant(False)\n",
"mH.setConstant(False)\n",
"mu.setVal(1)\n",
"mH.setVal(125)\n",
"\n",
"canvas = TCanvas()\n",
"frame = mass.frame()\n",
"data.plotOn(frame, ROOT.RooFit.Name(\"curve_data\"))\n",
"\n",
"legend = ROOT.TLegend(0.4, 0.5, 0.9, 0.9)\n",
"\n",
"print \"{:<20s}{:^4s}{:^12s}{:^12s}{:^17s}{:^17s}\".format(\"\", \"ndof bkg\", \"min nll\", \"chi/ndof\", \"mH\", \"mu\")\n",
"print \"=\" * 80\n",
"\n",
"mu_values, mu_values_err, mH_values, mH_values_err = [], [], [], []\n",
"models_for_fit = [m for m in models.keys() if (m != 'complicated' and m != 'landau')]\n",
"\n",
"for i, (color, pdfname) in enumerate(zip(colors, models_for_fit)):\n",
" mH.setVal(1)\n",
" mH.setVal(125)\n",
" pdf = ws.pdf('model_' + pdfname)\n",
" fit_result = pdf.fitTo(data, ROOT.RooFit.Save())#, ROOT.RooFit.Minos(True))\n",
" \n",
" curve_sb_name = \"curve_sb_%s\" % pdfname\n",
" pdf.plotOn(frame, ROOT.RooFit.LineColor(color), ROOT.RooFit.Name(curve_sb_name))\n",
" pdf.plotOn(frame, ROOT.RooFit.Components(\"bkg*\"), ROOT.RooFit.LineColor(color), ROOT.RooFit.LineStyle(ROOT.kDashed))\n",
" chi2 = frame.chiSquare(curve_sb_name, \"curve_data\", bkgs[pdfname][0])\n",
" print u\"{:<20s}{:>4d}{:>12.0f}{:>12f}{:>8.2f} ± {:<8.2f}{:>5.3f} ± {:<7.2f}\".format(pdfname, bkgs[pdfname][0], fit_result.minNll(), chi2, mH.getVal(), mH.getError(), mu.getVal(), mu.getError())\n",
" curve = frame.findObject(curve_sb_name)\n",
" curve.SetFillColor(0)\n",
" curve.SetMarkerSize(0)\n",
" legend.AddEntry(curve, pdfname)\n",
" \n",
" mu_values.append(mu.getVal())\n",
" mH_values.append(mH.getVal())\n",
" mu_values_err.append(mu.getError())\n",
" mH_values_err.append(mH.getError()) \n",
"\n",
"legend.SetNColumns(2)\n",
"frame.addObject(legend)\n",
"frame.SetTitle(\"s+b fit\")\n",
"frame.Draw()\n",
"\n",
"fig, ax = plt.subplots(figsize=(5, 8))\n",
"ax.errorbar(np.arange(len(models_for_fit)) + 0.5, mu_values, mu_values_err, fmt='.')\n",
"ax.set_ylim(1 - 3 * np.mean(mu_values_err), 1 + 3 * np.mean(mu_values_err))\n",
"ax.set_xticks(np.arange(len(models_for_fit)) + 0.5)\n",
"ax.set_xticklabels(models_for_fit, rotation='vertical')\n",
"ax.hlines(1., 0.5, len(models_for_fit), linestyles='dotted')\n",
"ax.set_ylabel(\"$\\hat{\\mu}$\")\n",
"plt.show()\n",
"\n",
"canvas"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"obs = mu\n",
"nbins = 40"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"mH.setVal(125)\n",
"mu.setVal(1)\n",
"\n",
"mH.setConstant(False)\n",
"mu.setConstant(False)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.999599047823 0.108023172576 125.000141218 0.167403437545\n"
]
}
],
"source": [
"true_model.fitTo(data)\n",
"mu_error = mu.getError()\n",
"print mu.getVal(), mu.getError(), mH.getVal(), mH.getError()\n",
"obs_min, obs_max = 1 - 8 * mu_error, 1 + 8 * mu_error"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"exp\n",
"poly1\n",
"poly2\n",
"poly3\n",
"poly4\n",
"poly5\n",
"exppol2\n",
"exppol3\n",
"exppol4\n",
"laurent2\n",
"laurent3\n",
"laurent_half\n",
"[<ROOT.RooCurve object (\"nll2P_exp\") at 0x5d50f80>, <ROOT.RooCurve object (\"nll2P_poly1\") at 0x5d170a0>, <ROOT.RooCurve object (\"nll2P_poly2\") at 0x5d719d0>, <ROOT.RooCurve object (\"nll2P_poly3\") at 0x5d7b290>, <ROOT.RooCurve object (\"nll2P_poly4\") at 0x5dc5750>, <ROOT.RooCurve object (\"nll2P_poly5\") at 0x5db3670>, <ROOT.RooCurve object (\"nll2P_exppol2\") at 0x5d892a0>, <ROOT.RooCurve object (\"nll2P_exppol3\") at 0x5e0ee80>, <ROOT.RooCurve object (\"nll2P_exppol4\") at 0x5e55110>, <ROOT.RooCurve object (\"nll2P_laurent2\") at 0x5e88dd0>, <ROOT.RooCurve object (\"nll2P_laurent3\") at 0x5ea6860>, <ROOT.RooCurve object (\"nll2P_laurent_half\") at 0x5c77610>]\n"
]
}
],
"source": [
"nbins = 40\n",
"\n",
"canvas_noc = TCanvas()\n",
"\n",
"canvas = TCanvas()\n",
"canvas.Divide(1, 2)\n",
"canvas.cd(1)\n",
"\n",
"frame = obs.frame(ROOT.RooFit.Bins(nbins), ROOT.RooFit.Range(obs_min, obs_max))\n",
"frame_noc = obs.frame(ROOT.RooFit.Bins(nbins), ROOT.RooFit.Range(obs_min, obs_max))\n",
"frame_control = obs.frame(ROOT.RooFit.Bins(nbins), ROOT.RooFit.Range(obs_min, obs_max))\n",
"\n",
"legend = ROOT.TLegend(0.7, 0.7, 0.9, 0.9)\n",
"\n",
"\n",
"for color, pdfname in zip(colors, models_for_fit):\n",
" ndof = bkgs[pdfname][0]\n",
" print pdfname\n",
" mH.setVal(125)\n",
" mu.setVal(1)\n",
"\n",
" mH_set = ROOT.RooArgSet(obs)\n",
" ws.pdf(\"model_\" + pdfname).fitTo(data)\n",
" nll = ws.pdf(\"model_\" + pdfname).createNLL(data)\n",
"\n",
" nll_min_value = nll.getVal()\n",
"\n",
" nllp = nll.createProfile(mH_set)\n",
" \n",
" nllp_list = ROOT.RooArgList(nllp, ROOT.RooFit.RooConst(nll_min_value))\n",
" nllp_list2 = ROOT.RooArgList(nllp, ROOT.RooFit.RooConst(nll_min_value), ROOT.RooFit.RooConst(ndof))\n",
" \n",
"\n",
" nllp2 = ROOT.RooFormulaVar('nllp2', \"-2 log #\\lambda\", '2 * (@0 + @1)', nllp_list)\n",
" nllp2P = ROOT.RooFormulaVar('nllp2P', \"-2 log#\\lambda_P\", '2 * (@0 +@1) + @2', nllp_list2)\n",
" nllp2.plotOn(frame_noc, ROOT.RooFit.LineColor(color), ROOT.RooFit.LineWidth(1), ROOT.RooFit.FillColor(0), ROOT.RooFit.Precision(1), ROOT.RooFit.LineStyle(ROOT.kDotted))\n",
" nllp2P.plotOn(frame, ROOT.RooFit.LineColor(color), ROOT.RooFit.LineWidth(1), ROOT.RooFit.FillColor(0), ROOT.RooFit.Precision(1), ROOT.RooFit.Name(\"nll2P_%s\" % pdfname))\n",
" \n",
"\n",
"curves = []\n",
"for pdfname in models_for_fit:\n",
" curve = frame.findObject(\"nll2P_%s\" % pdfname)\n",
" curves.append(curve)\n",
" curve.SetMarkerSize(0)\n",
" legend.AddEntry(curve, pdfname)\n",
"\n",
"print curves \n",
" \n",
"envelope = ROOT.TGraph()\n",
"envelope_index = ROOT.TGraph()\n",
"for i, x in enumerate(np.linspace(obs_min, obs_max, 50)): \n",
" values = [curve.Eval(x) for curve in curves]\n",
" min_index = np.argmin(values)\n",
" envelope_index.SetPoint(i, x, min_index)\n",
" envelope.SetPoint(i, x, values[min_index])\n",
"envelope.SetLineStyle(ROOT.kDashed)\n",
"legend.AddEntry(envelope, \"envelope\")\n",
"envelope.SetLineWidth(3)\n",
"envelope.SetMarkerSize(0)\n",
"envelope.SetFillStyle(0)\n",
"\n",
"frame.addObject(envelope)\n",
"frame.addObject(legend)\n",
" \n",
"frame.Draw()\n",
"\n",
"canvas.cd(2)\n",
"envelope_index.SetMarkerSize(0)\n",
"envelope_index.Draw(\"APL\")\n",
"\n",
"canvas_noc.cd()\n",
"frame_noc.addObject(legend)\n",
"frame_noc.Draw()"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"False"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ws.writeToFile(\"f.root\")"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"ws.Print()"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"canvas_noc.Draw()"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"<ROOT.TCanvas object (\"c1_n2\") at 0x5a5cb60>"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"canvas"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"64.85185093841042"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ws.var(\"bkg_laurent2_par3\").getVal()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"ws.Print()"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"<ROOT.RooFitResult object at 0x(nil)>"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ws.pdf('model_poly5').fitTo(data)"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"collapsed": false
},
"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.5"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment