Logo ROOT  
Reference Guide
rf801_mcstudy.py
Go to the documentation of this file.
1 ## \file
2 ## \ingroup tutorial_roofit
3 ## \notebook
4 ## Validation and MC studies: toy Monte Carlo study that perform cycles of event generation and fitting
5 ##
6 ## \macro_code
7 ##
8 ## \date February 2018
9 ## \authors Clemens Lange, Wouter Verkerke (C++ version)
10 
11 import ROOT
12 
13 
14 # Create model
15 # -----------------------
16 
17 # Declare observable x
18 x = ROOT.RooRealVar("x", "x", 0, 10)
19 x.setBins(40)
20 
21 # Create two Gaussian PDFs g1(x,mean1,sigma) anf g2(x,mean2,sigma) and
22 # their parameters
23 mean = ROOT.RooRealVar("mean", "mean of gaussians", 5, 0, 10)
24 sigma1 = ROOT.RooRealVar("sigma1", "width of gaussians", 0.5)
25 sigma2 = ROOT.RooRealVar("sigma2", "width of gaussians", 1)
26 
27 sig1 = ROOT.RooGaussian("sig1", "Signal component 1", x, mean, sigma1)
28 sig2 = ROOT.RooGaussian("sig2", "Signal component 2", x, mean, sigma2)
29 
30 # Build Chebychev polynomial pdf
31 a0 = ROOT.RooRealVar("a0", "a0", 0.5, 0., 1.)
32 a1 = ROOT.RooRealVar("a1", "a1", -0.2, -1, 1.)
33 bkg = ROOT.RooChebychev("bkg", "Background", x, ROOT.RooArgList(a0, a1))
34 
35 # Sum the signal components into a composite signal pdf
36 sig1frac = ROOT.RooRealVar(
37  "sig1frac", "fraction of component 1 in signal", 0.8, 0., 1.)
38 sig = ROOT.RooAddPdf(
39  "sig", "Signal", ROOT.RooArgList(sig1, sig2), ROOT.RooArgList(sig1frac))
40 
41 # Sum the composite signal and background
42 nbkg = ROOT.RooRealVar(
43  "nbkg", "number of background events, ", 150, 0, 1000)
44 nsig = ROOT.RooRealVar("nsig", "number of signal events", 150, 0, 1000)
45 model = ROOT.RooAddPdf(
46  "model", "g1+g2+a", ROOT.RooArgList(bkg, sig), ROOT.RooArgList(nbkg, nsig))
47 
48 # Create manager
49 # ---------------------------
50 
51 # Instantiate ROOT.RooMCStudy manager on model with x as observable and given choice of fit options
52 #
53 # The Silence() option kills all messages below the PROGRESS level, only a single message
54 # per sample executed, any error message that occur during fitting
55 #
56 # The Extended() option has two effects:
57 # 1) The extended ML term is included in the likelihood and
58 # 2) A poisson fluctuation is introduced on the number of generated events
59 #
60 # The FitOptions() given here are passed to the fitting stage of each toy experiment.
61 # If Save() is specified, fit result of each experiment is saved by the manager
62 #
63 # A Binned() option is added in self example to bin the data between generation and fitting
64 # to speed up the study at the expemse of some precision
65 
66 mcstudy = ROOT.RooMCStudy(
67  model,
68  ROOT.RooArgSet(x),
69  ROOT.RooFit.Binned(
70  ROOT.kTRUE),
71  ROOT.RooFit.Silence(),
72  ROOT.RooFit.Extended(),
73  ROOT.RooFit.FitOptions(
74  ROOT.RooFit.Save(
75  ROOT.kTRUE),
76  ROOT.RooFit.PrintEvalErrors(0)))
77 
78 # Generate and fit events
79 # ---------------------------------------------
80 
81 # Generate and fit 1000 samples of Poisson(nExpected) events
82 mcstudy.generateAndFit(1000)
83 
84 # Explore results of study
85 # ------------------------------------------------
86 
87 # Make plots of the distributions of mean, error on mean and the pull of
88 # mean
89 frame1 = mcstudy.plotParam(mean, ROOT.RooFit.Bins(40))
90 frame2 = mcstudy.plotError(mean, ROOT.RooFit.Bins(40))
91 frame3 = mcstudy.plotPull(mean, ROOT.RooFit.Bins(
92  40), ROOT.RooFit.FitGauss(ROOT.kTRUE))
93 
94 # Plot distribution of minimized likelihood
95 frame4 = mcstudy.plotNLL(ROOT.RooFit.Bins(40))
96 
97 # Make some histograms from the parameter dataset
98 hh_cor_a0_s1f = ROOT.RooAbsData.createHistogram(
99  mcstudy.fitParDataSet(), "hh", a1, ROOT.RooFit.YVar(sig1frac))
100 hh_cor_a0_a1 = ROOT.RooAbsData.createHistogram(mcstudy.fitParDataSet(),
101  "hh", a0, ROOT.RooFit.YVar(a1))
102 
103 # Access some of the saved fit results from individual toys
104 corrHist000 = mcstudy.fitResult(0).correlationHist("c000")
105 corrHist127 = mcstudy.fitResult(127).correlationHist("c127")
106 corrHist953 = mcstudy.fitResult(953).correlationHist("c953")
107 
108 # Draw all plots on a canvas
109 ROOT.gStyle.SetPalette(1)
110 ROOT.gStyle.SetOptStat(0)
111 c = ROOT.TCanvas("rf801_mcstudy", "rf801_mcstudy", 900, 900)
112 c.Divide(3, 3)
113 c.cd(1)
114 ROOT.gPad.SetLeftMargin(0.15)
115 frame1.GetYaxis().SetTitleOffset(1.4)
116 frame1.Draw()
117 c.cd(2)
118 ROOT.gPad.SetLeftMargin(0.15)
119 frame2.GetYaxis().SetTitleOffset(1.4)
120 frame2.Draw()
121 c.cd(3)
122 ROOT.gPad.SetLeftMargin(0.15)
123 frame3.GetYaxis().SetTitleOffset(1.4)
124 frame3.Draw()
125 c.cd(4)
126 ROOT.gPad.SetLeftMargin(0.15)
127 frame4.GetYaxis().SetTitleOffset(1.4)
128 frame4.Draw()
129 c.cd(5)
130 ROOT.gPad.SetLeftMargin(0.15)
131 hh_cor_a0_s1f.GetYaxis().SetTitleOffset(1.4)
132 hh_cor_a0_s1f.Draw("box")
133 c.cd(6)
134 ROOT.gPad.SetLeftMargin(0.15)
135 hh_cor_a0_a1.GetYaxis().SetTitleOffset(1.4)
136 hh_cor_a0_a1.Draw("box")
137 c.cd(7)
138 ROOT.gPad.SetLeftMargin(0.15)
139 corrHist000.GetYaxis().SetTitleOffset(1.4)
140 corrHist000.Draw("colz")
141 c.cd(8)
142 ROOT.gPad.SetLeftMargin(0.15)
143 corrHist127.GetYaxis().SetTitleOffset(1.4)
144 corrHist127.Draw("colz")
145 c.cd(9)
146 ROOT.gPad.SetLeftMargin(0.15)
147 corrHist953.GetYaxis().SetTitleOffset(1.4)
148 corrHist953.Draw("colz")
149 
150 c.SaveAs("rf801_mcstudy.png")
151 
152 # Make ROOT.RooMCStudy object available on command line after
153 # macro finishes
154 ROOT.gDirectory.Add(mcstudy)